Route to turbulence via oscillatory states in polar active fluid under confinement
Sora Shiratani,1, ∗ Kazumasa A. Takeuchi,1, 2, † and Daiki Nishiguchi1, 3, ‡arXiv:2304.03306v1 [cond-mat.stat-mech] 6 Apr 2023

1

Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
2
Institute for Physics of Intelligence (ipi), The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
3
PRESTO, Japan Science and Technology Agency, 4-1-8 Honcho, Saitama 332-0012, Japan(Dated: April 10, 2023)

We report a novel route to active turbulence, observed in numerical simulations of a polar active
fluid model under confinement. To deal with large-scale computations with arbitrary geometries, wedeveloped a GPU-based scheme that can be used for any boundary shape in a unified manner. For
the circular confinement, as the radius was increased, we found a series of transitions first from asingle stationary vortex to an oscillating pair of vortices, then through reentrant transitions between
oscillatory and chaotic dynamics before finally reaching the active turbulence. The first transition
turned out to be hysteretic, with the emergence of the oscillatory state consistent with the subcritical
Hopf bifurcation. In dumbbell-shaped boundaries composed of two overlapping circles, we observed
a transition comparable to the ferromagnetic-antiferromagnetic vortex-order transition reported inprevious experiments, but the transition point turned out to show a qualitatively different geometry
dependence.

I.

INTRODUCTION

Routes to chaos and turbulence have been one of
the central topics in fluid mechanics at high Reynolds
numbers [1, 2] and statistical physics. However, turbulent phenomena are widely observed beyond the highReynolds-number realm, with a notable example of active
turbulence [3] in low-Reynolds-number active matter systems, for which the route to turbulence remains largely
unexplored. Active matter refers to a collection of selfpropelled particles and it often exhibits collective motion
due to alignment interaction. While ordered collective
motion often arises in theoretical models [4] and also occasionally in experiments [5–7], it is not rare that destabilizing interactions also act and render the collective motion turbulent. Such active turbulence has indeed been
observed in various experimental systems, such as reconstituted cytoskeletal systems [8], electrokinetic Janus
particles [9], sperms [10] and bronchial epithelial cell cultures [11]. As demonstrated thereby, active turbulence ischaracterized by collective motion with many swirls andvortices despite the low Reynolds numbers, which has
been diagnosed, among other approaches, through scaling behavior of the power spectrum [3]. Hydrodynamic
descriptions were also proposed, which successfully reproduced dynamics and statistical properties of bulk active
turbulence [3, 12–14].
Besides these developments on bulk systems, the presence of boundaries and confinement have provided interesting new perspectives. While theoretical approaches
are often difficult because of the a priori unknown boundary condition for hydrodynamic descriptions and heavy∗ sora.shiratani@phys.s.u-tokyo.ac.jp† kat@kaztake.org‡ nishiguchi@noneq.phys.s.u-tokyo.ac.jp

computational costs of agent-based simulations including hydrodynamic effects, experiments have shown that
active turbulence often self-organizes into ordered states
under confinements. For example, an ordered vortex has
been observed in circular geometries and a directed flowin channels, in diverse systems such as bacterial suspensions [15, 16], epithelial cells [17], and reconstituted cytoskeletons [18, 19]. Pillars were reported to rectify bacterial active turbulence by pinning topological defects inthe flow field [20–22]. More elaborated geometries, suchas connected circular chambers, were also studied andreported to show transitions between antiferromagnetic
and ferromagnetic vortex order [23–25].
It is then natural to ask how ordered active flows under confinements are destabilized as the confinement is
weakened and end up in the bulk active turbulence, i.e.,
the route to active turbulence. In the literature, routes toactive turbulence have been studied more often withoutboundaries, typically by changing the activity, and different pathways were proposed for different symmetries [3]:While active nematic fluids tend to undergo transitions
akin to excitable systems [26], some compressible polar active fluids show oscillatory phenomena due to selfadvection of the polar order [27, 28]. In contrast, studies
of confined systems were limited so far to a few cases. For
channels, the transition to active nematic turbulence wasnumerically studied and reported to be in the directed
percolation universality class [29], somewhat similarly toshear-driven transitions to turbulence in Navier-Stokes
fluids [30, 31]. For circular confinements, experiments ofreconstituted nematic cytoskeletons showed that topological defects play crucial roles in dynamics and transitionsof vortices under strong confinements [19]. To our knowledge, such routes to turbulence have not been studied so
far for polar active systems.
In this paper, we determine the route to turbulencein polar active fluid under varying confinements, using a
prototypical model of polar active turbulence known as

2
the Toner-Tu-Swift-Hohenberg (TTSH) equation [3, 12–14] and the boundary conditions identified by previous
experiments on bacterial turbulence [21]. By increasing
the radius of the circular confinement, we found a series of transitions starting from a single stationary vortex, passing through intermediate oscillatory states andfinally reaching the active turbulence. The first transition is from a single stationary vortex to a periodically oscillating pair of vortices and turned out to behysteretic. This is followed by a sequence of transitions across periodic and chaotic oscillations, as well asa quasiperiodic one, which consists of two fundamental
frequencies with an irrational ratio, before reaching thefinal turbulent state. This anomalous route to turbulence is a novel scenario different from those known forthe conventional Navier-Stokes turbulence [1, 2] and for
active nematics in a channel [29]. This finding was madepossible by a GPU-based solver that we developed here,
which efficiently integrates the TTSH equation for arbitrary boundary shapes. GPU implementation allowed us
to achieve high-resolution computation within reasonable
time, which was essential to correctly identify intermediate oscillatory states and transitions. It also helped us tocharacterize detailed properties of the transitions, in particular the Lyapunov exponents, which require massive
extra computations. As our scheme works with arbitrary
boundary configurations in a unified manner, we expect
it to be a useful platform to predict behavior of active
turbulence under confinement, which can also be used fordesigning experimental setups for bacterial turbulence.
This paper is organized as follows. In Sec. II, we describe our calculation scheme, regarding how to fully automate calculations under arbitrary-shaped geometries,
possible artifacts caused by ordinary schemes, and advantages of GPU implementation. Then, the results arepresented in two sections. Section III is for circularconfinements, where we unveil the hysteretic stationaryoscillatory transition as well as subsequent transitions
between oscillatory and chaotic states. Section IV reports simulation results for dumbbell-shaped boundaries,
which consist of two overlapping circles, chosen here as atest case of the TTSH equation with complex boundary
shapes. In this geometry, we found a transition comparable to the ferromagnetic-antiferromagnetic vortex-ordertransition reported by previous experiments [24, 25], but
the transition point turned out to show a qualitatively
different geometry dependence. Section V is devoted to
discussions. Finally, we summarize the results and giveconcluding remarks in Sec. VI.

its dimensionless form [21], it reads:∇ · v = 0,

(1)

2
∂v
+ λv · ∇v = av − bv2 v − 1 + ∇2 v − ∇p,∂t(2)

