Vortex reversal is a precursor of confined bacterial turbulence
Daiki Nishiguchi,1, ∗ Sora Shiratani,1 Kazumasa A. Takeuchi,1, 2 and Igor S. Aranson3, 1
1

arXiv:2407.05269v1 [cond-mat.soft] 7 Jul 2024

Department of Physics, The University of Tokyo,
7–3–1 Hongo, Bunkyo-ku, Tokyo 113–0033, Japan
2
Institute for Physics of Intelligence, The University of Tokyo,
7–3–1 Hongo, Bunkyo-ku, Tokyo 113–0033, Japan3
Departments of Biomedical Engineering, Chemistry, and Mathematics,
The Pennsylvania State University, University Park, Pennsylvania 16802, USA(Dated: July 9, 2024)
Active turbulence, or chaotic self-organized collective motion, is often observed in concentrated
suspensions of motile bacteria and other systems of self-propelled interacting agents. To date, there
is no fundamental understanding of how geometrical confinement orchestrates active turbulence
and alters its physical properties. Here, by combining large-scale experiments, computer modeling,
and analytical theory, we have discovered a generic sequence of transitions occurring in bacterial
suspensions confined in cylindrical wells of varying radii. With increasing the well’s radius, weobserved that persistent vortex motion gives way to periodic vortex reversals, four-vortex pulsations,and then well-developed active turbulence. Using computational modeling and analytical theory, we
have shown that vortex reversal results from the nonlinear interaction of the first three azimuthalmodes that become unstable with the radius increase. The analytical results account for our key
experimental findings. To further validate our approach, we reconstructed equations of motion fromexperimental data. Our findings shed light on the universal properties of confined bacterial active
matter and can be applied to various biological and synthetic active systems.

Interacting self-propelled particles, often termed active matter, exhibit a remarkable tendency to selforganization and the onset of collective behavior. Being
intrinsically out-of-equilibrium, active matter systems
exhibit a slew of collective phenomena such as the spontaneous onset of long-range order [1–5], odd viscoelasticity [6], rectifications of chaotic flows [7–10], and reductionof the effective viscosity [11, 12]. One of the most visible manifestations of collective dynamics in active mattersystems is the emergence of self-sustained spatiotemporal chaotic flows termed active turbulence [13–18]. In
stark contrast to conventional Navier-Stokes turbulence,
active turbulence, occurring for essentially zero Reynolds
numbers, is characterized by the well-defined characteristic length scale. In the case of bacterial turbulence, thisscale corresponds to typical vortex size, which is about
40–50 µm [14, 15]. The existence of the typical vortex
size allows transforming bacterial motion into stable vortex arrays under geometrical confinements [19–24] or in
the presence of periodic obstacles [9, 10].
Experimental and computational studies of selforganization of bacterial and related active systems haveshown that strong confinement, e.g., a cylindrical well,
may suppress active turbulence and generate persistent
vortex motion [19, 20, 22–24]. However, a fundamental question on the nature of the transition from orderedstates under strong confinement to chaotic motion in unconstrained systems remains open. Answering this question will shed light on intricate fundamental mechanismsof self-organization in a broad class of active systems under confinement.
In the context of active nematics exemplified by
microtubules-motors assays, multiple experimental andnumerical studies interrogated a transition from ordered
quasi-stationary states to chaotic motion occurs under
the confinement in channels, rings, and wells [25–30].
The primary observation is that the instability of static
nematic configuration occurs via unbinding and subsequent chaotic motion of half-integer topological defects.
In polar active systems such as bacterial suspensions, experimental investigations have been hindered by the difficulty in resolving the detailed dynamics very close to
the transition point and the necessity of long-time measurements for evaluating the vortex stability.
Here, we examine the route to active turbulence by
combining large-scale experiments, high-resolution numerical modeling, and analytical theory. We focused ona well-characterized active system: suspensions of swimming bacteria [5]. We confined the suspensions into anarray of isolated cylindrical wells comparable to the sizeof individual vortices. We systematically varied the wells’
radii to characterize the transition from stabilized vortices to bacterial turbulence. Increasing the well radius,
we have detected reversals of vortex rotation as the first
instability from a stable vortex. The reversals were also
captured as periodic oscillations in our numerical simulations and analytical theory, unraveling a robust fundamental mechanism for the onset of polar active turbulence. It differs from the reversals caused by viscoelasticity of the suspending fluid [31]. Our analysis revealed
that the reversal originated from the nonlinear interaction of the three lowest azimuthal modes near the threshold of linear instability. To further validate our theoretical arguments, we reconstructed equations of motion
from experiential data. Our studies indicate that thevortex reversal is a generic precursor of turbulence-like2behavior in bacterial and related active systems. Ourfindings provide insights into how geometrical confinement orchestrates spatiotemporal organization in a broadclass of active systems.
Results

Experiment. We conducted experiments with suspensions of swimming bacteria confined in cylindrical wells,
Fig. 1(a,b). The height of the wells was set to 30 µm,
which is smaller than the typical length scale of collective motion, ensuring effectively two-dimensional dynamics within each well. Experiments were conducted simultaneously in an array of isolated wells of different radii(≈ 400 wells in total), see Fig. 1(d), Fig. S1, and Supplementary Movies 1–7. We observed stabilized vortices
with steady rotational directions within the wells withsmall radii. For larger radii, the vortices exhibited a transition to unsteady configurations with reversing rotation
directions. This observation is exemplified by the instantaneous vorticity field ω(r, t) = ẑ · [∇ × v(r, t)] shown in
Fig. 1(c), where ẑ is the unit vector in the z-direction. As
one sees from Fig. 1(c), the smaller wells hosted a single
stabilized vortex with persistent rotation, with the velocity and vorticity profiles shown in Fig. 1(e) and Fig. 2(a).
In contrast to previous studies on bacterial suspensions
confined in water-in-oil droplets [19, 20], we did not observe any counter-rotating edge flows, suggesting different boundary conditions for collective motion.
To quantify the vortex rotation direction, we defined aspin variable for each well as,Xẑ ·(r − ri ) × v(r, t)
Si (t) :=r∈i-th wellXr∈i-th well|r − ri |,(1)

where ri is the center of the i-th well, and the summations
run over the area of the i-th well. As shown in Fig. 1(e,f),
the spins for the small wells stayed almost constant andrarely flipped their signs over time, while the spins for thelarger wells persistently alternated their signs, reflecting
the reversals of vortices. The spin probability distribution for such a well with reversals exhibits a bimodal
distribution, indicating the presence of two states withclockwise (CW, Si < 0) and counterclockwise (CCW,Si > 0) rotations [Fig. 2(b)]. Contrary to Ref. [24], oursetup uses two symmetric surfaces for the top and bottom to compensate for systematic bias in the rotation
direction. Thus, the fraction of CW rotations as a function of the well’s radius was always ≈ 0.5, Fig. 2(d). The
absence of bias is crucial for characterizing the vortex reversals. The transition from a single stabilized vortex toreversing vortices was inspected through the spin correlation time, defined as the time at which the autocorrelation function of the spins decayed to 1/e. The correlation

time has successfully captured the transition at the radiiof approximately 46 to 48 µm [Fig. 2(c)]. For large wells’
radii, four-vortex pulsating states were observed as well,Fig. 1(c). The pulsation was characterized by the kinetic
energy of azimuthal modes corresponding to 2n vortices
within a well (see Supplementary Note 4),
Z
Z R
2
1
exp
(2)
dθe−inθ v(r, θ) .
mn =
drr
2π
0
By this mode analysis, we probed anti-phase oscillationof the modes n = 1 and n = 2, see Fig. 1(g).
Computational modeling. We performed numericalsimulations using a phenomenological active fluid model,
the Toner-Tu-Swift-Hohenberg equation (TTSHE) [5, 10,
15–17]. The TTSHE qualitatively captures the bulk
properties of polar active turbulence. It can describe thetransformation of bacterial turbulence into stable vortex
arrays in the presence of periodic obstacles [9, 10]. In thevorticity representation, the dimensionless TTSHE is of
the form[10]:∂ω+ λv · ∇ω = aω − b∇ × |v|2 v∂t2− 1 + ∇2 ω − γv ∇ × [K(r)v] − γω K(r)ω,(3)

where λ, a, and b are constants, K(r) ≥ 0 is a scalar fieldthat dampens v and ω outside the well (K ≃ 1) without
affecting the inside (K ≃ 0), and γv,ω > 0 are damping
coefficients. In this dimensionless form, the vortex characteristic size is 2π. Following Ref. [10], we adopt theparameter values (λ, a, b, γv , γω ) = (9, 0.5, 1.6, 40, 4) and
impose three boundary conditions on well’s wall,v = 0, ω = 0 at r = R.

(4)

Compared with the Navier-Stokes equation, the extra
boundary condition ω = 0 is imposed due to the higherorder differential operator (∇4 ) in Eq. (3). We solved
Eq. (3) with the above boundary conditions in two dimensions by the pseudospectral method, see Methods.
Our simulations successfully reproduced the entire sequence of transitions observed in experiments, Fig. 3(a)and Supplementary Movies 8–13. We have found a single stable vortex for small radii. As the radius increases,
the vortex becomes destabilized and yields a periodically reversing two-vortex state, see Fig. 3(b,c). It wasdemonstrated by the time series of the spin variable, seeFig. 3(b). Increasing the radius, the reversing two-vortex
state transforms into a pulsating four-vortex state, similarly to the experiment, Fig. 3(d).
Weakly-nonlinear analysis. We examined the linear
stability of Eq. (3) around v = 0, yielding∂t ω = aω−(1+P∞∇2 )2 ω. Its solution is of the form ω = −∞ exp(λn t)ωn ,ωn = (Gn+ Jn (kn+ r) + Gn− Jn (kn− r)) exp(inθ),(5)

where λn are the growth rates of the corresponding azimuthal modes, Gn± are constants, Jn are the Bessel

3

(a)

(b)

z y

2.5

~300 µm