where v(x, t) is a coarse-grained velocity field of bacteria, p is an effective pressure that ensures the incompressibility, and λ, a, b are model parameters that do
not depend on x, t, v or p. Here, following Ref. [21],we express the dimensionless coordinates x and time t in
the unit of the characteristic length and time scales, respectively, determined by the Swift-Hohenberg-like term
2
1 + ∇2 v, which has coefficients otherwise. As a result,the characteristic wavenumber |k∗ | excited by the SwiftHohenberg-like term is equal to unity, so that the typicalsize of vortices is 2π in real space. This corresponds to≈ 105 µm for bacterial turbulence of Bacillus subtilis [21].
In the present work, we focus on two-dimensional systems, to make our results comparable with experimental
observations.
The boundary conditions for the TTSH equation were
determined experimentally in Ref. [21], to be v ≡(vx , vy ) = 0 and the vorticity ω ≡ ∂x vy − ∂y vx = 0,
at least for the experimental condition employed therein.
Then, the existence of boundaries (as well as obstacles)
can be incorporated by adding damping terms into the
TTSH equation in the form of the vorticity equation [21]:∇ · v = 0,∂ω
+ λv · ∇ω = aω − b∇ × v2 v∂t
2
− 1 + ∇2 ω − γv ∇ × [K(x)v] − γω K(x)ω,

A.

(4)

where K(x) is a nonnegative scalar field such that
K(x) ≈ 0 inside the system and K(x) ≈ 1 outside, andγv , γω are positive parameters representing the damping
strengths. The first damping term in Eq. (4) amounts
to adding the damping −γv K(x)v to Eq. (2) to preventv from growing in the area outside the system (K ≈ 1,
referred to as masked area) without affecting the area
inside the system (K ≈ 0, unmasked area). The second damping term −γω K(x)ω does the same for thevorticity. In the following, we set (a, b, λ, γv , γω ) =(0.5, 1.6, 9, 40, 4), which had been reported to quantitatively reproduce the experimental results in Ref. [21](see Table I in Appendix B for the list of the parametervalues used in this paper).
B.

II.

(3)

computation

METHODmodel

We use the TTSH equation, which describes bulk
behavior of turbulent collective motion spontaneously
formed in dense bacterial suspensions [3, 12–14, 21]. In

Although we use the standard pseudospectral method
to integrate the TTSH equation, the calculation is notan easy task because large-scale computations are required by nature. There are two specific reasons for thedifficulty. First, as previously stated, we use the virtual boundaries represented by the damping terms. This

3
(b)

(a)

(c)
pre
post

K(0, y)

1.00
0.75
0.50
0.25
0.00
−19.00−18.75−18.50y−18.25−18.00

FIG. 1: Emoji boundary. (a) The mask K(x) generated by our automatic mask generator. (b) A snapshot of thevorticity field from a simulation using the generated mask. Note that the color scale displayed here is used in all
vorticity snapshots in the present paper. See also Movie S1 [32]. (c) Emoji mask before and after the low-passfiltering in the automatic mask generator, sliced at x = 0.

method is sometimes problematic because v and ω penetrate into the damped area outside the system. To avoid
artifacts, damping walls need to be sufficiently thick sothat the penetrating fields decay to almost zero. Second,
the spatial discretization interval ∆x needs to be small
because the results may be influenced by small changes
in the mask field K(x) as addressed in the next section.
Fortunately, this model is not so sensitive to the temporal
discretization interval ∆t if implemented properly. This
allowed us to use the Euler method and obtain converging
results efficiently. Detailed descriptions of the algorithm
are provided in Appendix C. Care was taken to makefull use of the performance of GPU (see Appendix E).
Taking these properties into account, we discretized themodel on an 8192 × 8192 square lattice with ∆x = 0.005
and simulated its time evolution with ∆t = 0.01, unlessotherwise stated. Note that most of the calculation was
done for the masked region (typically 90 % in area), toconstruct a thick barricade to ensure the decay of thepenetrating fields.

C.

mask generation

Although the physical meaning of the mask field K(x)is straightforward, due care should be taken to design
a suitable K(x) for computation, due to constraints
originating from the pseudospectral method (see Appendix D). One such constraint is the 12 rule for antialiasing in the case of the cubic nonlinearity (see Appendix C). In short, K(x) should not contain highwavenumber modes in order to properly integrate thenonlinear terms of the TTSH equation, and therefore
a naive step-function-like binary mask cannot be used.
An intuitive treatment is to interpolate the gap between
K = 0 (unmasked) and K = 1 (masked) by a slowlyvarying function such as tanh, but this is not an ideal
solution because tanh contains high-wavenumber modes
even though their amplitudes are significantly smaller

than the original binary mask. Moreover, this approachrequires us to design and encode the mask explicitly as apart of the program before the calculation. This task is
not straightforward when the geometry is complex, andthere is no general way of predicting whether aliasing
noise associated with the mask is actually tolerable ornot. Therefore, even though the TTSH equation has been
actually calculated by using such a mask without breaking down, possibly thanks to the Swift-Hohenberg-like2term 1 + ∇2 v which damps high-wavenumber modes,
it is desirable to design and use a mask that is free of any
aliasing noise.
To overcome these difficulties, we adopted a completely
different approach, which automatically generates a mask
K(x) that strictly satisfies all the requirements by using
a low-pass filter (see Appendix D for details). As a result,the only required input is a binary scalar field indicatingwhere to mask. Figure 1 displays a toy example usinga binarized unicode emoji [33] as the boundary (see also
Movie S1 [32]). As a result of the low-pass filtering, thegenerated mask K(x) oscillates as displayed in Fig. 1(c).
Although this side effect is more or less inevitable, theoscillation can be suppressed by introducing additional
low-pass filtering operations, at the price of sharpness.
We, however, did not take this option intentionally because we prioritized the sharpness and we did not see any
discernible artifact.

III.

CONFINEMENT IN CIRCLE

A.

motivations and background

Circle is one of the simplest geometries parametrized
only by its radius R, and its simplicity has long been
arousing experimental interests as surveyed in the introduction section. In particular, it was reported that bacterial flow can be stabilized into a single vortex if confined

4in a sufficiently small circular geometry[15, 24, 25], while
the opposite, unconfined limit corresponds to turbulence.
This led us to perform numerical simulations in circularareas for various radii R and investigated the route to active turbulence in this case. Unless otherwise stated, we
started from a random initial state and discarded transients to ensure that the system is in a steady state. Notethat high-resolution calculation realized by GPU is crucial here (see Appendix E), because it turned out that
small changes in R may affect qualitative features of the
flow.
B.

overview of results

First we outline the results obtained in this geometry. Our calculations revealed a rich phase diagram withdifferent numbers of vortices [Fig. 2(a)] and different dynamical states [Fig. 2(b) and Movies S2-S8 [32]] depending on the radius R. While no vortex is observed forR . 5.1 due to strong damping, around R = 5.2 we observe a single stationary vortex generated at the centerof the circle (blue box in Fig. 2 and Movie S3 [32]). As
R is increased, the number of vortices increases to two(orange box), four (green), and more (purple). Dynamics
becomes non-stationary as soon as multiple vortices are
generated. In particular, the first non-stationary state
observed around R = 5.4 consists of an oscillating pairof vortices, as shown in Fig. 3 and Movie S4 [32]. The
transition between the single-vortex stationary state andthe vortex-pair oscillatory state will be characterized indetail, in Sec. III D.
To characterize the changes in the dynamics morequantitatively, we measure the temporal power spectral
density (PSD) of the vorticity field (Fig. 4), defined by
2

1 X X
S(Ω) ≡ 2
ω(x, t)e−iΩt .
N x
t

(5)

At R = 5.4 (Movie S4 [32]), we confirm the periodic oscillation by the PSD that consists of a single fundamental
frequency and its harmonics. At R = 5.6 (Movie S5 [32]),
the first chaotic region suddenly appears while the oscillatory dynamics is still observed. At R = 5.8 (Movie S6[32]), the system regains regularity but another fundamental frequency emerges and the oscillation becomesquasiperiodic (with the fundamental frequencies being
0.01125 and 0.0175). Simultaneously, the system now
accommodates four vortices (see Fig. 2, Movie S2 [32]).
From R = 6.0 to R = 7.4, the system goes back and forth
between the chaotic and oscillatory states (see Movie S7[32] for R = 7.2) until it finally falls into the chaotic state
at R = 7.6 (Movie S8 [32]) and never returns. FromR = 7.6, the system has more than four vortices and thenumber increases with R.
In the following, we will characterize the single-vortex
stationary state (Sec. III C) and the transition to thevortex-pair oscillatory state (Sec. III D). Then we willcomment on the route to turbulence observed in this geometry (Sec. III E).
C.

single-vortex stationary state

In this state observed around R = 5.2 (blue box in
Fig. 2), a single vortex is generated and located at thecenter of the circle. The sign of the vorticity is determined by the initial condition that we generated at random. This is reminiscent of the single-vortex stationary
state reported in experiments [15, 24, 25, 34], even though
the boundary condition may differ and no edge currents
are observed in our simulations.
It is useful to compare the velocity and vorticity fields
in our single-vortex stationary state with those of an analytical solution to the linearized version of the TTSH
equation [21]. As described in Appendix A, the general stationary solution to the TTSH equation without
nonlinear terms, expressed in terms of polar coordinates(r, θ), is given by∞ X
X

Cn± ωn± (r, θ),(6)

ωn± (r, θ) ≡ Jn (k± r) cos Θn ,(7)

ω(x) =

n=0 ±

andv(x) =∞ X
X

Cn± vn± (r, θ),

(8)

n=0 ±

Jn0 (k± r)
nJn (k± r)
sin Θn cos θ,
cos Θn sin θ −
2r
k±
k±
(9)
0
(k
r)
J
nJ
(k
r)
±
n ±
±
vn,y
(r, θ) = n
sin Θn sin θ,cos Θn cos θ +
2r
k±
k±
(10)
±
vn,x
(r, θ) =±±being the x and y components, reand vn,ywith vn,x±spectively,ofv,J
the Bessel function of the first kind,np√nk± ≡ 1 ± a, and Θn ≡ nθ + const. Since the numerically observed single-vortex state is isotropic (i.e.,
independent of θ), we are led to compare with the analytic solutions with n = 0. The result is displayed inFig. 5(a). Interestingly, we found that both the tangential component of the velocity field vθ (r, θ) and the vorticity field ω(r, θ) have qualitative features in common
with a linear combination of those analytic solutions, especially near the center, despite the existence of the nonlinear terms in our simulations. This may partly be because the advection term λv · ∇ω of Eq. (4) vanishes forthe isotropic solution. However, since the analytic solutions cannot satisfy the boundary conditions for v and ω
simultaneously, they cannot describe the numerical observation precisely. We may argue that the existence ofthe cubic term −bv2 v of the TTSH equation may serve,
effectively, for the Bessel-type solution to adjust itself toreconcile with the required boundary conditions.

5

(a)

(b)

FIG. 2: Phase diagram for the circular confinement, showing snapshots of the vorticity field for different R (samecolor scale as Fig. 1(b)). The frames of the snapshots are given different colors corresponding to the number ofvortices (blue: one, orange: two, green: four, purple: more). (a) Typical snapshot chosen for each R. (b) Time seriesof the vorticity field. See also Movie S2 [32].

6

1040

R
5.4

1035

5.6

1030

5.8
6.0

S [a.u.]

1025

6.2
1020

6.4
6.6

1015

6.8
7.0

1010

7.2
105
FIG. 3: Time series of the vorticity field in theoscillating vortex pair state at R = 5.4. Same color
scale as Fig. 1(b). See also Movie S4 [32].

D.

7.4
7.6

100
0.0

0.1

0.2

transition to vortex-pair oscillatory state

As R is increased from the single-vortex stationary
state, the system undergoes a transition to the vortexpair oscillatory state (orange box in Fig. 2). This state
consists of two oscillating vortices with vorticities of different signs, as already shown in Figs. 2(b) and 3, as wellas in Movies S2 and S4 [32].
The transition between the single-vortex stationary
state and the vortex-pair oscillatory state turned out tobe hysteretic, as illustrated in Fig. 6(a). A useful quantity to capture this transition is the vortex order parameter [15]P|eθ (x) · v(x)|
1
2
xP
Ψ=
−
,
(11)
π
1 − π2x |v(x)|where eθ (x) is the azimuthal unit vector and the sum
is taken over the entire space. By construction, Ψ =1 if the velocity field is completely azimuthal, whereasΨ = 0 if it is completely disordered. To investigate thehysteresis, we ran simulations sequentially, as follows. To
go up, starting from the steady state at R = 5.20, we
increased R by ∆R = 0.02 and measured Ψ typically over
5-10 periods after the system reached the steady state.
We repeated this step until R = 5.34. To go down, we
did the same in the opposite direction. The result isshown in Fig. 6(b), which clearly shows the hysteresis.
From this, we approximately estimated the lower andupper transition points at Rcdown ' 5.23 and Rcup ' 5.31,
respectively.

0.3

0.4

0.5

Ω
2π
FIG. 4: Temporal PSD of the vorticity field [Eq. (5)]measured for different R. Each PSD is multiplied by aconstant to avoid overlapping.
(a)

(b)

0.6

R
0.4

ω
ω fit
vθ
vθfit

0.6
0.3

0.2

0.0

0.0
0

ω(k+)
ω(k−)
vθ (k+)vθ (k−)

0.9

1

2

3

r

4

5

6

−0.3

0

1

2

3

r

4

5

6

FIG. 5: Comparison between the single-vortex
stationary state and analytic solutions for the linearized
TTSH equation. (a) Radial dependence of vθ (blue) andω (orange) for the single-vortex stationary state at
R = 5.2 (solid) and those for the best-fit linear
combination of the analytical solutions with n = 0(dashed). The fitting was carried out based on both v
and ω normalized by their mean absolute values. The
same set of C0± is used for vθ and ω. (b) Two±
independent modes of the analytical solutions, v0,θand±
ω0 for k = k± .
To quantitatively characterize the dynamic aspect of
this transition, we measured the Lyapunov exponents,
i.e., the exponential growth rates of infinitesimal pertur-

7
(b)

(a)

Rcup

−0.010−0.015−0.020−0.025λ1
λ2

5.20 5.22 5.24 5.26 5.28 5.30 5.32 5.34

R

−0.1
−0.2
λ1
λ2
λ3
λ4

−0.3

−0.4
5.20 5.22 5.24 5.26 5.28 5.30 5.32 5.34