2
1.5

z

x

(c)

(d)
30 µm

1

~500 µm

0.5

xR = 44.6 µm

R = 46.7 µm

0

R = 48.8 µm

[s-1]

-0.5
-1
-1.5

100 µm

-2
-2.5

(e)

20 µm
10
5
0
-5
-10
0

20 µm

R = 44.0 µm

radius R (0.5 µm step)(f)

R = 46.7 µm

R = 44.6 µm
50

100

150

10
5
0
-5
-10
0

20 µm

50

R = 50.5 µm

100

150

20 µm

(g)
R = 48.8 µm

20 µm

FIG. 1. Transitions from a stabilized vortex to reversing vortices and a four-vortex state. 3D schematics (a) andside view (b) of the experimental setup. (c) Typical vorticity profiles of a single stabilized vortex (R = 44.6 µm), reversing
vortices (R = 46.7 µm), and a four-vortex state (R = 48.8 µm). Vorticity field ω is overlaid on the experimental snapshots. The
color scales of the vorticity fields in all panels are identical and are indicated by the color bar in (d). (Supplementary Movies
2–4). (d) Experimental snapshot overlaid with the instantaneous vorticity field. Wells with the same radius are arranged
vertically, with the radius increasing from left to right. All the 119 wells within this image out of ∼ 400 wells within the whole
field of view (Fig. S1, Supplementary Movie 1) were used for analysis, see Supplementary Note 1 and Fig. S3 for the selection
criteria. (e,f) Time series of spins for the wells with the radii of 44.6 µm (e) and 46.7 µm (f), respectively. The instantaneous
velocity and vorticity fields are shown below the time series, with the colors of the rectangle corresponding to the time points
highlighted by colored circles in the time series (Supplementary Movie 2, 3, 5, 6). (g) Anti-phase relation of mode amplitudes
mexp
and mexpof the four-vortex state at R = 48.8 µm. The instantaneous vorticity fields are shown on the right of the time
1
2
series, with the colors of the rectangle corresponding to the time points highlighted by colored circles in the time series, see
Supplementary Note 4 and Supplementary Movies 3, 7.

p√functions, kn± = 1 ± a − λn . Applying the boundary conditions to ωn and solving the characteristic equation (see Supplementary Note 4), one finds the growth
rates λn vs radius R. The results are shown in Fig. 4(a).
For small enough R, all λn are negative, so that no vortex is excited. For R ⪆ 4.2, λ0 becomes positive, corresponding to the onset of the steady-state vortical motion
observed in computational modeling and experiment, seeFig. S7(b) for quantitative agreement between the nu-

merical and analytical solutions. Then, with the gradualincrease in R, higher rotational modes become unstable.
We find that vortex reversal occurs when the first twomodes with n = 0, ±1 are unstable, and the mode withn = ±2 is still stable but close to the threshold, compare
Figs. 4(a) and (e).
The radial vorticity and velocity profiles predicted by
the linear analysis, Eq. (5), are in excellent agreementwith the numerical solutions of Eq. (3), Figs. 4(b,c). Fur-

4

FIG. 2. Characterization of vortex states. (a) Timeaveraged velocity v (blue) and vorticity ω (red) profiles of astabilized vortex shown in Fig. 1(d). The vertical dashed line
corresponds to the radius detected from the image analysis,see Supplementary Note 1. The dashed lines are fits to theanalytical solutions, Eq. (5), for tangential velocity vθ (blue)and vorticity ω (red) fields. (b) Spin probability density function for the reversing vortex state shown in Fig. 1(f). (c) A
scatter plot of the spin correlation time and the well radii detected from the image analysis. The red line represents themoving median of the scatter plot. The horizontal dashedline corresponds to the experimental duration, 150 seconds,
see Supplementary Note 1 for details. (d) Fraction of CW
rotations as a function of the well radius. Error bars are thestandard errors.

thermore, fitting the theoretical expression, Eq. (5), forn = 0 to experimental vorticity and velocity profiles ofa stable vortex provides an excellent approximation as
well, see Fig. 2(a) and Supplementary Note 1.
Next, we approximate the solution to Eq. (3) as a sumof the three lowest azimuthal modes with n = 0, ±1, ±2,
ω = C(t)ω0 (r) + [A1 (t)eiθ ω1 (r) + A2 (t)e2iθ ω2 (r) + c.c.](6)Here ω0 , ω1 , ω2 are the eigenfunctions obtained from linear stability analysis. For definiteness, the eigenmodes
are normalized by their kinetic energy, see Supplementary Note 4, Eqs. (S35),(S36). C(t), A1 (t), A2 (t) areslowly-varying amplitudes that are derived from the corresponding orthogonality leading to a set of normal form
equations (7)-(9), see Methods.
Equations (7)-(9) faithfully reproduce the numericalresults from Eq. (3) without further approximation, seeFig. 4. Specifically, for small radii, Eqs. (7)-(9) reproduce
a stable vortex solution as shown in Fig. 3(a). Then,with the increase in R, the bifurcation to a limit cycle is faithfully captured. Furthermore, even the detailsof the time dependence of each azimuthal mode closely
agree with those of the numerical solutions of Eq. (3),see Figs. 4(d), S7(c,d). In the reversing vortex state, dis-played in Fig. 4(e), all three amplitudes C, A1 , A2 are
non-zero, see Supplementary Movie 14. With the further
increase in R, a transition from a reversing state to a pulsating four-vortex solution occurs, see Fig. 4(f) and Supplementary Movie 15. Here, the zero mode, n = 0, is suppressed, and the first and second modes A1 , A2 pulsate
in anti-phase. As shown in Fig. 1(g), this anti-phase relation was indeed observed experimentally, further demonstrating the quantitative agreements among the experimental, numerical, and analytical results. The normal
form analysis indicates that the transition to vortex reversals and other time-dependent states is a result of resonant nonlinear interaction among the three lowest azimuthal modes. This behavior only exists for sufficiently
large values of the nonlinear advection term λv · ∇ω in
Eq. (3), which controls the resonant three-mode interaction. No limit cycles were found for λ ⪅ 3.75.
Validating equations of motion. The use of the
TTSHE was validated through regression analysis of ourexperimental data; see Supplementary Note 2 for details.
Similar approaches were used in Refs. [32, 33]. In addition to the TTSHE, we tested another model for bacterial
turbulence, the Nikolaevskiy equation, which includes∇6 v term but no cubic nonlinearity |v|2 v nor linear termv [34–36]. The TTSHE outperformed the Nikolaevskiy
equation in terms of the residuals, justifying our numerical and theoretical approaches, see Figs. S4, S5, S6 andTables. S1, S2, S3. The regression for the two-vortex reversing state shown in Fig. 1(f) yields λdim = 1.69 ± 0.38
for the dimensional TTSHE, proving the presence of theadvection term with λ > 1, larger than λ = 1 for the
Naiver-Stokes equation. Transforming the TTSHE into
the form of Eq. (3) with characteristic values in the unconstrained bacterial turbulence (velocity V ≈50 µm/s,
length scale L ≈40 µm, and time scale T ≈0.5 s) yields
λnondim = VLT λdim ≈ 4.2 in the dimensionless TTSHE. It
is consistent with our theoretical prediction of λ ⪆ 3.75
for the onset of oscillations.
Concluding remarks

We observed a generic route to active turbulence inconfined suspensions of swimming bacteria: a single
steady vortex gives way to a reversing vortex pair, four
pulsating vortices, and then to a well-developed spatiotemporal chaos. The fact that the entire bifurcation
sequence is reproduced by a generic phenomenologicalmodel for active turbulence reveals the universal fundamental mechanism governing the transition: resonantinteraction of the three lowest azimuthal modes associated with cylindrical confinement. Furthermore, the onset of the periodic reversal relies on the finite value of
the Navier-Stokes-like advection term in the phenomenological model of active turbulence [10, 15, 17, 37]. The
regression of experimental data also reliably corroborates

5

(a)

(b)

(c)

(d)

FIG. 3. Computational modeling using the TTSHE. (a) Vorticity profiles obtained in the numerical simulations. Typicalsnapshots of a single stabilized vortex (R = 5.2, Supplementary Movie 8), periodically reversing two-vortex state (R = 5.6,
Supplementary Movie 9), and a pulsating four-vortex state (R = 6.4, Supplementary Movie 10) are shown. The color bar in
all panels is the same. (b) Time series of the spin for the reversing two-vortex state (R = 5.35, Supplementary Movie 11). The
instantaneous vorticity fields are shown as insets, with the colors of the rectangle corresponding to the time points highlighted
by colored circles in the time series. For computational convenience, instead of using the spin defined in Eq. (1), we plot theamplitude C of the zeroth azimuthal mode defined in Eq. (6), because C is proportional to the spin (see Supplementary Note 4).
(c) Azimuthal mode decomposition of the instantaneous vorticity field
R for the reversing two-vortex Rstate shown in the middle
panel of (a) (R = 5.6, Supplementary Movie 12), which is defined as ω(r, θ) dθ /2π for n = 0 and e−inθ ω(r, θ) dθ /2π + c.c.
otherwise. (d) Snapshots of the pulsating four-vortex state (R = 6.2, Supplementary Movie 13).