R

0.0

800

−0.1

600

103

T

0.0

T

−0.005(f)(e)
Lyapunov exponents

Lyapunov exponents

0.000

Lyapunov exponents

(d)

(c)

−0.2
λ1
λ2
λ3
λ4

−0.4
5.20 5.22 5.24 5.26 5.28 5.30 5.32 5.34

R

200

1
3
1
2
10−3

400

−0.3

102.5

10−2

R − Rcfit

Rcfit

5.23 5.24 5.25 5.26 5.27 5.28 5.29

R

FIG. 6: Hysteretic transition between the single-vortex stationary state and the vortex-pair oscillatory state. (a)Snapshots of the vorticity field (same color scale as Fig. 1(b)) as R is increased (top) or decreased (bottom). (b)Vortex order parameter Ψ. The time-averaged values are shown. The shades indicate the temporal standard
deviations. (c)(d) Lyapunov exponents for increasing R. The two and four largest exponents are shown in (c) and(d), respectively. The first two exponents increase linearly with R and cross zero at R = Rcup , consistently with thecase of the subcritical Hopf bifurcation. (e) Four largest Lyapunov exponents for decreasing R. (f) Oscillation
period T above the lower transition point Rcdown . Its standard deviation is indicated by error bars, visible only for
R = 5.225 (open symbol) where the oscillation was fluctuating. Blue solid lines indicate results of fitting by−p
T ∝ R − Rcfit , for several choices of data points to use (hence many blue lines are drawn). Inset: The same datain the log-log scale. Two dashed lines are guides for the eyes corresponding to power laws with exponent 1/2 (green)and 1/3 (orange).

bations to the solution, which can be a direct clue todetermine the type of bifurcation underlying this transition [35]. To calculate them, following the standard
method [36, 37], we simulated the time evolution of independent perturbations δωi (x, t) (i: index) along withω(x, t) and measured the exponential growth rates using the QR decomposition. The Lyapunov exponents λiwere then obtained in ascending order. This procedure
is even more costly than the main calculation and tookdays even though we had access to cutting-edge GPUs.
The result is shown in Fig. 6(c) for increasing R, with
∆R = 0.01. This shows that the two largest exponents
λ1 , λ2 , which are negative for R < Rcup as expected, increase linearly with R and reach zero at R = Rcup . Thisbehavior is consistent with the subcritical Hopf bifurcation, which indeed shows a hysteretic transition to an
oscillatory state [35]. It is also notable that the oscillatory state involves two vanishing exponents, despite theabsence of quasiperiodic behavior. We consider that therotational symmetry of the system introduces the second
vanishing exponent, in addition to the one corresponding
to the time translation symmetry. Regarding the thirdand fourth exponents λ3 , λ4 , they also increase linearlywith R up to Rcup [Fig. 6(d)]. Interestingly, extrapolationof this linear dependence suggests that it would cross zeroat R ' 5.65, which is close to the transition point to thestate with four oscillating vortices.
For decreasing R, the result is shown in Fig. 6(e).
Given that the transition at Rcup was consistent withthe subcritical Hopf bifurcation, the standard scenario
from the low-dimensional dynamical systems theory suggests that the lower transition may be described by thesaddle-node bifurcation, in which case the first negative
1/2
exponent increases as ∝ R − Rcdown
when approaching Rcdown from above [35]. However, it was not our case:
The first negative exponent λ3 does not approach zerobut remains at ' −0.2 as displayed in Fig. 6(e). We
also measured the oscillation period T as a function of R
[Fig. 6(f)], using the half smaller discretization interval∆x to reduce discretization effect (with the edge length
N ∆x kept unchanged; see Table I in Appendix B). If thetransition were the saddle-node bifurcation, we would ex−1/2
pect T ∝ R − Rcdown. However, while our data indeed show seemingly power-law divergence near Rcdown ,8FIG. 7: Time series of the vorticity field for R ≥ 8 (same color scale as Fig. 1(b)). Note that the results presented
here are affected by relatively strong effects of the periodic boundary, because R is close to L/2 and the vorticity
field does not decay sufficiently at the periodic boundary.
−pi.e., T ∝ R − Rcfit , we estimated the exponent p at
p = 0.35(1). This was obtained by varying the rangeof fitting, on which the result hardly depends, exceptthat we obtained p = 0.39(1) if the point closest totransition, at which the oscillation seemed less stable,
was excluded. This suggests that the bifurcation we observed may not be understood within the framework of
the low-dimensional dynamical systems theory. In otherwords, this hints at a hitherto unknown bifurcation in
high-dimensional dynamical systems.

E.

route to turbulence

Finally we comment on the route to turbulence observed in this circular confinement. As already describedin Sec. III B and Fig. 2, the first nontrivial state observed in this geometry is the single-vortex stationary
state (R = 5.2). As R is increased, the system first undergoes a hysteretic transition to the vortex-pair oscillatory
state (5.23 . R . 5.31). This state is periodic for R close
to the transition, while chaotic modulation may be added
for larger R. From R = 5.8, the system has four vortices
and shows multiple transitions among different dynamical states (periodic, quasiperiodic, and chaotic states) ina reentrant manner. The number of vortices increases
further from R = 7.6 and the system now stays in the
chaotic state (Fig. 7). This state continuously shifts toactive turbulence in the bulk limit (R → ∞).

When compared to the route to turbulence for the
Navier-Stokes turbulence, while the emergence of periodic and quasiperiodic states is also seen in the RuelleTakens-Newhouse scenario [2, 38], the rest of the observations do not correspond to any well-known scenario.
Instead, we note that similar reentrant transitions among
periodic, quasiperiodic, and chaotic states were observedin the Kuramoto-Sivashinsky equation [39], as well as in anumerical simulation of the Navier-Stokes equation whena so-called high-symmetry condition was imposed on theflow field [40]. Although we did not impose such a condition, our highly symmetric circular geometry may be
relevant to the reentrant behavior that characterizes theobserved route to turbulence.

IV.

CONFINEMENT IN DUMBBELL
A.

motivations and background

Here we turn our eyes to a dumbbell-shaped confinement as shown in Fig. 8(a), which consists of two overlapping circles of radius R with the centers separated
by distance D. This test case serves as a benchmark todemonstrate the ability of our method to implement acomplex geometry without handcrafting the mask K(x).
In the case of bacterial turbulence enclosed in a PDMS
device [24, 25], it was reported that a pair of vortices
was formed in dumbbell-shaped boundaries, each vortic-

9
(a)

B.

results

(b)
1.00

3.4
3.5
3.6
3.7

0.75

Ψ2

R
D

R

0.50
0.25
0.00

2.8

3.0

3.2

3.4

3.6

D

(c)

FIG. 8: Confinement in dumbbell. (a) Outline of thedumbbell-shaped boundary made of two overlapping
circles and the associated parameters. (b) Orderparameter Ψ2 vs the center-to-center distance D fordifferent radii R. By construction, Ψ2 = 1 for the
ferromagnetic flow and Ψ2 = 0 for the antiferromagnetic
one. (c) Typical vorticity snapshots taken after
relaxation to the stationary state, for different radii R
and center-to-center distances D (same color scale as
Fig. 1(b)). The background color indicates the value ofΨ2 according to the color bar. For simplicity, somesnapshots are displayed after mirror reflection so thatthe vorticity is always positive in the left half and thatthe stronger vortex has the positive vorticity (the lattercan be realized by applying y 7→ −y).

ity peak located near the center√ of the circle. With athreshold distance Dc (R) = 2R, the signs of the two
vorticity peaks were identical (ferromagnetic vortex order) for D < Dc and opposite (antiferromagnetic vortex
order) for D > Dc [24]. It is therefore of interest to see
if a similar transition takes place in our hydrodynamic
setup.
Our simulations on the dumbbell confinement were
performed as follows. Similarly to the bacterial experiments [24, 25], we set the cavity radius near the characteristic size of a single vortex, R ≈ π. Note that, although even a single vortex could not appear at suchsmall radii R in the case of the circular confinement (seeSec. III B), in the dumbbell confinement we do observe
vortices, presumably because of the relatively larger total
area of the confinement. In the simulations reported below, we started from a random initial state and discarded
transients to ensure that the system is in a steady state,for each choice of R and D.

Figure 8(c) displays our simulation results. At each
fixed R, for small D we observed a single vortex formednear the center of the dumbbell, and for large D a pairof vortices is formed near the centers of the two cavities.
In either case, the vortices are stationary in the steadystate. The latter, vortex-pair state corresponds to theantiferromagnetic state observed in previous experiments[24, 25]. Regarding the former, single-vortex state, it isanalogous to the ferromagnetic state in the sense thatthe vorticity has the same sign in the entire region of thedumbbell, but we did not observe a split of vortices asin the experiments. This difference may be attributed
to the presence of the edge flow under the conditions of
Beppu et al.’s experiments [24, 25], while our choice of theboundary conditions is based on another experimental
setup by Nishiguchi and coworkers [20, 21].
A difference is also seen in the transition between thetwo states. In our simulations, as shown in Fig. 8(c), thethreshold distance Dc (R) decreases√ with increasing R,
while it was increasing as Dc (R) = 2R in Beppu et al.’s
experiments [24]. To quantify this observation, we measured an order parameter specialized
 for this case, defined∗∗by Ψ2 = ω(xleft ) + ω(xright ) / |ω|(x∗left ) + |ω|(x∗right )with x∗left/right = argmaxleft/right |ω|. Here, the subscript“left” and “right” stands for the area within D/2 fromthe center of either cavity. By construction, Ψ2 = 1 forthe ferromagnetic flow and Ψ2 = 0 for the antiferromagnetic one. The results in Fig. 8(b) clearly demonstratethat Ψ2 transitions from ≈ 1 to ≈ 0 as D is increased,with Dc (R) decreasing with increasing R. Therefore,
our Dc (R) is ruled by a law different from the scaling Dc (R) ∝ R observed in Beppu et al.’s experiments[24, 25]. √Beppu et al. [24] accounted for the relation
Dc (R) = 2R on the basis of edge currents, i.e., tangential flow on the boundary, which is incompatible withthe boundary conditions chosen here (see discussions in
Sec. V A). It is therefore reasonable to consider that theabsence of edge current in our setup may be responsible
for this difference.

V.
A.

DISCUSSIONS
edge current

Here we discuss the edge current, i.e., tangential flow
along the boundary, typically seen in experiments undercircular confinements [15, 24, 25, 34]. The edge current
involves non-vanishing tangential velocity at the boundary, vθ 6= 0 in the case of the circular confinement. In
the aforementioned experiments, bacteria actually swim
along the boundary. The formation of counter-rotating
double layers has also been reported in droplet suspensions of B. subtilis [15], while no such counter-rotating
layers were reported in the case of E. coli confined in a

10
microfluidic device [24, 25]. In contrast, in our simulations, vθ continuously decays to zero due to the damping
terms, without edge current nor counter-rotating layer.
This difference clearly results from our choice of theboundary conditions, v = 0 and ω = 0, which werededuced from the experimental observation of dense B.
subtilis suspensions placed on a substrate with microfabricated pillars and bordered by a liquid-air interface[20, 21]. Although counter-rotating layers were not observed in the low-magnification microscopy carried out inthese experiments, the possible existence of such layers
may have led to our boundary conditions, as discussed inRef. [21]. Therefore, we expect that our numerical simulations provide predictions for experimental conditions
similar to those in Refs. [20, 21], rather than those in
the existing experiments of circularly confined bacteria[15, 24, 25, 34].
To extend the model to deal with the case with edge
currents, we may (1) introduce a slip velocity to allow anonvanishing vθ at r = R or (2) use another governing
equation that does not suppress high-wavenumber variations, thus allowing the existence of current near the
no-slip boundary. We tested the approach (1) by using
a damping scheme that only removes the radial velocity component vr . However, we were unable to carryout physically sound simulations in this case, because
the tangential component vθ penetrated deep inside thedamped zone. In passing, we note that this boundary
condition allowed us to set the circle radius R smaller
than the minimum value reported in Sec. III to generatea vortex. This may be related
H for the noR to the fact that,slip condition, the identity r<R ωdxdy = r=R v · dx = 0
guarantees that positive and negative vorticities must exist in the same amount inside the confinement, while forthe slip case ω can escape from the confinement despite
the damping. Regarding the approach (2), we may remove the ∇4 v term of the TTSH equation. The equation then reduces to the incompressible version of the
Toner-Tu equation, for which vortices were reported toappear in bulk systems [41] unlike the original, compressible Toner-Tu equation [42–44]. These approaches may
deepen our understanding of effect of edge current on
structures and dynamics of active fluids in confined geometries.

B.

implementation of the boundary conditions

One may wonder if the boundary conditions, whether
slip or no-slip, can be imposed more directly, without resorting to the damping scheme. In this case, we have todeal with arbitrary-shaped boundaries directly. In theliterature, such a method has been pursued in broader
contexts. For example, in computational fluid dynamics [45], an arbitrary-shaped boundary is typically realized by a tailored mesh, and the model is discretized
and integrated on it. For the TTSH equation, however,the discretization of the fourth-order derivative and theboundary condition of vanishing vorticity is not straightforward on such a mesh. More specifically, one may want
to use the second-order derivative of v (to use the computationally efficient Laplacian) on the boundary, but thisis incompatible with the discrete fourth-order derivativeof the time evolution equation. More sophisticated algorithms may be devised, but we suspect that it is difficult
to avoid uncontrollable approximations and that such algorithms are more costly than the pseudospectral method
we adopted.

VI.