the presence of the advection term with its coefficientλ > 1 in the effective equation. These findings suggestthat the observed transitions should also occur in a broadclass of active self-propelled systems under confinement.
This robust mechanism is presumably responsible for the
onset of reversing edge currents numerically observed inRef. [38] and is not sensitive to the details of boundary
conditions or geometry [37]. Furthermore, the observed
traditions occur in a Newtonian fluid environment. Viscoelasticity or anisotropy may only affect the details of
the transition [31, 39, 40]. This generic mechanism isbased on the three-mode resonant interaction and should
be relevant for the variety of biological and synthetic active systems, e.g., Janus colloids [3, 4, 41].
Another intriguing aspect is the effect of chirality.
Since bacteria are chiral objects due to counter-rotationof the body and helical flagella [42–45], there could be
an asymmetry between CW/CCW rotating vortices [24].
In this work, a sustained effort was undertaken to makethe upper/bottom surfaces of the wells as identical as
possible to suppress the asymmetry. While a minor chiral shift does not affect the transition sequence, it could
introduce slightly different thresholds for the onsets ofvortex oscillations of opposite chirality.
The current experiment is unavoidably susceptible tostrong fluctuations discarded in the theoretical description. For example, the number of bacteria within a sin-gle microscopic well is about ∼ 104 bacteria/well. The
dynamics of such a small bacterial population is intrinsically stochastic. Therefore, understanding how noise
affects the nature of transitions and exploring ways totame and control the fluctuating active dynamics would
be of interest to future studies.
Finally, the controls and rectification of vortices in confined active matter open up new possibilities for engineering active matter. For instance, weak coupling between neighboring wells may realize a “bacterial latticeclock”, in which reversing vortex pairs synchronize andexhibit higher regularity and persistence. The reversing
or pulsating vortices may be useful for mixing at low
Reynolds numbers. Taming the fluctuations in active
systems based on the fundamental instability uncoveredin this work provides novel design principles for functioning active devices, such as biosensors or microrobotic
swarms for targeted drug delivery, precision surgery, or
detoxification [46, 47].

Methods

Experimental details. Bacteria Bacillus subtilis(strain: 1085) were grown in Terrific Broth (T9179,
Sigma-Aldrich) growth medium until optical density(OD) achieved OD600 nm ≈1. After concentrating the

6

(a)

(b)

(c)
n

0.5

1
2

0.2

0.4

0
1
2

0.2

−0.4

0.3

vθn

−0.2

vrn

λn

0.0

n

0.4

0.2

−0.2

n

−0.6
−0.8
−1.0

4.5

5.0

(d)

5.5

6.0

0
1
2
6.5

0.1
−0.4

0.0

1

2

3

4

5

6

7

R

(e)

R

0.0

1

2

3

(f)

4

5

6

7

R

1.00

|A2|

theory
num.

0.4

0.75

0.3

0.2

0.2

0.0

0.1
0.0
−0.4

0.50

−0.2
−0.2

0.0

C

0.2

0.4

0.45
0.40
0.35
0.30
0.25
|A1|

−0.4
0

C
|A1|
|A2|

500

1000

1500

0.25
0.00
0

C
|A1|
|A2|

50

100

150

200

t

t

FIG. 4. Analytical results. (a) The growth rates λn vs radius R for n = 0, 1, 2. (b,c) Comparison of velocity profiles forthe azimuthal modes with n = 0, 1, 2 obtained from the linear theory (solid), Eq. (5), and the TTSHE simulations (dashed).
(d) Comparison of trajectories in 3D phase space obtained by the solution of Eqs. (7)-(9) and the TTSHE. (e,f) Amplitudes
C, A1 , A2 vs time obtained from Eqs. (7)-(9) for R = 5.9 (e) and R = 7.0 (f), see Supplementary Movies 14, 15.

suspension 180-fold, it was sandwiched between two thinpolydimethylsiloxane (PDMS) membranes to facilitate
sufficient oxygen supply for sustaining high bacterial
motility. The bottom PDMS membrane was patternedwith 30-µm-deep multiple microscopic wells with theradius ranging from 44 to 51 µm with 0.5 µm increments [Figs. 1(a,b,d), S1, S2]. To overcome systematic
errors arising from different preparations of bacterial
cultures and slight density variations caused during
the confinement process, we simultaneously observed∼ 400 wells (19 radii, ∼ 20 wells for each radius) ina single field of view by using an inverted microscope
equipped with a large-sensor sCMOS camera (Kinetix,
Teledyne Photometrics, 3200×3200 pixels) and a 10×
objective lens, realizing the 2.1 mm × 2.1 mm fieldof view (Fig. S1). It allowed resolving the bacterial
dynamics very close to the transition point. We captured
the movies at 50 fps for 150 seconds and analyzed thebacterial velocity fields v(r, t) = (vx , vy ) using theparticle image velocimetry (PIV), see Supplementary
Note 1 for the detailed protocols.
Computational details. Equation 3 was solved by thepseudospectral method in a two-dimensional periodic
40.96 × 40.96 domain, discretized as the 8192 × 8192
square lattice. Spatial derivatives were handled by thefast Fourier transform, see Supplementary Note 3. Time
update was performed in the Fourier space, with the time
step ∆t = 0.01. To accelerate simulations, we performed
the whole computation on GPUs (NVIDIA RTX A6000

or A100).
The damping wall implemented in Eq. (3) with the kernel K(r) permits some leakage outside of the well radius
R. We calibrated R to account for the leakage and defined the effective radius Reff where the velocity
R and vorticity vanish. Reff is calculated as the root of vθ (r, θ) dθ(the zeroth azimuthal mode), Supplementary Fig. S7. We
obtained Reff − R ≈ 0.5. The zeroth mode amplitude
C in Eq. (6) provides a convenient measure of the spin
(Eq. (1)) up to a certain prefactor. For the details of
the numerical mode decomposition and related quantities, see Supplementary Note 4.
Normal form equations. Substituting Eq. (6) into
Eq. (3), and implementing the orthogonality conditions,
we obtain the set of equations for amplitudes C, A1 , A2
∂t C = λ0 C − c1 C 3 − c2 C|A1 |2 − c3 C|A2 |2− 2c4 ReA2 A2∗
1

2

2∂t A1 = λ1 A1 − b1 A1 |A1 | − b2 A1 C − b3 A1 |A2 |(7)
2− b4 CA2 A∗1 + δ1 A1 C + γ1 A2 A∗1(8)− a4 CA21 + δ2 A2 C + γ2 A21(9)∂t A2 = λ2 A2 − a1 A2 |A2 |2 − a2 A2 C 2 − a3 A2 |A1 |2

λ0,1,2 are the linear growth rates; other coefficients areintegrals over the nonlinearities, Supplementary Note 4.
Acknowledgments. This work is supported in partby JSPS KAKENHI Grant Numbers JP19H05800,
JP20K14426, JP23K25838, JP24K00593, JP24KJ0890,
JST PRESTO Grant Number JPMJPR21O8, Japan, and7
the JSPS Core-to-Core Program “Advanced core-to-core
network for the physics of self-organizing active matter(JPJSCCA20230002)”. The research of I.S.A. was supported by the NSF award PHY-2140010. I.S.A. also acknowledges the visiting professorship for the Global Science Graduate Course -B, at The University of Tokyo,
where main part of this work was done. S.S. acknowledges support from FoPM, WINGS Program, The University of Tokyo. The computation in this paper wasdone partly by the nekoya/ai cluster at the Institute forPhysics of Intelligence, The University of Tokyo, and thefacilities of the Supercomputer Center, the Institute for
Solid State Physics, The University of Tokyo.
Competing interests. The authors declare no competing interests.
Data availability. All the experimental and numerical data, as well as the relevant codes and scripts, areavailable upon request.
Author contributions. D.N. conceived and initiatedthe project, performed the experiments, and analyzed theexperimental data. S.S. performed the numerical simulations and analyzed the numerical data. I.S.A. developed
the analytical theory. All authors discussed and interpreted the results, and wrote the manuscript.

∗
nishiguchi@noneq.phys.s.u-tokyo.ac.jp[1] D. Nishiguchi, K. H. Nagai, H. Chaté, and M. Sano,
Long-range nematic order and anomalous fluctuations insuspensions of swimming filamentous bacteria, PhysicalReview E 95, 020601(R) (2017).
[2] H. Chaté, Dry Aligning Dilute Active Matter, Annual
Review of Condensed Matter Physics 11, 189 (2020).
[3] J. Iwasawa, D. Nishiguchi, and M. Sano, Algebraic correlations and anomalous fluctuations in ordered flocks ofJanus particles fueled by an AC electric field, PhysicalReview Research 3, 043104 (2021).
[4] D. Nishiguchi, Deciphering Long-Range Order in Active
Matter: Insights from Swimming Bacteria in Quasi-2Dand Electrokinetic Janus Particles, Journal of the Physical Society of Japan 92, 121007 (2023).
[5] I. S. Aranson, Bacterial active matter, Reports on
Progress in Physics 85, 076601 (2022).
[6] M. Fruchart, C. Scheibner, and V. Vitelli, Odd Viscosity
and Odd Elasticity, Annual Review of Condensed Matter
Physics 14, 471 (2023).
[7] A. Sokolov, M. M. Apodaca, B. A. Grzybowski, and I. S.
Aranson, Swimming bacteria power microscopic gears.,
Proceedings of the National Academy of Sciences of theUnited States of America 107, 969 (2010).
[8] C. O. Reichhardt and C. Reichhardt, Ratchet Effects in
Active Matter Systems, Annual Review of Condensed
Matter Physics 8, 51 (2017).
[9] D. Nishiguchi, I. S. Aranson, A. Snezhko, and A. Sokolov,
Engineering bacterial vortex lattice via direct laser
lithography, Nature Communications 9, 4486 (2018).
[10] H. Reinken, D. Nishiguchi, S. Heidenreich, A. Sokolov,
M. Bär, S. H. Klapp, and I. S. Aranson, Organizing bac-terial vortex lattices by periodic obstacle arrays, Communications Physics 3, 76 (2020).
[11] A. Sokolov and I. S. Aranson, Reduction of viscosity insuspension of swimming bacteria, Physical Review Letters 103, 148101 (2009).
[12] V. A. Martinez, E. Clément, J. Arlt, C. Douarche,
A. Dawson, J. Schwarz-Linek, A. K. Creppy, V. Škultéty,
A. N. Morozov, H. Auradou, et al., A combined rheometry and imaging study of viscosity reduction in bacterial
suspensions, Proceedings of the National Academy of Sciences 117, 2326 (2020).
[13] A. Sokolov, I. S. Aranson, J. O. Kessler, and R. E. Goldstein, Concentration dependence of the collective dynamics of swimming bacteria, Physical Review Letters 98,
158102 (2007).
[14] A. Sokolov and I. S. Aranson, Physical properties of collective motion in suspensions of bacteria, Physical Review Letters 109, 248109 (2012).
[15] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher,
R. E. Goldstein, H. Löwen, and J. M. Yeomans, Mesoscale turbulence in living fluids, Proceedings of the National Academy of Sciences of the United States of America 109, 14308 (2012).
[16] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink,
M. Bär, and R. E. Goldstein, Fluid Dynamics of Bacterial
Turbulence, Physical Review Letters 110, 228102 (2013).
[17] J. Dunkel, S. Heidenreich, M. Bär, and R. E. Goldstein,
Minimal continuum theories of structure formation indense active fluids, New Journal of Physics 15, 045016(2013).
[18] K. Qi, E. Westphal, G. Gompper, and R. G. Winkler,
Emergence of active turbulence in microswimmer suspensions due to active hydrodynamic stress and volume
exclusion, Communications Physics 5, 49 (2022).
[19] H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler,
and R. E. Goldstein, Confinement stabilizes a bacterial
suspension into a spiral vortex, Physical Review Letters
110, 268102 (2013).
[20] E. Lushi, H. Wioland, and R. E. Goldstein, Fluid flows
created by swimming bacteria drive self-organizationin confined suspensions, Proceedings of the National
Academy of Sciences 111, 9733 (2014).
[21] H. Wioland, E. Lushi, and R. E. Goldstein, Directedcollective motion of bacteria under channel confinement,
New Journal of Physics 18, 075002 (2016).
[22] H. Wioland, F. G. Woodhouse, J. Dunkel, and R. E.
Goldstein, Ferromagnetic and antiferromagnetic order inbacterial vortex lattices, Nature Physics 12, 341 (2016).
[23] 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).
[24] K. Beppu, Z. Izri, T. Sato, Y. Yamanishi, Y. Sumino, and
Y. T. Maeda, Edge current and pairing order transitionin chiral bacterial vortices, Proceedings of the National
Academy of Sciences 118, e2107461118 (2021).
[25] K.-T. Wu, J. B. Hishamunda, D. T. Chen, S. J. DeCamp, Y.-W. Chang, A. Fernández-Nieves, S. Fraden,
and Z. Dogic, Transition from turbulent to coherent flows
in confined three-dimensional active fluids, Science 355,
1284 (2017).
[26] A. Opathalage, M. M. Norton, M. P. Juniper, B. Langeslay, S. A. Aghvami, S. Fraden, and Z. Dogic, Selforganized dynamics and the transition to turbulence ofconfined active nematics, Proceedings of the National