CONCLUDING REMARKSIn this work, we realized numerical simulations of the
TTSH equation with arbitrary-shaped boundaries, andpresented results for the two representative test cases,
namely the circular and dumbbell confinements.
From the computational perspective, first we emphasize that GPU implementation makes calculations fast
and affordable, without the need to use supercomputersor to wait for weeks. Indeed, this paper contains several
simulations that would take months or even years using workstation-class CPU (typically 101 ∼ 102 threads,
3 GHz).
Let us conclude this paper by discussing physical implications of the work. Our numerical investigation with the
TTSH equation has succeeded in reproducing the emergent vortex order reported in the experiments, at least
qualitatively. At the same time it has highlighted thatthe slight difference in the confinement geometry can giverise to quantitatively different vortex structures. Specifically, we have explored a novel route to chaos and turbulence under the circular confinement in the TTSH equation. This route starts with a hysteretic transition froman ordered vortex to an oscillating pair of vortices, consistent with the subcritical Hopf bifurcation. This is followed by reentrant transitions across periodic, quasiperiodic, and chaotic oscillations, until the system finallyreaches the active turbulent state. Since our boundary
conditions are the ones inferred from the specific experimental realization reported in Ref. [21], it is an important future task to investigate how robust our findings
are, for other models and boundary conditions that may
be more suitable for other experimental setups. In thiscontext, it is interesting to note that a Hopf bifurcation
was also reported for another polar active fluid model
without confinement but variable activity [28], suggesting some extent of universality in the route to polar active
turbulence via oscillatory states. We anticipate that our
results will contribute to the fundamental understanding of how turbulent structure develops in active mattersystems.

11
Appendix A: analytic solution for the linearized
TTSH equation

Following Ref. [21], here we describe analytic solutions for the linearized TTSH equation, which amounts
to setting (b, λ) = (0, 0) in Eq. (2). Note that the
dropped nonlinearity cannot be regarded as a perturbation. Therefore, the analytic solutions described inthis section may not necessarily represent the numerical observations even approximately. Nevertheless, theseanalytic solutions help interpret some numerical resultspresented in Sec. III.
The linearized TTSH equation for the vorticity, withthe stationarity condition, reads:
2
0 = aω − 1 + ∇2 ω.

(A1)This can be easily solved in the polar coordinate system(r, θ), on the basis of the real-valued eigenfunctions of
the Laplacian, Jn (kr) cos (nθ + const.), where Jn is theBessel function of the first kind, k > 0, and n = 0, 1, . . . .
The corresponding eigenvalue is −k 2 . Therefore, a solution to Eq.p(A1) needs to satisfy 0 = a − (1 − k 2 )2 , i.e.,√k = k± ≡ 1 ± a. The general solution is then given
by linear combinations of them:∞ X
X

Cn± ωn± (r, θ),(A2)

ωn± (r, θ) ≡ Jn (k± r) cos Θn ,(A3)

ω(x) =

n=0 ±with Θn ≡ nθ + const.
Moreover, in the case of the two-dimensional incompressible fluid studied in this work, the velocity field v(x)can be reconstructed from the vorticity field ω(x) by using the stream function. First we obtain the stream function ψ(x) by solving ∇2 ψ = −ω, then it follows thatvx = ∂y ψ and vy = −∂x ψ. This gives
v(x) =∞ X
X

Cn± vn± (r, θ),

(A4)

n=0 ±

Jn0 (k± r)
nJn (k± r)
cos Θn sin θ −sin Θn cos θ,
2r
k±
k±
(A5)
0
J (k± r)nJn (k± r)±vn,y(r, θ) = n
cos Θn cos θ +sin Θn sin θ,
2r
k±
k±
(A6)

±
vn,x
(r, θ) =±±using the same Cn± as Eq. (A2). Here, vn,xand vn,yare±the x and y components, respectively, of vn .
In relation to our results for the circular domain withfinite radius R (Sec. III), we first note that, if any special
solutions ωn± (r, θ) and vn± (r, θ) satisfy ωn± (R, θ) = 0
and vn± (R, θ) = 0 simultaneously, this can be regarded
as a solution for the circular domain case, aside fromthe ignored nonlinearity. However, these two conditions
are actually never satisfied simultaneously, because forall x > 0 the roots of Jn (x) and Jn0 (x) do not coincide.
We may make ω(x) and v(x) closer to zero at r = R bycombining the two modes with k = k± , but the boundary conditions cannot be exactly satisfied. This point iscommented in Sec. III C.
Appendix B: notation

Table I summarizes the parameters and symbols usedin the paper.
TABLE I: Table of parameters and symbols.
symboldescription
λ=9advection strengtha = 0.5
activity parameterb = 1.6
nonlinearity parameter
γv = 40
velocity damping strengthγω = 4
vorticity damping strength
N = 8192a
number of grid points per dimension∆x = 0.005a
spatial discretization interval∆t = 0.01
temporal discretization interval
L = N ∆xedge length of the calculation areaˆ∗, DFT[∗]
discrete Fourier transform of ∗iDFT[∗]
inverse discrete Fourier transform of ∗

 π
π 2
, ∆x
DFT wavenumberk ∈ − ∆x∆k = 2π
discretization interval of kLK(x)preprocessed damping maska(N, ∆x) = (16384, 0.0025) were used for the oscillation
period T shown in Fig. 6(f).

12Appendix C: calculation scheme

Algorithm 1 is our integration scheme for the TTSH
equation based on the pseudospectral method and theEuler method.
Algorithm 1 integration scheme
Input: (v(t), ω̂(t))
1: rhs1 = −bv2 v − γv Kv
2: rhs2 = −λv · iDFT[ikω̂] − γω Kω
ˆ = ik × rhs
ˆ2ˆ 1 + rhs
3: rhshn2 o iˆ
4: ω̂naive = exp
a − 1 − k2∆t ω̂ + rhs∆t
5: for all k do
6:
if k = 0 then7:ω̂new (0) = 0
P
1(v + rhs1 ∆t)
8:v̂new (0) = exp [(a − 1)∆t] · x √N
2

2
1
π 2
9:
else if k > 2 · ∆x then10:ω̂new (k) = 0
11:v̂new (k) = 0
12:
else
13:
ω̂new (k) = 12 {ω̂naive (k) + ω̂naive (−k)∗ }(k)14:
ψ̂(k) = ω̂new
k2
15:v̂new (k)|x = iky ψ̂(k), v̂new (k)|y = −ikx ψ̂(k)16:
end if
17: end for
Output: (vnew , ω̂new ) = (v(t + ∆t), ω̂(t + ∆t))

The key operations are the discrete Fourier transP1−ik·x
form (DFT) defined by fˆ(k) = √N2
x f (x)e
and its inverse (iDFT) similarly defined by f (x) =
P ˆ
ik·x
√1
. Thanks to them, we can replace allk f (k)eN2the ∇ by −ik. Note that the periodic boundary condition is implicitly assumed, i.e., f (x, y) = f (x ± L, y) =f (x, y ± L), whose artifacts are negligible as long as thewalls are thick enough so that v(x) and ω(x) are sufficiently small at the periodic boundary.
As the first step, GPU calculates all the nonlinear
terms using element-wise operations and DFT, which areˆ A pitfall here is aliasing, i.e., artifacts
then stored to rhs.
caused by the expansion of the range of k due to multiplication in real space. To deal with it, we need toeliminate high-wavenumber modes before and after thisstep, though the “before” operation can actually be optimized out. For the TTSH equation containing the cubicπterm bv2 v, the appropriate cutoff is 12 · ∆x(the 12 cutoffrule).
Next, GPU computes the time evolution of ω̂ by the
Euler method and stores the result to ω̂naive . The replacement ∇ → −ik plays a crucial role here, because it
allows us to integrate the linear terms of the time evolution equation separately and exactly by using the exponential multiplier. Without this technique, calculation
will diverge unless we use ridiculously small ∆t ∼ ∆x4 .
As the last step, GPU performs post-calculation operations. For k = 0, ω̂ is corrected to the exact value 0 andv̂ is directly calculated from the time evolution equationof v, because this spatially uniform component cannot