8
Academy of Sciences 116, 4788 (2019).
[27] J. Hardoüin, R. Hughes, A. Doostmohammadi, J. Laurent, T. Lopez-Leon, J. M. Yeomans, J. Ignés-Mullol, and
F. Sagués, Reconfigurable flows and defect landscape ofconfined active nematics, Communications Physics 2, 121(2019).
[28] J. Hardoüin, C. Doré, J. Laurent, T. Lopez-Leon,
J. Ignés-Mullol, and F. Sagués, Active boundary layers
in confined active nematics, Nature Communications 13,
6675 (2022).
[29] C. D. Schimming, C. Reichhardt, and C. Reichhardt,
Friction-mediated phase transition in confined active nematics, Physical Review E 108, L012602 (2023).
[30] M. Ravnik and J. M. Yeomans, Confined active nematic
flow in cylindrical capillaries, Physical Review Letters
110, 026001 (2013).
[31] S. Liu, S. Shankar, M. C. Marchetti, and Y. Wu, Viscoelastic control of spatiotemporal order in bacterial active matter, Nature 590, 80 (2021).
[32] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the
National Academy of Sciences 113, 3932 (2016).
[33] R. Supekar, B. Song, A. Hastewell, G. P. Choi, A. Mietke, and J. Dunkel, Learning hydrodynamic equations
for active matter from particle simulations and experiments, Proceedings of the National Academy of Sciences
120, e2206994120 (2023).
[34] I. A. Beresnev and V. N. Nikolaevskiy, A model for nonlinear seismic waves in a medium with instability, Physica
D: Nonlinear Phenomena 66, 1 (1993).
[35] M. I. Tribelsky and K. Tsuboi, New scenario for transition to turbulence?, Physical Review Letters 76, 1631
(1996).
[36] J. Słomka and J. Dunkel, Spontaneous mirror-symmetry
breaking induces inverse energy cascade in 3D active fluids, Proceedings of the National Academy of Sciences
114, 2119 (2017).

[37] S. Shiratani, K. A. Takeuchi, and D. Nishiguchi, Routeto turbulence via oscillatory states in polar active
fluid under confinement, arXiv preprint arXiv:2304.03306(2023).
[38] H. Matsukiyo and J.-i. Fukuda, Oscillating edge current
in polar active fluid, Physical Review E 109, 054604(2024).
[39] W. Liao and I. S. Aranson, Viscoelasticity enhancescollective motion of bacteria, PNAS Nexus 2, pgad291(2023).
[40] S. Zhou, A. Sokolov, O. D. Lavrentovich, and I. S. Aranson, Living liquid crystals, Proceedings of the National
Academy of Sciences of the United States of America
111, 1265 (2014).
[41] D. Nishiguchi and M. Sano, Mesoscopic turbulence andlocal order in Janus particles self-propelling under an acelectric field, Physical Review E 92, 052309 (2015).
[42] W. R. DiLuzio, L. Turner, M. Mayer, P. Garstecki,
D. B. Weibel, H. C. Berg, and G. M. Whitesides, Escherichia coli swim on the right-hand side, Nature 435,
1271 (2005).
[43] Marcos, H. C. Fu, T. R. Powers, and R. Stocker, Bacterial rheotaxis, Proceedings of the National Academy of
Sciences 109, 4780 (2012).
[44] G. Jing, A. Zöttl, É. Clément, and A. Lindner, Chiralityinduced bacterial rheotaxis in bulk shear flows, Scienceadvances 6, eabb2012 (2020).
[45] D. Grober, I. Palaia, M. C. Uçar, E. Hannezo, A. Šarić,
and J. Palacci, Unconventional colloidal aggregation inchiral bacterial baths, Nature Physics 19, 1680 (2023).
[46] J. Li, B. Esteban-Fernández de Ávila, W. Gao, L. Zhang,and J. Wang, Micro/nanorobots for biomedicine: Delivery, surgery, sensing, and detoxification, Science robotics
2, eaam6431 (2017).
[47] J. Li and M. Pumera, 3d printing of functional microrobots, Chemical Society Reviews 50, 2794 (2021).

Supplementary Material for“Vortex reversal is a precursor of confined bacterial turbulence”
Daiki Nishiguchi,1, ∗ Sora Shiratani,1 Kazumasa A. Takeuchi,1, 2 and Igor S. Aranson3, 1
1

arXiv:2407.05269v1 [cond-mat.soft] 7 Jul 2024

Department of Physics, The University of Tokyo,
7–3–1 Hongo, Bunkyo-ku, Tokyo 113–0033, Japan
2
Institute for Physics of Intelligence, The University of Tokyo,
7–3–1 Hongo, Bunkyo-ku, Tokyo 113–0033, Japan3
Departments of Biomedical Engineering, Chemistry, and Mathematics,
The Pennsylvania State University, University Park, Pennsylvania 16802, USA(Dated: July 9, 2024)I.

SUPPLEMENTARY NOTE 1: EXPERIMENTAL DETAILSA.

Bacterial culture

Bacteria Bacillus subtilis (strain 1085) were first grown on an LB (Luria-Bertani) agar plate (BD Difco LB BrothMiller + 1.5 w% agar) at 37◦ C. The grown bacterial colony was transferred into 13 mL of Terrific Broth (T9179,
Sigma-Aldrich) medium in a glass test tube, sealed with a cap, and cultured at 37◦ C at 200 rpm. When it reachedOD600 nm ≈0.8–1.0, the test tube was taken out from the incubator and placed at room temperature (23◦ C) for 45
minutes so that the bacteria became adjusted to the temperature before the experiments. The bacterial suspension
was centrifuged at 3000 rpm (1057 rcf) for 2 minutes and concentrated 180-fold by removing the supernatant. The
concentrated bacterial suspension was confined into the PDMS device for observation.
B.

Microfabrication

The fabrication of the microfluidic device was performed with the standard soft lithography. The photoresist SU-8
3050 was spread on a silicon wafer at the thickness of 30 µm, and circular patterns with multiple diameters weredrawn with a maskless aligner µMLA (Heidelberg Instruments) and developed following the standard protocol. Then,
dimethylpolysiloxane (PDMS, Sylgard 184, Dow) was poured onto the patterned silicon wafer to form the layer of the
thickness ∼ 500 µm. We also poured PDMS onto another non-patterned bare silicon wafer to make a ∼ 300-µm-thickflat PDMS membrane to be used as the top plate. The cured PDMS membranes were peeled off from the siliconwafers. The plasma cleaning was applied to the membranes to make the surface hydrophilic so that the bacterial
suspensions can smoothly enter the small circular wells.
We placed the concentrated bacterial suspension on the patterned PDMS membrane and then sandwiched thesuspension by placing the non-patterned PDMS membrane from the top. The use of the thin PDMS membranes forboth top and bottom substrates realized sufficient oxygen supply to the dense bacterial suspension, enabling longer
observation of highly active bacterial dynamics and the evaluation of vortex stability.
C.

Microscopy

The bright field images of the bacterial dynamics were captured with an inverted microscope (IX83, Evident)
equipped with 10x objective lens (UPLSXAPO, NA=0.4) and a large-field-of-view sCMOS camera (Kinetix, Teledyne
Photometrics, 3200×3200 pixels). With this setup, we captureed a field of view of about 2.1 mm × 2.1 mm with thespatial resolution of 0.65 µm/pixel at 50 fps for 150 seconds, enabling the simultaneous observation of multiple wellswith different radii with 0.5 µm increment (Fig. S1). This allowed us to avoid conducting multiple experiments fordifferent radii using different bacterial cultures, thereby preventing uncontrolled variability in the sample states.
In addition, simultaneous observation of ∼ 20 wells with the same designed radii Rdesign was crucial for our
purpose. The bacterial density within each well cannot be uniform due to the confinement process, which can be∗ nishiguchi@noneq.phys.s.u-tokyo.ac.jp

2

FIG. S1. Snapshot of the whole field of view taken by the camera. The wells aligned in each column have the same designed
radii and their radii are increasing from left to right. The wells whose centers were located within the green rectangular region
were used for analysis. This rectangular region is displayed in Fig.1(c) in the main text.

visually confirmed by the brightness of the bright field images (Fig. S1). The slight variations of the density largely
affect the dynamics of the bacterial collective motion very close to the transition point. This is the reason why, evenamong the wells with the same designed radii Rdesign , we observed variations of the spin correlation times [Fig. 2(c)]and the CW-CCW bias in the rotational directions (Fig. S3).

3

50

48

46

44

42
42

44

46

48

50

FIG. S2. Comparison between the designed radii Rdesign and the detected radii Rdetect of the wells. The wells in the field ofview (Fig. S1) contained wells with Rdesign ranging from 41.0 to 50.0 µm with 0.5 µm increment. The detected radii Rdetectare slightly larger than the designed radii Rdesign , with the average difference of 0.98 ± 0.05 µm (±: standard errors). Bluecircles are detection results. A black dashed line is the line Rdetect = Rdesign to guide the eye. Error bars: standard errors.
D.

Image processing

The obtained movie was processed by using the open-source MATLAB plugin PIVlab [1] to obtain the velocity field
vraw (x, y, t). The PIV interrogation box size was 16×16 pixels (10.4 µm ×10.4 µm, sufficiently smaller than the lengthscale of bacterial collective motion) with the step size of 8 pixels (5.2 µm, 50% overlap). Since the use of the PDMS
membranes for both substrates instead of glass substrates resulted in the deflection of the device
reduced opticalih  and221x2
resolutions, we applied a three-dimensional Gaussian filter, G(x, y, t) = (2π)3/2 σx σy σt exp − 2σx 2 + 2σyy 2 + 2σt t 2 ,to the vraw (x, y, t) to remove the noise. The standard deviations of the Gaussian were chosen as σx = σy = 16 pixels(2 PIV grids) and σt = 0.04 s (2 frames), which gave the best results in our experimental
conditions. All the data
RRR
presented in the main text were based on the filtered velocity field, v(x, y, t) =vraw (x′ , y ′ , t′ )G(x − x′ , y − y ′ , t −t′ ) dx′ dy ′ dt′ .
The positions and the radii of the wells were detected by the circular Hough transform using the MATLAB function
imfindcircles. A background image was first obtained by averaging all the 7500 frames of the movie. We applied anadaptive thresholding (contrast limited adaptive histogram equalization, CLAHE) to the background image to cancel
the inhomogeneity of the illumination to ensure the accurate detection of the wells. The circular Hough transformwas then applied to the processed image to detect the positions and the radii of the wells.
E.

Analyzed region and bias in rotational directions

Due to the confinement process, in some regions the bacterial suspensions were not completely confined within thewells and leaked out in between the top and bottom PDMS membranes, which we excluded from the analysis. Since
the well patterns were fabricated only on the bottom PDMS membrane, the leakage breaks the up-down symmetry,
resulting in weak bias in CW rotations. This is the reason for our choice of the green rectangular region used for the
analysis (Fig. S1).
The bias in CW rotation in the leaked regions looking from the above is exemplified in Fig. S3, in which the wellswith the smaller radii on the left side of the green rectangular region in Fig. S1 are also plotted. The bias in therotational directions was quantified by calculating the fraction of time of CW rotation, tCW /(tCW + tCCW ), where

4
1.2
1
0.8
0.6
0.4
0.2
0
-0.2
42

44

46

48

50

FIG. S3. The fraction of time that the spins rotated in the CW direction, tCW /(tCW + tCCW ), as a function of the average
detected radius ⟨Rdetect ⟩ of the wells with the same designed radius Rdesign . Error bars: standard deviations (blue) and standard
errors (red). The wells with the radii < 44 µm showed a significant bias in the CW direction, looking from above, due to theleakage of the bacterial suspension, which was excluded from the analysis.

tCW and tCCW are the durations of CW and CCW rotations judged by the sign of the spin Si (t), respectively. We
calculated this quantity for each well and then calculated the mean value of the wells with the same designed radius
Rdesign . The absence of the CW-CCW bias was used as a criterion for determining the regions with good confinement
to be used in the further analysis.
F.

Correlation time of the spins

Correlation times of the spins were calculated to demonstrate the transition from the stabilized vortex state to theoscillatory state. For each well, we calculate the spin Si (t), where i denotes the index of the well. Then, we calculated
the correlation function,
Ci (τ ) =⟨Si (t)Si (t + τ )⟩t,⟨Si (t)2 ⟩t(S1)where ⟨ ⟩t means temporal averaging. We defined the correlation time of the spin at which Ci (τ ) becomes smaller than
1/e for the first time. Within the experimental observation period of 150 s, some wells did not reverse the sign of Si (t)and exhibited quite stable behavior. Therefore, we could not define the correlation time for such wells and instead
we plotted these wells at the vertical axis denoted as > 150 in Fig. 2(c). Since taking the average is not possible
due to these data points with diverging correlation times, we calculated the moving median with the window widthof 2 µm. To highlight the discontinuity between the diverging data points at > 150 and those with finite correlation
times (< 150), the red line connecting these two regions is plotted with a dashed line in Fig.2(c).
G.

Velocity & vorticity profiles and fit to the analytical solution

Quantitative comparisons between the experimental data and theoretical solutions were performed by examining
the velocity and vorticity profiles for the single stabilized vortex, shown in Fig. 2(a). Due to the PIV interrogation
box size 16×16 pixels and the standard deviation of the Gaussian filter σx = 16 pixels, the obtained velocity andvorticity fields can leak out from the radius R detected by the image analysis by 32 pixels, corresponding to 20.8 µm.
This length scale is consistent with what is observed in Fig. 2(a), where the experimental data extends beyond theradius R. These penetrating fields may justify to the numerical implementation of the damping terms for the velocity
and vorticity fields.

5
The velocity and vorticity profiles in the single stabilized vortex, shown in Fig. 2(a), were fitted by the theoreticalcurves Eq. (S15) and Eq. (S19) for n = 0. The best fits yield G0+ = 1.51 s−1 , G0− = 0.588 s−1 , k0+ = 0.0914 µm−1 ,k0− = 0.0453 µm−1 for v and G0+ = 1.31 s−1 , G0− = 0.633 s−1 , k0+ = 0.0920 µm−1 , k0− = 0.0483 µm−1 for ω.
II.

SUPPLEMENTARY NOTE 2. REGRESSION TO THE TTSHEA.

Least squares method

We inferred the parameters in the TTSHE by regressing the experimental data to the TTSHE. Specifically, weapplied the least squares method to the experimental data of the vorticity field ω(x, y, t) to find the parameters
A, C, Γ0 , Γ2 in the TTSHE. By assuming the vorticity equation for the TTSHE,λ(v · ∇)ω + Aω + C∇ × |v|2 v − Γ0 ∇2 ω + Γ2 ∇4 ω = −∂t ω,(S2)

we constructed its matrix representation by using the experimental data. For a single regression, we used the datafrom time t to t + (N − 1)∆t with ∆t being the time interval between the frames and N being the number of frames
used for the regression. The matrix representation is given by,
ΦP = Ψ,

(S3)

with,


−∂t ωtr1−∂t ωtr2
..
.













 


λ
−∂t ωtrMr1 −∂t ωt+∆t
A


 
r2
P =  C  , Ψ =  −∂t ωt+∆t
,


Γ 
..


0


.
Γ2


r
M

 −∂t ωt+∆t


..




.
rM
−∂t ωt+(N−1)∆t(vtr1 · ∇)ωtr1(vtr2 · ∇)ωtr2
..
.

ωtr1
ωtr2
..
.

∇ × |vtr1 |2 vtr1
∇ × |vtr2 |2 vtr2
..
.

(S4)

−∇2 ωtr1−∇2 ωtr2..
.

∇4 ωtr1∇4 ωtr2
..
.













rM
rM
rM
rM
rM
rM 2 rM2
4
(vt · ∇)ωtωt
∇ × |vt | vt
−∇ ωt
∇ ωt


r1
r1
r1
r1
2 r1
2 r1
4 r1
(vt+∆t · ∇)ωt+∆tωt+∆t∇ × |vt+∆t | vt+∆t−∇ ωt+∆t∇ ωt+∆t 



r2
r2
r2
r2
r2
r2
r2
2
2
4
Φ=
(vt+∆t · ∇)ωt+∆tωt+∆t∇ × |vt+∆t | vt+∆t−∇ ωt+∆t∇ ωt+∆t  ,


..
..
..
..
..




.
.
.
.
.
 r1

rM
r1
r1
2 r1
2 r1
4 r1
(vt+(N −1)∆t · ∇)ωt+(N

ω
∇
×
|v
|
v
−∇
ω
∇
ω
−1)∆t
t+(N −1)∆tt+(N −1)∆tt+(N −1)∆tt+(N −1)∆tt+(N −1)∆t 



..
..
..
..
..


.
.
.
.
.
rM
rM
rM
rM
2 rM
2 rM
4 rM
(vt+(N −1)∆t · ∇)ωt+(N −1)∆t ωt+(N −1)∆t ∇ × |vt+(N −1)∆t | vt+(N −1)∆t −∇ ωt+(N −1)∆t ∇ ωt+(N −1)∆t(S5)where r1 to rM denote the positions within the region of interest (ROI) used for the regression, and vtrk = v(rk , t)and ωtrk = ω(rk , t). Each column of Φ and Ψ contains data from all the M positions within the ROI and all the time
points between t and t + (N − 1)∆t, totaling M N elements. This overdetermined system was solved by the leastsquares method using the pseudo-inverse matrix of Φ defined as,
Φ+ = (ΦT Φ)−1 ΦT ,