be recovered from ω. Specifically, we used the following
equation, obtained by taking the spatial average of bothsides of the TTSH equation [Eq. (2)]:∂hvi= (a − 1)hvi − bhv2 vi − γv hK(x)vi,
∂t

(C1)

where the brackets denote spatial averaging. Note that∇p, ∇2 v, and v · ∇v = ∇ · (vv) vanish because of the
π 2
periodic boundary condition. For k2 > 21 · ∆x
, every
π 2,thing is set to 0 to avoid aliasing. For 0 < k2 ≤ 12 · ∆x∗ω̂ is corrected to satisfy ω̂(k) = ω̂(−k) and v is recovered from the corrected ω̂. The correction is essential forstable calculation; without that, accumulated numerical∗errors violate ω̂naive (k) = ω̂naive (−k) , which is necessary
for real-valued ω. For v, GPU solves the equation for the
stream function ψ(x), ∇2 ψ = −ω, and obtains v from ψ.
This technique is valid for the two-dimensional incompressible flow studied here and ∇ · v = 0 is automaticallysatisfied. Without using this, one may need to solve Poisson’s equation for the pressure, which is severely timeconsuming, but here we can bypass it just by dividing ψ̂by k2 .
Appendix D: processing mask K(x)As described in the previous section, to avoid aliasing, it is necessary to apply a low-pass filter after realspace multiplications. This requirement also applies tothe mask K(x). Thus, we have to design an appropriate mask field K(x) such that it has no high-wavelength
modes, and at the same time it is nonnegative, K ' 0
inside the region of interest and K ' 1 otherwise. Thiscan be automated as follows:
Algorithm 2 mask processing scheme
Input: Knaive
1: K̂naive = DFT[Knaive ]2: for all k do
π 2
3:
if k2 > 14 · ∆xthen
4:K̂naive (k) = 0
5:
else
6:
modify K̂naive (k) (optional)
7:
end if
8: end for
9: K = iDFT[K̂naive ]2
Output: processed K(x)To begin with, Algorithm 2 computes DFT of Knaive ,followed by a low-pass filtering at a 41 cutoff. This cutoff
ensures that the output K(x) does not violate the 21 cutoff rule for antialiasing (we remark below why 14 ). Then,
optionally, one may include additional low-pass operations (or any other processing) in the else block, as long∗as K̂naive (k) = K̂naive (−k) is satisfied. In the last step,
K is substituted by squared iDFT of K̂naive . This ensuresK ≥ 0, without changing the essential shape of the mask.

elapsed time [ms/step]

13

104

unchanged. The color scale is identical to that used inFig. 1(b).

CPU
GPU

Movie S1: Vorticity field from a simulation with theemoji boundary. See Fig. 1 of the main text.
Movie S2: Vorticity fields for the circular confinementwith radii R = 5.2, 5.4, . . . , 8.0.

103
102
101
102.5

103.0

103.5

104.0

N

Movie S3: Velocity (arrows) and vorticity fields (color)
for R = 5.2 (single-vortex stationary state). This
video includes the transient from a random initialcondition.

FIG. 9: Time required for an iteration of the Euler
method versus lattice size N , compared between CPU(Intel Xeon W-2295) and GPU (NVIDIA A6000). Thisshows that GPU is approximately 60 times faster forlarge N .

Movie S4: Velocity (arrows) and vorticity fields (color)
for R = 5.4 (vortex-pair oscillatory state).

It is because of this square that we chose the 14 cutoff for
Knaive , to satisfy the 21 rule in the finally obtained K(x).
We used a binary mask for the input. For the circularcase, we set Knaive (x) = 1 for kxk ≥ R and 0 otherwise.
For the dumbbell case, Knaive was obtained by logicalAND of two circles. We did not use any optional operations because any additional cutoff makes the difference
between the designed geometry and the obtained K(x)larger.

Movie S6: Velocity (arrows) and vorticity fields (color)
for R = 5.8 (quasiperiodic state).

Appendix E: why GPU?

The core observation is that any operation in the calculation scheme explained above reduces to an element-wisearithmetic operation either in the real or Fourier space.
Element-wise operations can easily be parallelized, andusing cuFFT library the DFT becomes massively faster.
Note that any if statement should be purged from GPU
codes for better performance, and in the present scheme
it was possible to eliminate all the conditional branchesfrom the performance-critical components.
As a result, although the computational cost of ourlarge-scale simulations was far more massive comparedto the previous studies, GPU was able to deal with them
fast enough as shown in Fig. 9.

Movie S5: Velocity (arrows) and vorticity fields (color)
for R = 5.6 (chaotic state).

Movie S7: Velocity (arrows) and vorticity fields (color)
for R = 7.2 (periodic state).
Movie S8: Velocity (arrows) and vorticity fields (color)
for R = 7.6 (chaotic state).

Appendix G: code availability

The codes used in this work and associated information
are available upon request.

ACKNOWLEDGMENTS

In Movies S2-S8, scales of the arrows for the velocity field and the colors for the vorticity field are kept

We are grateful to Yusuke T. Maeda, Kazusa Beppu,
Lailai Zhu for helpful discussions, Susumu Goto, Yuji
Hattori, and Soshi Kawai for discussions from the perspective of fluid mechanics, Henning Reinken and Ryuna
Nagayama for debugging, and Taichi Taniguchi forrefactoring. The computation in this paper has been
done partly by Kugui GPU supercomputer (ISSP, Univ.
Tokyo) and nekoya/ai cluster (ipi, Univ. Tokyo). This
work is supported in part by JSPS KAKENHI Grant
Numbers JP19K23422, JP19H05800 and JP20K14426,
and JST PRESTO Grant Number JPMJPR21O8, Japan.

[1] L. Landau and E. Lifshitz, Fluid Mechanics, 2nd ed. (Elsevier, 1987).

[2] J. P. Eckmann, Roads to turbulence in dissipative dynamical systems, Rev. Mod. Phys. 53, 643 (1981).

Appendix F: movie captions