(S6)

6
TABLE S1. Estimated parameters of the TTSHE from the regression of the numerical data.
λ
A
C
Γ0
Γ2
Ground truth
9
0.5
1.6
-2
1
8.844 0.4599 1.2033 -1.8204 0.9231
AverageStandard error 0.012 0.0026 0.0046 0.0073 0.0030
Standard deviation 0.054 0.012 0.021 0.034 0.014

where ΦT is the transpose of Φ. Then, the optimal parameters Popt were estimated as,
Popt = Φ+ Ψ.

(S7)Since the TTSHE is a model for the bulk flow, we excluded the region close to the confining boundary from the ROI
so that the boundary does not affect the parameter estimations. The time window N ∆t was determined based on
the spin correlation time. We repeated this regression at time t and then shifted the time window to obtain the timeseries of estimated parameters. In this computation, we used central differences except for ∇2 ω and ∇4 ω that were
calculated by using the fast Fourier transform.
B.

Verification with numerical data set

To verify this method, we first applied this regression to the numerical data obtained from the TTSHE. We usedthe data for R = 5.30, which shows the two oscillating vortices with the period of 117 in the numerical unit. A circular
ROI centered at the well with half of the well’s radius R/2 = 2.65 was chosen to diminish the boundary effects. The
time window was chosen to be N ∆t = 22, which is the correlation time of this oscillatory state. To be consistent withexperimental data, we decreased the spatial resolution of our high-resolution numerical data by 8 folds in both x andy directions, resulting in spatial resolution of ∆x = 0.32, and ∆t was set to 1. These settings result in M = 217 andN = 22, and the duration of the numerical data was 500. To estimate the standard errors of the average values, weused the effective number of independent samples as Nsample = 500/22 = 22.7, which we calculated on the basis of
the time window N ∆t = 22 s and the duration of the data 500.
As shown in Fig. S4 and in Table S1, the time series of the estimated parameters fluctuate in time but theirtemporal averages gave reasonable values. Remarkably, the most accurate estimation was achieved for the mostimportant parameter λ, which determines the strength of the nonlinear advection and the onset of the limit cycle.
The error of the estimated λ was as small as 2%, while those of other parameters were 8% to 25%. Therefore, thisestimation method can reliably be applied at least to estimating λ. The estimations of the other parameters havedifficulties arising from their nonlinearities and higher-order derivatives.
C.

Application to experimental data

Now we apply this regression to our experimental data. We used the velocity and vorticity field data within thewell with oscillatory behavior presented in Fig. 1(e). We excluded the region within 32 pixels (20.8 µm) from theboundary of the well, which composed the ROI for the regression. The choice of 32 pixels was based on the sum of the
PIV interrogation box size 16×16 pixels and the standard deviation of the Gaussian filter σx = 16 pixels. Above thislength scale, the estimated PIV velocity field is not affected by the boundary. With this choice, our ROI contained
M = 107 positions.
The time window N ∆t used for the regression was determined on the basis of the spin correlation time. Since thespin correlation time for this well was 14.2 s, which corresponds to 710 frames, we chose N = 710.
We obtained the time series of estimated parameters shown in Fig. S5. Based on the time window N ∆t = 14.2 s andthe duration of the movie 150 s, we estimate the effective number of independent samples as Nsample = 150/14.2 = 10.6,
which is then used for evaluating the standard errors of the estimated parameters. The average values and thestandard errors of the estimated parameters are shown in Table S2. As expected, the regression gave all the signs of
the estimated parameters consistent with the numerical parameters in the main text.
Although the standard errors of the estimation were relatively large for the nonlinear term C and the fourth-orderderivative term Γ2 , the advection term λ is rather reliably estimated, as expected from our control calculation with
the numerical data set. The value λ = 1.69 ± 0.12, which is larger than unity, is indeed a characteristic of the active
turbulence model of pusher-type microswimmers such as Bacillus subtilis. To compare with the nondimensional form

7
11
10
9
8
7
0

50

100

150

200

250

300

350

400

450

500

0

50

100

150

200

250

300

350

400

450

500

0

50

100

150

200

250

300

350

400

450

500

0

50

100

150

200

250

300

350

400

450

500

0

50

100

150

200

250

300

350

400

450

500

0.6
0.5
0.4

2
1.8
1.6
1.4
1.2
-1.5
-2
-2.5

1.2
1
0.8

FIG. S4. Time series of the parameters in the TTSHE estimated from the regression of the numerical data. Blue curves are
the results of the regression, and the red dashed lines are the ground truth values used for generating the numerical dataset.
The vertical axes in all the plots ranges from -30% to +30% of the ground truth values.

TABLE S2. Estimated parameters in the TTSHE from the regression of the experimental data. The brackets [ ] indicate
the units of the parameters. The time-averaged values have signs consistent with our numerical parameters. Note that λ isnondimensional.
λ A [s−1 ] C [µm−2 · s] Γ0 [µm2 · s−1 ] Γ2 [µm4 · s−1 ]Average
1.69 0.124
0.00011-13.8
2.5
Standard error 0.12 0.047
0.00022
3.9
23Standard deviation 0.38 0.14
0.00068
11.9
71

of the TTSHE, we need to choose typical length L, time T and velocity V . After nondimensionalization, all thecoefficients scale as λ → VLT λ, A → T A, C → V 2 T C, Γ0 → LT2 Γ0 , and Γ2 → LT4 Γ2 . Therefore, if we choose a large
value for V , e.g. based on an ideal unconstrained bacterial turbulence without confinement to the well, the value of
λ can be in the range where we found a limit cycle in our analytical calculation as we describe in the main text.

D.

Regression to other equations

The use of the TTSHE was justified on the basis of the regression to another model, Nikolaevskiy equation [2, 3],also used as a model for the fluid motion of dense bacterial suspensions. Nikolaevskiy equation is, together with the

8
3
2
1
0

50

100

150

50

100

150

0

50

100

150

0

50

100

150

0

50

100

150

0.6
0.4
0.2
0
-0.2
0
10 -3
1
0
-1
-2

0
-20
-40

200
100
0
-100

FIG. S5. Estimated parameters in the TTSHE from the regression of the experimental data as functions of time.

incompressibility, given by,∇ · v = 0,
2

(S8)
4

6

∂t v + (v · ∇)v = −∇p + Γ0 ∇ v − Γ2 ∇ v + Γ4 ∇ v.

(S9)

Γ0 ∇2 ω − Γ2 ∇4 ω + Γ4 ∇6 ω = ∂t ω + (v · ∇)ω,(S10)

RT = ∥ΦPopt − Ψ∥2 ,(S11)We rewrote these equations into a vorticity equation by taking the rotation,and reconstructed a matrix representation for Nikolaevskiy equation for regression, in the same way as we did for the
TTSHE. We obtained the optimal parameters as shown in Table S3, with physically reasonable signs, Γ0 > 0, Γ2 < 0,and Γ4 > 0, for active turbulence with a characteristic length scale.
We evaluated the performance of the two models by comparing the residuals of the regression, RT for the TTSHE
and RN for Nikolaevskiy equation, respectively. Specifically, we calculated,
where ∥ ∥ denotes the norm of a vector. We compute these residuals for each time point t. The same procedure wasapplied also to Nikolaevskiy equation, and the time series of both models are plotted in Fig. S6. As a result, RN isalways larger than RT , suggesting that the TTSHE outperforms Nikolaevskiy equation in terms of their ability todescribe bacterial turbulence, at least, in our setup.
We note that it would be ideal to try modern machine learning techniques such as SINDy [4] to infer the hydrodynamic description for bacterial turbulence by preparing the library of the possible terms. However, this is beyond
the scope of our current work.

9
TABLE S3. Estimated parameters in the Nikolaevskiy equation from the regression of the experimental data. The brackets [ ]indicate the units of the parameters.
Γ0 [µm2 · s−1 ] Γ2 [µm4 · s−1 ] Γ4 [µm6 · s−1 ]Average
11.5-803.72808.8
Standard error
2.6
135.0
662.5
Standard deviation
8.0
417.4
2048.7

10 5

(a)

(b)
10000

1.6

8000
1.4

6000
4000

1.2
1

TTSH

2000

Nikolaevskiy
0

50

100

150

0

0

50

100

150

FIG. S6. Residuals of the regressions. (a) Time series of the residuals for TTSHE RT (blue) and Nikolaevskiy equation RN(red). (b) Time series of the difference of the residuals RN − RT .
III.

SUPPLEMENTARY NOTE 3. DETAILS OF COMPUTATIONAL MODELINGA.

Model with Boundary

To incorporate the effect of boundary, we use the following equation [5, 6]:


2
∂ω
+ λv · ∇ω = aω − b∇ × |v|2 v − 1 + ∇2 ω − γv ∇ × [K(r)v] − γω K(r)ω,
∂t

(S12)

where K(r) is a positive-valued scalar field introduced to damp v, ω outside the well, and γv , γω > 0 are the associateddamping strengths. From its design, K is defined to be ≃ 0 inside the well while ≃ 1 for the outside. Similar to themodel parameters, we also need to determine the values of γv , γω based on experiments. In this paper, we adopt the
choice (λ, a, b, γv , γω ) = (9, 0.5, 1.6, 40, 4) suggested in Ref. [6].
B.

Simulation Scheme

To solve Eq. (S12), we use the pseudospectral method [6], which goes back and forth between the real and Fourier
space to compute the spatial derivatives. In this algorithm, the real space coordinates are defined on a periodic gridr ∈ ∆x · Z2[0,N −1] where N is the number of grid points per dimension and ∆x is the grid spacing. For the Fourier
2

space, wavenumber vector k resides in (π/L) · [−N, N ] where L = N ∆x. Each time we encounter spatial derivatives,
we first Fourier transform quantities using the discrete Fourier transform (DFT), substitute ∇ with ik·, and thenpush the result back to the real space by the inverse DFT (iDFT). The time evolution is done by the Euler method in
the Fourier space since we can exactly handle the linear terms in a stable manner (operator splitting). After updatingω (from time t to t + ∆t), we compute the associated incompressible velocity field v by the streamfunction method,which first solves the Poisson equation ∇2 ψ = −ω and then computes v̂x = iky ψ̂, v̂y = −ikx ψ̂ where ˆ∗ denotes theFourier transform. As this method does not fix the average velocity v̂(0), we need to compute the dynamics of themseparately.

10
In addition to the procedures described above, we need to take care of two simulation-specific issues. First, wecannot compute the full dynamics for all the k because real-space multiplication doubles the range of k. To deal with
that, we manually discard the high-k modes before and after real-space multiplications. Similar treatment is also
applied to K. Specifically, we use the following procedure to generate appropriate K from the naive boolean mask
Knaive (r) = 0 if |r| ≤ R, 1 otherwise:
Algorithm 1 smoothing K
Input: Knaive
1: K̂naive = DFT[Knaive ]2: for all k do

π 2
then
3:
if k2 > 14 · ∆x
4:K̂naive (k) = 0
5:
end if
6: end for
7: K = iDFT[K̂naive ]2
Output: processed K(r)Note that the positivity of K is guaranteed by the square operation. Second, numerical errors can lead to the violation∗of ω(k) = ω(−k) , which is a necessary condition for the real-valued ω. To enforce this condition, we symmetrize ωeach time after the Euler update.
Putting everything together, the simulation proceeds as follows:
Algorithm 2 simulation 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
8:v̂new (0) = exp [(a − 1)∆t] · r √ 1 2 (v + rhs1 ∆t)

2

N

π
else if k2 > 12 · ∆xthen
ω̂new (k) = 0v̂new (k) = 0
else
ω̂new (k) = 21 {ω̂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))

9:
10:
11:
12:
13:

In typical simulations, we set N = 8192, ∆x = 0.005, and ∆t = 0.01. All the computations were performed on
GPUs (NVIDIA RTX A6000 or A100) in single precision, which achieves ∼ 60 times speedup compared to CPU
implementation with standard CPUs (e.g. Intel Xeon W-2295) [7].
IV.

SUPPLEMENTARY NOTE 4. ANALYTICAL THEORYA.

Linear Theory

We start with the TTSHE in the form∂t ω + λv · ∇ω = aω − (1 + ∇2 )2 ω − b∇ × |v|2 v.

(S13)∂t ω = aω − (1 + ∇2 )2 ω.

(S14)Consider a linearized Eq. (S13) in a disk domain of a radius R

11
A generic solution to Eq. (S14) can be written in the form ω =P∞n=−∞ ωn ,ωn = exp(λn t) (Gn+ Jn (kn+ r) + Gn− Jn (kn− r)) exp(inθ),(S15)

where λn is the growth rate, Gn± arbitrary constants, Jn are the Bessel functions, and kn± are given by the equation
q
p
(S16)
kn± = 1 ± a − λn .

First, we calculate the stream function ψ satisfying the condition ∇2 ψ = −ω. In polar coordinates, for eachazimuthal mode n we obtain1n2∂r (r∂r ψn ) − 2 ψn = −ωn .
r
r

(S17)

The solution is
ψn =




Gn+
Gn−
n
exp(inθ) + c.c.
2 Jn (kn+ r) + k 2 Jn (kn− r) + Gn0 r
kn+
n−

(S18)
2

Here, the term Gn0 rn is a solution to the Laplace equation in polar coordinates 1r ∂r (r∂r ψn ) − nr2 ψn = 0.
Correspondingly, we obtain for the velocity components


Gn−
Gn+
(Jn−1 (kn+ kr ) − Jn+1 (kn+ r)) +(Jn−1 (kn− r) − Jn+1 (kn− r)) + nGn0 rn−1 exp(inθ) + c.c.,vθ = −∂r ψn = −
2kn+
2kn−
(S19)


in Gn+Gn−1nexp(inθ) + c.c.
(S20)
vr = ∂θ ψn =2 Jn (kn+ r) + k 2 Jn (kn− r) + Gn0 r
r
r
kn+
n−
Now, to satisfy the b.c., we have the following conditions,
Gn+ Jn (kn+ R) + Gn− Jn (kn− R) = 0,
Gn+
Gn−
n
2 Jn (kn+ R) + k 2 Jn (kn− R) + Gn0 R = 0,
kn+
n−
Gn+
Gn−
(Jn−1 (kn+ R) − Jn+1 (kn+ R)) +(Jn−1 (kn− R) − Jn+1 (kn− R)) + nGn0 Rn−1 = 0.
2kn+
2kn−

(S21)
(S22)
(S23)

We can set Gn+ = 1 due to system linearity. Then, the above 3 equations have 3 unknowns (λn , Gn− , Gn0 ). It
provides the following characteristic equation to determine the growth rates λn for arbitrary n:
Jn−1 (kn+ R)
Jn−1 (kn− R)
2n
2
2
−
−
(1/kn+− 1/kn−)=0
kn+ Jn (kn+ R) kn− Jn (kn− R)
R
B.

(S24)

Orthogonality condition

Equation for the eigenmodes and eigenvaluesL̂ωn = (a − (1 + ∇2n )2 )ωn = λn ωn(S25)

has rather subtle features (here ∇2n = r−1 ∂r r∂r−n2 /r2 is the radial Laplacian). Namely, while the differential operator
is symmetric, it is not formally self-adjoint since the orthogonality condition can not be applied to the vorticity modes
ωn . For the proper self-adjointness, the boundary conditions (b.c.) need to be explicit functions of the vorticity ω
and its derivatives. Because part of the b.c. v = 0 is expressed in terms of the velocity, the partial integration in thescaler product of the eigenfunctions leads to non-vanishing boundary terms. By formulating the equations in termsof the stream function, the b.c can be reformulated in terms of the derivatives of the stream function, which leads tozero boundary terms in the scalar product of the eigenfunctions. If we consider a scalar product,⟨ωnk ωnl ⟩ =
̸ 0,

(S26)

12
the eigenmodes are not orthogonal. Here ωnl,k are different modes corresponding to the same n. It can be seen fromexamining the product ⟨ωnk L̂ωnl ⟩ that generates non-vanishing boundary terms after integration by parts. In thestandard situation, the boundary terms vanish due to the boundary conditions on ω. But we have instead boundary
conditions for the stream function ψ. It implies that the orthogonality should be implemented differently, for the
stream function ψn rather than the vorticity.
Therefore, we formulate the eigenvalue problem in terms of the stream functionL̂1 ψn = (a − (1 + ∇2n )2 )∇2n ψn = λn ∇2n ψn .

(S27)

L̂1 w = (a − (1 + ∇2n )2 )∇2n w = −∂t ω − ∇ × B,

(S28)

The b.c. are ψn = 0, ∂r ψn = 0, ∂r2 ψn = 0 for r = R. With this b.c., the operator L̂1 is self-adjoint, and the boundary
terms vanish identically. A little subtlety here is that it is a generalized eigenvalue problem because instead of λn ψnwe have λn ∇2n ψn .
Now consider the orthogonality condition for each ψn mode. The inhomogeneous equation is written as
where B = −b|v|2 v−λv·∇v. The solvability means that the r.h.s. is orthogonal to the zero eigenmodes ψn . Applying
the solvability condition, after partial integration we obtain
Z R
Z R
∗
dθrψ ∇ × Bdr =rv∗ Bdr,
(S29)
0

0

where v = (vr , vθ ), where vr = ∂θ ψ/r, vθ = −∂r ψ. Correspondingly, for the normalization coefficients, we obtain,
after partial integration,
Z R
Z R
∗
dθrψ ωdr =r|v|2 dr.
(S30)
0

C.

0

Weakly-nonlinear analysis

Here we consider a weakly-nonlinear solutions including only azimuthal modes with n = 0, ±1, ±2,
ψ = C(t)ψ0 (r) + [A1 (t) exp(iθ)ψ1 + A2 (t) exp(2iθ)ψ2 + c.c.] + w.

(S31)

Here ψ0 , ψ1 , ψ2 are the eigenfunctions obtained from linear stability analysis, C(t), A1,2 (t) are slowly-varying functions, and w is small correction to the solution. Functions C(t), A1 (t), A2 (t) are obtained from the corresponding
orthogonality conditions guaranteeing that w does not grow.
Substituting solution Eq. (S31) into Eq. (S13), and retaining terms for each azimuthal harmonics, we obtain after
applying the orthogonality conditions,∂t C = λ0 C − c1 C 3 − c2 C|A1 |2 − c3 C|A2 |2 − 2c4 ReA2 A2∗
1 ,
2
2
2
∂t A1 = λ1 A1 − b1 A1 |A1 | − b2 A1 C − b3 A1 |A2 | − b4 CA2 A∗1 + δ1 A1 C + γ1 A2 A∗1 ,∂t A2 = λ2 A2 − a1 A2 |A2 |2 − a2 A2 C 2 − a3 A2 |A1 |2 − a4 CA21 + δ2 A2 C + γ2 .A21

(S32)
(S33)
(S34)

Here λ0,1,2 are the linear growth rates, and all other coefficients are integrals over the eigenfunctions. In the following,
we can consider A1,2 complex.
D.

Calculations of the nonlinear coefficients

We substitute the approximate solution Eq (S31) into equation (S13) and apply the solvability conditions. The
coefficients m0,1,2 are the integrals of the eigenmode squares
Z R
m0,1,2 =drr|v0,1,2 (r)|2 ,(S35)
0

where the normalized radial eigenmodes are given by Eq. (S18)


Gn+
Gn−
1
n
J
(k
r)
+
J
(k
r)
+
G
r
.
ψn = √
n n+
n n−
n0
2
2
mn kn+
kn−