14[3] R. Alert, J. Casademunt, and J.-F. Joanny, Active turbulence, Annu. Rev. Condens. Matter Phys. 13, 143 (2022).
[4] H. Chaté, Dry aligning dilute active matter, Annu. Rev.
Condens. Matter Phys. 11, 189 (2019).
[5] A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot,and D. Bartolo, Emergence of macroscopic directed motion in populations of motile colloids, Nature 503, 95(2013).
[6] D. Nishiguchi, K. H. Nagai, H. Chaté, and M. Sano,
Long-range nematic order and anomalous fluctuations
in suspensions of swimming filamentous bacteria, Phys.
Rev. E 95, 020601 (2017).
[7] J. Iwasawa, D. Nishiguchi, and M. Sano, Algebraic correlations and anomalous fluctuations in ordered flocks ofjanus particles fueled by an ac electric field, Phys. Rev.
Res. 3, 043104 (2021).
[8] T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann,and Z. Dogic, Spontaneous motion in hierarchically assembled active matter, Nature 491, 431 (2012).
[9] D. Nishiguchi and M. Sano, Mesoscopic turbulence andlocal order in janus particles self-propelling under an acelectric field, Phys. Rev. E 92, 052309 (2015).
[10] A. Creppy, O. Praud, X. Druart, P. L. Kohnke, and
F. Plouraboué, Turbulence of swarming sperm, Phys.
Rev. E 92, 032722 (2015).
[11] C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos,L. Giomi, and P. Silberzan, Turbulent dynamics of epithelial cell cultures, Phys. Rev. Lett. 120, 208101 (2018).
[12] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher,
R. E. Goldstein, H. Löwen, and J. M. Yeomans, Mesoscale turbulence in living fluids, Proc. Natl. Acad. Sci.
USA 109, 14308 (2012).
[13] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink,
M. Bär, and R. E. Goldstein, Fluid Dynamics of Bacterial
Turbulence, Phys. Rev. Lett. 110, 228102 (2013).
[14] J. Dunkel, S. Heidenreich, M. Bär, and R. E. Goldstein, Minimal continuum theories of structure formationin dense active fluids, New J. Phys. 15, 045016 (2013),
arXiv:1208.4464.
[15] H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler,
and R. E. Goldstein, Confinement stabilizes a bacterial
suspension into a spiral vortex, Phys. Rev. Lett. 110,
268102 (2013).
[16] H. Wioland, E. Lushi, and R. E. Goldstein, Directedcollective motion of bacteria under channel confinement,
New J. Phys. 18, 075002 (2016).
[17] K. Doxzen, S. R. K. Vedula, M. C. Leong, H. Hirata, N. S.
Gov, A. J. Kabla, B. Ladoux, and C. T. Lim, Guidanceof collective cell migration by substrate geometry, Integr.
Biol. 5, 1026 (2013).
[18] K. T. Wu, J. B. Hishamunda, D. T. Chen, S. J. DeCamp, Y.-W. Chang, A. Fernández-Nieves, S. Fraden,
and Z. Dogic, Transition from turbulent to coherent flows
in confined three-dimensional active fluids, Science 355,
1284 (2017).
[19] A. Opathalage, M. M. Norton, M. P. N. Juniper,
B. Langeslay, S. A. Aghvami, S. Fraden, and Z. Dogic,
Self-organized dynamics and the transition to turbulenceof confined active nematics, Proc. Natl. Acad. Sci. USA
116, 4788 (2019).
[20] D. Nishiguchi, I. S. Aranson, A. Snezhko, and A. Sokolov,
Engineering bacterial vortex lattice via direct laser
lithography, Nat. Commun. 9, 4486 (2018).

[21] H. Reinken, D. Nishiguchi, S. Heidenreich, A. Sokolov,
M. Bär, S. H. L. Klapp, and I. S. Aranson, Organizing bacterial vortex lattices by periodic obstacle arrays,
Commun. Phys. 3, 76 (2020).
[22] N. Figueroa-Morales, M. M. Genkin, A. Sokolov, and I. S.
Aranson, Non-symmetric pinning of topological defects inliving liquid crystals, Commun. Phys. 5, 301 (2022).
[23] H. Wioland, F. G. Woodhouse, J. Dunkel, and R. E.
Goldstein, Ferromagnetic and antiferromagnetic order inbacterial vortex lattices, Nat. Phys. 12, 341 (2016).
[24] K. Beppu, Z. Izri, J. Gohya, K. Eto, M. Ichikawa, and
Y. T. Maeda, Geometry-driven collective ordering of bacterial vortices, Soft Matter 13, 5038 (2017).
[25] K. Beppu, Z. Izri, T. Sato, Y. Yamanishi, Y. Sumino, and Y. T. Maeda, Edge current and pairing
order transition in chiral bacterial vortices, Proc.
Natl. Acad. Sci. USA 118, e2107461118 (2021),
https://www.pnas.org/doi/pdf/10.1073/pnas.2107461118.
[26] L. Giomi, L. Mahadevan, B. Chakraborty, and M. F. Hagan, Excitable Patterns in Active Nematics, Phys. Rev.
Lett. 106, 218101 (2011).
[27] L. Giomi, M. C. Marchetti, and T. B. Liverpool, Complex
Spontaneous Flows and Concentration Banding in Active
Polar Films, Phys. Rev. Lett. 101, 198101 (2008).
[28] L. Giomi and M. Cristina Marchetti, Polar patterns in
active fluids, Soft Matter 8, 129 (2012).
[29] A. Doostmohammadi, T. N. Shendruk, K. Thijssen, andJ. M. Yeomans, Onset of meso-scale turbulence in active
nematics, Nat. Commun. 8, 15326 (2017).
[30] M. Sano and K. Tamai, A universal transition to turbulence in channel flow, Nat. Phys. 12, 249 (2016).
[31] G. Lemoult, L. Shi, K. Avila, S. V. Jalikop, M. Avila, andB. Hof, Directed percolation phase transition to sustained
turbulence in couette flow, Nat. Phys. 12, 254 (2016).
[32] See Supplemental Material at [URL] for Movies S1-Sx.
[33] TwitterOSS, twemoji, https://github.com/twitter/
twemoji.
[34] E. Lushi, H. Wioland, and R. E. Goldstein,
Fluid
flows
createdbyswimming
bacteria
drive self-organization in confined suspensions,
Proc. Natl. Acad. Sci. USA 111, 9733 (2014),
https://www.pnas.org/doi/pdf/10.1073/pnas.1405698111.
[35] S. H. Strogatz, Nonlinear Dynamics and Chaos: With
Applications to Physics, Biology, Chemistry, and Engineering, Studies in Nonlinearity (Westview Press, 2001).
[36] I. Shimada and T. Nagashima, A Numerical Approachto Ergodic Problem of Dissipative Dynamical Systems,
Prog. Theor. Phys. 61, 1605 (1979).
[37] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn,
Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method forcomputing all of them. Part 1: Theory, Meccanica 15, 9(1980).
[38] D. Ruelle and F. Takens, On the nature of turbulence,
Commun. Math. Phys. 20, 167 (1971).
[39] J. M. Hyman and B. Nicolaenko, The kuramotosivashinsky equation: A bridge between pde’s and dynamical systems, Physica D 18, 113 (1986).
[40] S. Kida, M. Yamada, and K. Ohkitani, A route to chaos
and turbulence, Physica D 37, 116 (1989).
[41] M. Besse, H. Chaté, and A. Solon, Metastability ofConstant-Density Flocks, Phys. Rev. Lett. 129, 268003(2022).

15[42] J. Toner and Y. Tu, Long-range order in a twodimensional dynamical XY model: How birds fly together, Phys. Rev. Lett. 75, 4326 (1995), arXiv:adaporg/9506001.
[43] J. Toner and Y. Tu, Flocks, herds, and schools: A quantitative theory of flocking, Phys. Rev. E 58, 4828 (1998).

[44] J. Toner, Reanalysis of the hydrodynamic theory of fluid,polar-ordered flocks, Phys. Rev. E 86, 031918 (2012),
arXiv:1204.4527v1.
[45] H. K. Versteeg and W. Malalasekera, An introduction tocomputational fluid dynamics: the finite volume method(Pearson education, 2007).