(S36)

13
The nonlinear terms in Eqs. (S32), (S33), (S34) appear from the quadratic term v·∇ω and the cubic term ∇×|v|2 v.
The coefficients for the quadratic terms δ1,2 , γ1,2 are given by the expressions (after partial integration)Z
λ∗dθdrrv1∗ exp(−iθ)v · ∇v,(S37)
δ1 A1 C + γ1 A2 A1 = −2πZ
λδ2 A2 C + γ2 A21 = −dθdrrv2∗ exp(−2iθ)v · ∇v.
(S38)
2π
(S39)
Here vn = (vrn , vθn ) = (inψn /r, −∂r ψn ). Note that in polar coordinates,v · ∇v = (vr ∂r vr + vθ ∂θ vr − vθ2 /r, vr ∂r vθ + vθ ∂θ vθ + vθ vr /r).
The advection term does not generate any contribution to Eq. (S32) due to the symmetry.
Correspondingly, for a1,2,3,4 , b1,2,3,4 , c1,2,3,4 we obtainZbdθdrrv0∗ v|v|2 ,c1 C 3 + c2 C|A1 |2 + c3 C|A2 |2 + 2c4 ReA2 A2∗
=
1
2π
Z
b
2
2
2
∗
b1 A1 |A1 | + b2 A1 C + b3 A1 |A2 | + b4 CA2 A1 =dθdrrv1∗ exp(−iθ)v|v|2 ,2πZba1 A2 |A2 |2 + a2 A2 C 2 + a3 A2 |A1 |2 + a4 CA21 =dθdrrv2∗ exp(−2iθ)v|v|2 .
2π

(S40)

(S41)
(S42)
(S43)

All linear and nonlinear coefficients are calculated in Mathematica and imported directly into the normal form
equations (S32), (S33), (S34). It appears that γ1 = γ2 , and c4 = a4 = b4 /2. The resulting equations generate a limit
cycle without any further adjustments.
E.

Bridging theory and numerics

To perform faithful comparison, mode decomposition of numerical data was done by spline interpolation on theradial grid followed by the trapezoidal quadrature based on the same formula as in the analytical theory. The upperbound of r integral was set larger than R (typically ≃ 9) to contain the well plus the leakage. The azimuthal modes(Fig. 3(c), Fig. 4(c,d)) were based on the following formulae:
Z
1
v∗n =
dθ e−inθ v∗ (r, θ),
(S44)
√
2π mn
where ∗ = r, θ. Similarly the associated mode amplitudes (Fig. 3(b), Fig. 4(b)) were calculated by the following:sR
2
dr r|vθ0,num. |(S45)|C| = R2,
dr r|vθ0,theory |sR
2
dr r|vθ1,num. ||A1 | = R
(S46)
2,
dr r|vθ1,theory |sR
2
dr r|vθ2,num. |(S47)|A2 | = R2.
dr r|vθ2,theory |R
As for the sign of C, as vθ0,theory > 0, we can naturally identify it as the sign of dr vθ0,num. . Moreover,
RR beyond thesign consistency,
thespinvariabledefined
in
Eq.
(1)
is
proportional
toC,as
readilyverified:S∝r dr dθ rvθ =
i
h
i
R
R
R
P
2
inθ
2
dr r dθ Cvθ0 + n≥1 An vθn e + c.c. = C dr r vθ0 ∝ C, where An are the mode amplitudes and theintegrals run over the area of the i-th well. This proportionality was used in Fig. 3(b). With the help of mode
decomposition and Reff plotted in Fig. S7(a), we can almost directly compare the theory and simulation for somequantities, such as max |ω| shown in Fig. S7(b).
In the large-R region, where we cannot compute reliable estimates of Reff due to oscillation, we still see good
agreements between the theory and simulation as in Fig. S7(c,d).

14

(b)

(a)

0.7

5.25

0.6

max |ω|

Reff

5.00

0.5

4.75

0.4

4.50

0.3

3.75 4.00 4.25 4.50 4.75 5.00 5.25

4.25

R
(c)

4.75

1.00

0.75

1.00

0.2

0.50

|A2|

0.50

0.1

0.0

0.4

0.8

0.0

0.6

0.8

0.4

0.6

0.2

0.4

0.0

0.2

−0.2

0.0
−0.4

−0.2

0.0

0.2

0.00
−0.50

0.2

0.3

0.4

0.5

−0.4

0.4

C

−0.2

0.0

C

−0.25

0.00

0.4
0.25

C

0.2

0.4

0.25

0.6

|A1|

|A1|

C

0.2

0.25

0.2
0.50 0.0

|A1|

0.00

0.8

0.75

0.6

0.50

0.25

0.2

0.00

0.0

0.2

0.4

0.6

0.8

|A1|

1.00

|A2|

−0.2

0.45
0.40
0.35
0.30
0.25
|A1|

|A1|

0.1

|A2|

0.75

|A2|

|A2|

0.2

5.00

Reff

theory
num.

0.3

0.3

|A2|

4.50

(d)
theory
num.

0.0
−0.4

theory
num.
5.25

0.2

4.25

0.4

−0.4

−0.2

0.0

0.2

0.4

C

0.0
−0.4

−0.2

0.0

0.2

0.4

C

FIG. S7. (a) Radius R used for K (damping mask) and effective radius Reff estimated from the simulation. Generally we have
Reff > R due to leakage. (b) R-dependence of max |ω| near the lower bound radius having λ0 > 0. (c,d) Three-dimensional
phase space plot and projections, in the two-vortex state ((c), R = 5.35 for simulation and R = 5.9 for theory) and four-vortex
state ((d), R = 6.4 for simulation and R = 7.0 for theory).

V.

SUPPLEMENTARY MOVIE CAPTIONS

Movie 1: Movie of the whole field of view (2.1 mm × 2.1 mm), corresponding to Fig. S1. The initial 20 seconds ofthe data are presented at twice the real speed (2×). The spatial resolution is lowered to reduce the file size.
Movie 2: The vorticity and velocity fields of the single stabilized vortex (R = 44.6 µm). The color bar of the vorticity
field is the same as in Fig. 1. The movie is played at real-time speed. The scale bar represents 20 µm.
Movie 3: The vorticity and velocity fields of the reversing vortices (R = 46.7 µm). The color bar of the vorticity
field is the same as in Fig. 1. The movie is played at real-time speed. The scale bar represents 20 µm.
Movie 4: The vorticity and velocity fields of the four-vortex state (R = 48.8 µm). The color bar of the vorticity fieldis the same as in Fig. 1. The movie is played at real-time speed. The scale bar represents 20 µm.
Movie 5: The velocity field of the single stabilized vortex (R = 44.6 µm). The velocity vectors are colored based on
the sign of the local vorticity: yellow for positive and green for negative. The movie is played at 4× real-time
speed. The scale bar represents 20 µm.

15Movie 6: The velocity field of the reversing vortices (R = 46.7 µm). The velocity vectors are colored based on thesign of the local vorticity: yellow for positive and green for negative. The movie is played at 4× real-time speed.
The scale bar represents 20 µm.
Movie 7: The velocity field of the four-vortex state (R = 48.8 µm). The velocity vectors are colored based on thesign of the local vorticity: yellow for positive and green for negative. The movie is played at 4× real-time speed.
The scale bar represents 20 µm.
Movie 8: Vorticity field of the single stabilized vortex (Fig. 3(a), left), numerically computed at R = 5.2.
Movie 9: Vorticity field of the reversing two-vortex state (Fig. 3(a), middle), numerically computed at R = 5.6.
Movie 10: Vorticity field of the pulsating four-vortex state (Fig. 3(a), right), numerically computed at R = 6.4.
Movie 11: Vorticity field of the reversing two-vortex state(Fig. 3(b)), numerically computed at R = 5.35.
Movie 12: Azimuthal Fourier decomposition of the vorticity field of the reversing two-vortex state at R = 5.6(Fig. 3(c)). The original vorticity field and the modes n = 0, 1, 2 are shown.
Movie 13: Vorticity field of the pulsating four-vortex state (Fig. 3(d)), numerically computed at R = 6.2.
Movie 14: Vorticity field and its azimuthal Fourier decomposition of the analytically calculated reversing two-vortex
state at R = 5.9 (Fig. 4(e)). The definitions of the modes are the same as in Supplementary Movie 12 andFig. 3(c).
Movie 15: Vorticity field and its azimuthal Fourier decomposition of the analytically calculated pulsating fourvortex state at R = 7.0 (Fig. 4(f)). The definitions of the modes are the same as in Supplementary Movie 12
and Fig. 3(c).

[1] E. Stamhuis and W. Thielicke, PIVlab – Towards User-friendly, Affordable and Accurate Digital Particle Image Velocimetryin MATLAB, Journal of Open Research Software 2, 10.5334/jors.bl (2014).
[2] I. A. Beresnev and V. N. Nikolaevskiy, A model for nonlinear seismic waves in a medium with instability, Physica D:
Nonlinear Phenomena 66, 1 (1993).
[3] J. Słomka and J. Dunkel, Spontaneous mirror-symmetry breaking induces inverse energy cascade in 3D active fluids,
Proceedings of the National Academy of Sciences 114, 2119 (2017).
[4] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear
dynamical systems, Proceedings of the National Academy of Sciences 113, 3932 (2016).
[5] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Meso-scale turbulencein living fluids, Proceedings of the National Academy of Sciences of the United States of America 109, 14308 (2012).
[6] H. Reinken, D. Nishiguchi, S. Heidenreich, A. Sokolov, M. Bär, S. H. Klapp, and I. S. Aranson, Organizing bacterial vortex
lattices by periodic obstacle arrays, Communications Physics 3, 76 (2020).
[7] S. Shiratani, K. A. Takeuchi, and D. Nishiguchi, Route to turbulence via oscillatory states in polar active fluid underconfinement, arXiv preprint arXiv:2304.03306 (2023).