# Scaling transition of active turbulence from two to three dimensions

Da Wei  Beijing National Laboratory for Condensed Matter Physics, Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China  Yaochen Yang  CAS
Key Laboratory for Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China  School of Physical
Sciences, University of Chinese Academy of Sciences, 19A Yuquan Road, Beijing
100049, China  Xuefeng Wei  CAS Key Laboratory for Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  School of Physical Sciences, University of Chinese Academy of Sciences,
19A Yuquan Road, Beijing 100049, China  Wenzhou Institute, University of
Chinese Academy of Sciences, Wenzhou, Zhejiang 325000, China  Ramin
Golestanian  Max Planck Institute for Dynamics and Self-Organization (MPIDS),
D-37077 Göttingen, Germany  Rudolf Peierls center for Theoretical Physics,
University of Oxford, Oxford OX1 3PU, United Kingdom  Ming Li  Beijing
National Laboratory for Condensed Matter Physics, Institute of Physics,
Chinese Academy of Sciences, Beijing 100190, China  Songshan Lake Materials
Laboratory, Dongguan, Guangdong 523808, China  Fanlong Meng
<EMAIL_ADDRESS>CAS Key Laboratory for Theoretical Physics, Institute
of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, 19A
Yuquan Road, Beijing 100049, China  Wenzhou Institute, University of Chinese
Academy of Sciences, Wenzhou, Zhejiang 325000, China  Yi Peng
<EMAIL_ADDRESS>Beijing National Laboratory for Condensed Matter Physics,
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, 19A
Yuquan Road, Beijing 100049, China

###### Abstract

Turbulent flows are observed in low-Reynolds active fluids. They are
intrinsically different from the classical inertial turbulence and behave
distinctively in two- and three-dimensions. Understanding the behaviors of
this new type of turbulence and their dependence on the system dimensionality
is a fundamental challenge in non-equilibrium physics. We experimentally
measure flow structures and energy spectra of bacterial turbulence between two
large parallel plates spaced by different heights $H$. The turbulence exhibits
three regimes as H increases, resulting from the competition of bacterial
length, vortex size and H. This is marked by two critical heights ($H_{0}$ and
$H_{1}$) and a $H^{0.5}$ scaling law of vortex size in the large-$H$ limit.
Meanwhile, the spectra display distinct universal scaling laws in quasi-two-
dimensional (2D) and three-dimensional (3D) regimes, independent of bacterial
activity, length and $H$, whereas scaling exponents exhibit transitions in the
crossover. To understand the scaling laws, we develop a hydrodynamic model
using image systems to represent the effect of no-slip confining boundaries.
This model predicts universal 1 and -4 scaling on large and small length
scales, respectively, and -2 and -1 on intermediate length scales in 2D and
3D, respectively, which are consistent with the experimental results. Our
study suggests a framework for investigating the effect of dimensionality on
non-equilibrium self-organized systems.

###### pacs:

Turbulence is ubiquitous in nature, from galaxy formation, ocean current to
human breathing. Besides the classical turbulence at high Reynolds number
(Re), turbulent flows are also discovered in active fluids at low Re, such as
bacterial suspensions [1, 2, 3], sperm swarms [4], mixtures of microtubules
and motor proteins [5, 6, 7], epithelial cells [8, 9], and artificial motile
colloids [10]. These active fluids are constituted of a large number of self-
driven agents that inject kinetic energy individually into the systems. These
agents generate active stress and drive the formation of turbulence. The
internal-driven and self-organized nature distinguish this new class of
turbulence from the classical ones. In this light, they are referred to as
active turbulence [11].

Classical inertial turbulence exhibits universal scaling in kinetic energy
spectra. In three dimensions (3D), turbulent flows are driven externally on
macroscopic length scales and kinetic energy cascades towards small scales
where it dissipates. On intermediate scales, vortices exhibit a scale-
invariant structure, giving rise to a universal scaling exponent of -5/3 [12].
When the systems are confined to be quasi-two-dimensional (2D), like large-
scale motions in the atmosphere and oceans, dual cascade emerges with varied
scaling because the vortex stretching is suppressed [13, 14].

Figure 1: Bacterial turbulence in samples of increasing confinement size (a)
Schematics of bacterial suspension confined between two parallel walls spaced
by $H$. Bright field microscopy is performed at a distance $d$ from the bottom
wall. (b-c) Typical flow patterns sampled at $d_{c}=H/2$ for different $H$.
(b) The time-elapsed images. The elapsing time of each panel corresponds to
$\sim 5$µm/$\langle v\rangle$, where $\langle v\rangle$ is the mean flow field
speed. (c) Corresponding vorticity field and streamlines. $\tilde{\omega}$
represents the vorticity scaled by its mean absolute value over the field
($\langle|\omega|\rangle$). Scale bars: 20 µm.

Compared with the better understood classical turbulence, whether universal
scaling also exists in active turbulence is still under debate [15, 16, 5, 17,
18]. Phenomenological models [15, 19, 20] extending Toner-Tu equations for
polar flocking predict non-universal scaling laws in energy spectra [16]. On
the contrary, theories considering the system as active liquid crystals
suggest universal scaling [21, 17] for 2D active nematic turbulence, which is
experimentally confirmed in the quasi-2D regime [18]. Evidently,
dimensionality and confinement play crucial roles in shaping the behaviors of
active turbulence. So far, it has been shown that active turbulence bears
different statistics in thin and thick samples [15]; and confinement
facilitates the coherent collective motion of active agents [22, 23, 24, 6,
25, 26, 27]. Despite the consensus on the crucial role of dimensionality and
confinement, distinguishing the 2D, 3D, and the crossover regimes remains to
be a challenging task - let alone resolving the flow structures and scaling
behaviors in these regimes.

Here we combine experiments and theories to resolve how active turbulence
evolves from 2D to 3D and whether it exhibits universal statistical
properties. We experimentally characterize turbulent flows in bacterial
suspensions by measuring spatial velocity correlations, structure functions,
and kinetic energy spectra with varying confinement heights. As a result, we
uncover two critical heights that classify the turbulence into three regimes:
the 2D, 3D, and an intermediate one. The two heights emerge from the
competition between confinement and the cell size, and that between the
confinement and the vortex size. Our results show that bacterial turbulence
follows universal scaling laws in the 2D and 3D limits, independent of cell
activity and length, and the system height. A hydrodynamic model is developed
to understand the universal scaling in the 2D and 3D limits, as well as the
transition between them, which is consistent with our experiments.

## The structures of active turbulence

We employ Escherichia coli as our model system to study active turbulence (see
Methods). Bacteria are suspended in a minimum motility buffer in which
bacteria cannot grow further. To examine the effect of dimensionality, we
inject bacterial suspensions into chambers confined in $z-$direction with
different heights $H$ (2-400 µm, dimensions in $xy$-plane are 5-10 mm, Fig.
6a). The chamber is sealed to avoid external flows and oxygen gradient. The
bacteria concentration $\rho=3.2\times 10^{10}$ cell/ml, which corresponds to
bacterial volume fraction $\phi\approx 6\%$. At this concentration, bacteria
interact with their neighbours mainly by the hydrodynamic interactions [2]. We
use video microscopy to image bacterial flow in $xy$-plane at a distance ($d$)
from the bottom of the chamber (Fig. 1a). We perform particle image
velocimetry (PIV) to extract the velocity field $\textbf{v}(\textbf{r})$ and
compute the vorticity field $\omega(\textbf{r})$ (Fig. 1b-c, see also
Methods).

Figure 2: The velocity structures of bacterial turbulence. (a) $L_{vv}$ as a
function of $H$. Data for $H=400$µm (hollow symbol) are experimentally
estimated, see Methods. Inset: representative $C_{vv}$ as a function of
$\delta R$; and $L_{vv}$ is defined as where $C_{vv}(\delta R)=1/e$. Dashed
line: best fit with $L_{vv}=\sqrt{a(H-H_{0})}$ [26]. (b)
$\tilde{L}_{vv}=L_{vv}/\sqrt{a(H-H_{0})}$ as a function of $\tilde{d}=2d/H$.
Left inset: measured in samples with $H\leq 40~{}$µm; right inset: with
$H>40~{}$µm. (c) Effective vortex diameter $D_{v}$ as a function of $H$.
Inset: typical $S^{2}_{\perp}$; and $D_{v}$ marks the peak position of
$S^{2}_{\perp}$. Symbols and error bars in this figure represent
mean$\pm$standard deviation. Data in (a) and (c) are taken at $d_{c}=H/2$.

Bacterial turbulence exhibits distinctive flow patterns from 2D to 3D (Fig. 1,
see also Supplementary Mov. 1-3). To characterize the evolving patterns, we
first calculate the spatial velocity correlation functions, $C_{vv}(\delta
R)=\langle\textbf{v}(\textbf{r})\cdot\langle\textbf{v}(\textbf{r}+\delta\textbf{R})\rangle/\langle\textbf{v}(\textbf{r})\cdot\textbf{v}(\textbf{r})\rangle$,
where $\delta R$ is the spatial distance and the angular brackets
$\langle\cdot\rangle$ represents averaging over space and time. Representative
$C_{vv}$ measured at the central height $d_{c}=H/2$ in samples of different
$H$ are displayed in Fig. 2a inset. The correlation length $L_{vv}$ is
measured as the distance where $C_{vv}$ drops to $1/e$ [28, 18]. For a
particular $H$, $L_{vv}$ is robust against variations in bacterial activity
after the onset of turbulence [29, 20], see also Fig. 9. $L_{vv}$ increases
negligibly when $H$ is below a critical height $H_{0}$, whereas it soars up as
$L_{vv}=\sqrt{a(H-H_{0})}$ above $H_{0}$ (Fig. 2a; $a=12.3~{}$µm,
$H_{0}=9.6~{}$µm). Previous experimental measurements agree with our results
[19, 30, 31], which yet spans over a wider range of the system height. This
allows us to resolve the critical length $H_{0}$ in the 2D samples, and to
discover the scaling $L_{vv}\sim H^{0.5}$ in the 3D limit - which does not
saturate up to $H=400~{}$µm (Methods). To understand this scaling, we conduct
stability analysis on active polar systems confined by two no-slip walls (see
Supplementary Information). We find that the confining boundaries make the
system deviate from the long-wavelength instability predicted by the kinetic
theory [32, 33]. The wavelength corresponding to the fastest growth rate
increases as $H^{0.5}$ in thick systems, which reflects the characteristic
lengths of bacterial turbulence and underlies the observed trend in
$L_{vv}(H)$.

At a certain $H$, $L_{vv}$ varies symmetrically with $d$ (Fig. 1a) with
respect to $d_{c}=H/2$ (Fig. 2 and Fig. 8). We make $L_{vv}$ and $d$
dimensionless: $\tilde{L}_{vv}=L_{vv}/\sqrt{a(H-H_{0})}$ and $\tilde{d}=2d/H$
(Fig. 2b). When $H$ is slightly above $H_{0}$, $\tilde{L}_{vv}$ depends weakly
on $\tilde{d}$, e.g., when $H=20$ µm $\tilde{L}_{vv}$ varies less than $10\%$
over the entire $z$ domain (Fig. 2b inset). As $H$ increases, the dependence
of $\tilde{L}_{vv}$ on $\tilde{d}$ becomes stronger (Fig. 2b inset).
Interestingly, above a critical height $H_{1}\approx 40~{}$µm,
$\tilde{L}_{vv}$($\tilde{d}$) collapses to a master curve
$\tilde{L}^{2}_{vv}+\tilde{d}^{\,2}=1$, which is $H$-independent (Fig. 2b
insets). Such universal relation reveals scale-invariant turbulence structures
for $H\geq H_{1}$.

We employ velocity structure functions to further resolve the structure of
bacterial turbulence [34, 14, 15]. They average the two-point velocity
difference
$\delta\textbf{u}=\textbf{u}(\textbf{r}+\delta\textbf{R})-\textbf{u}(0)$ on
orthogonal directions of $\delta\textbf{R}$ over the flow domain:
$S_{i}^{n}=\langle[\delta\textbf{u}\cdot\textbf{e}_{i}]^{n}\rangle\ \
(i=\perp,{\mkern 3.0mu\vphantom{\perp}\vrule depth=0.0pt\mkern 2.0mu\vrule
depth=0.0pt\mkern 3.0mu})$, where $n$ is the order and
$\textbf{e}_{\perp,{\mkern 3.0mu\vphantom{\perp}\vrule depth=0.0pt\mkern
2.0mu\vrule depth=0.0pt\mkern 3.0mu}}$ are the unit vectors normal and
parallel to $\delta\textbf{R}$ respectively. The maxima of $S_{\perp}^{2}$
correspond to an effective vortex size $D_{v}$ (Fig. 2c inset) [34]. $D_{v}$
stays around 10 µm for $H\leq H_{0}$ and bears similar $\sim H^{0.5}$ scaling
above $H_{0}$ (Fig. 2c). Importantly, $D_{v}(H)$ intersects with $H=D_{v}$ at
$H_{1}$ (Fig. 2c), suggesting that $H_{1}$ emerges from the interplay between
the vortex size and the confinement, e.g. due to the suppression of $z$-flows
by the walls when $H<D_{v}$. The near-wall hydrodynamics is long known to lead
rod-shaped bacteria to swimming parallel to the wall [35], and confinement
smaller than vortex sizes rectify bacterial flow from turbulent to coherent
(i.e. with uniform direction) [22, 23]. To identify this, we measure flow
velocity in all three dimensions simultaneously near the focal plane with an
enhanced particle tracking velocimetry (Methods and Fig. 11). A strong
suppression of $z$-flow is confirmed in thin samples, revealed by the large
ratio of the $xy$-velocity to the velocity in $z$, $\langle
u_{xy}/u_{z}\rangle$. The suppression is relaxed after $H\geq 40$µm $\approx
H_{1}$ \- as $\langle u_{xy}/u_{z}\rangle$ has reduced to the isotropic value
$\sqrt{2}$ (Fig. 11). Altogether, by analyzing velocity correlations and
structures of active turbulence in real space, we classify the turbulent flows
into three regimes: 2D ($H\leq H_{0}$), 2D-3D crossover ($H_{0}<H\leq H_{1}$),
and 3D ($H>H_{1}$).

## Kinetic energy spectra and scalings

Figure 3: Kinematic energy spectra of bacterial turbulence in different
regimes. Energy spectra for different activities ($\langle v\rangle$) at
representative heights: $H\approx H_{0}$ (a), $H\approx H_{1}$ (b), and
$H>H_{1}$ (c). $\langle v\rangle_{\rm max}$ is highest mean speed observed
over multiple samples of the same height. The shown wavenumbers are scaled
with $k_{b}=2\pi/L_{b0}$, where $L_{b0}=3$ µm is the bacterial length. Energy
spectra evolving with $H$ for $H\lesssim H_{0}$ (d), $H_{0}\lesssim H\lesssim
H_{1}$ (e), and $H>H_{1}$ (f). Every spectrum is averaged over the activity
range $\langle v\rangle/\langle v\rangle_{\rm max}>0.5$. (g) Characteristic
lengths extracted. $L^{*}=2\pi/k^{*}$, with $k^{*}$ defined as shown in the
schematics in f right inset. Left inset of (f): spectra vertically shifted for
better visualization and with $k^{*}$ marked marked as diamonds. (h) Scaling
exponents $\gamma_{1}$ and $\gamma_{2}$ as functions of $H$. Dotted and the
dash-dotted lines: theoretical prediction of Eq. (3). Blue, yellow, and red
shading mark data obtained from 2D, 2D-3D crossover, and 3D regimes,
respectively.

We next study how kinetic energy of active turbulence is distributed on
different length scales, by measuring the kinetic energy spectrum $E(k)$
defined as $\langle\textbf{u}^{2}\rangle/2=\int E(k)dk$ in the
wavenumber($k$)-space. $E(k)$ is reported for each height over multiple
samples at varying bacterial activities. At a specific $H$, the dimensionless
energy spectra $E(k)/\frac{1}{2}\langle u\rangle^{2}$ are independent of the
activity (Fig. 3a-c), showing the same scaling behaviors. Therefore, we report
the dimensionless spectrum averaged over samples and activities for a certain
$H$ (Fig. 3d-f).

In 2D ($H\leq H_{0}$), the energy spectrum is featured by two scaling
regimes:$E(k)\sim k$ at low-$k$ (large length scales) and $\sim k^{-2}$ at
high-$k$ (small length scales) (Fig. 3a, d). Note that the rising tail at
$k\gtrsim 0.3k_{b}$ results from noise in velocimetry and microscopy in
thicker samples. These spectra shift to lower $k$ (Fig. 3d-f) as $H$
increases, which can be measured by the turning point $k^{*}$ between two
regimes of different scaling exponents (Fig. 3f right inset). The
corresponding lengths $L^{*}=2\pi/k^{*}$ are displayed in Fig. 3g. They both
evolves as $H^{0.5}$ in the 3D limit, same as $L_{vv}$ and $D_{v}$. All these
three quantities reflect the averaged length scales of unstable modes in
bacterial turbulence at varying $H$, yet in different definitions or averaging
methods.

The scaling laws transition from 2D to 3D in two steps: the high-$k$ scaling
exponent $\gamma_{2}$ drops first from -2 to -4 at $H_{0}$, and then the
low-$k$ exponent $\gamma_{1}$ changes from +1 to -1 when $H>H_{1}$ (Figure
3h). The previous measurements on energy spectra are well in line with our
results [15, 31]. Altogether, from the kinetic energy spectra we obtain the
same statistical properties of active turbulence: the $H^{0.5}$ scaling of
vortex size (Fig. 3g and Fig. 2a) and the two critical heights $H_{0}$ and
$H_{1}$ (Fig. 3h and Fig. 2a-b).

Figure 4: Effect of bacterial length on the turbulence (a) The distributions
of cell lengths under different concentrations of cephalexin [Cephalexin].
Dotted and solid short lines respectively mark the median and mean values.
Inset: mean value of $L_{b}$ for different [Cephalexin]. Error bars represent
standard deviations. (b) $L_{vv}$ as a function of $H$ for different $L_{b}$.
Inset: zoom in of the marked region. (c) $D_{v}$ as a function of $H$ for
different $L_{b}$. (d-e) Energy spectra obtained with the most elongated cells
($3.6L_{b0}$). The blue-yellow transition corresponds to $H_{0}$. (f) Spectral
scaling $\gamma_{1}$ and $\gamma_{2}$ obtained with elongated bacteria. Black
dotted and dash-dotted lines: model predictions. Inset: $H_{0}$ and $H_{1}$
for different $L_{b}$. Error bars show standard deviations of $L_{b}$
($x$-axis) and uncertainties in measuring $H$ ($y$-axis). All measurements are
obtained at $d_{c}=H/2$ with $\phi=0.06$.

## The effect of cell length on active turbulence

Besides the system height and the vortex size, the cell length $L_{b}$ is
another fundamental length scale of the active turbulence system. To clarify
how $L_{b}$ affects the turbulent behaviors, we vary bacterial length by
adding cephalexin, an antibiotics that prohibits cell division [28, 36]
(Methods). $L_{b}$ increases monotonically with the cephalexin concentration
[Cephalexin], and meanwhile no obvious changes are found in the size
distributions (Fig. 4a and Fig. 7) and the motility of a bacterium [28, 36].
We experiment with bacteria whose mean length $\langle L_{b}\rangle\approx
2.0L_{b0}$ (5.8 µm, [Cephalexin]=30 µg/ml) and $3.6L_{b0}$ (10.7 µm,
[Cephalexin]=$50~{}$µg/ml) at the same bacterial volume fraction (6%).

Figure 4 shows the correlation length, vortex size, and energy spectra for
different bacterial lengths. Varying $L_{b}$ does not alter the $H^{0.5}$
scaling of $D_{v}$ or $L_{vv}$ (Fig. 4b-c). Energy spectra measured with
elongated cells show the same scaling exponents in low-$k$ and high-$k$
regimes, with $(\gamma_{1},\gamma_{2})=(+1,-2)$ in 2D and (-1, -4) in 3D,
respectively (Fig. 4d-e). This supports the universality of the scaling laws
in bacterial turbulence. However, the critical heights increase for longer
bacteria. $H_{0}$ follows an empirical relation $H_{0}=2.1L_{b}+2.3~{}$µm and
$H_{1}$ increases similarly with $L_{b}$ (Fig. 4b, f). Additionally, $H_{0}$
and $H_{1}$ measured in the representation of velocity structures ($L_{vv}$
and $D_{v}$) and in the representation of energy spectra display excellent
equivalence (Fig. 10).

## Hydrodynamic theory and universal scaling

To understand the scaling behaviors of 2D and 3D active turbulence, we develop
a hydrodynamic model for active fluids confined by two no-slip boundaries. The
model considers the active fluid as a continuum force field, and represents
the effects of the upper and the lower walls by two image systems. A Green
function with the first two image reflections is employed to approximate the
images for a stokeslet with two no-slip boundaries [37, 38, 39], see Fig. 5a
and Supplementary Information for more details. The energy spectrum $E(k)$
written in the form of the tensor is:

$E(k)\sim
k\langle|\hat{\textbf{u}}(\textbf{k})|^{2}\rangle=k\bar{G}_{\alpha\beta}\bar{G}_{\alpha\gamma}\langle
F_{\beta}F^{*}_{\gamma}\rangle,$ (1)

where $\bar{G}_{\alpha\beta}$ ($\alpha,\beta,\gamma=x,y$) is the Green
function in Fourier representation, and $F$ is the active force density. The
asymptotic behavior of $\bar{G}_{\alpha\beta}$ is:

$\bar{G}_{\alpha\beta}\sim\begin{dcases}(\delta_{\alpha\beta}-\frac{k_{\alpha}k_{\beta}}{2k^{2}})\cdot\frac{-H}{k},&kH\ll
1;\\\
(\delta_{\alpha\beta}-\frac{k_{\alpha}k_{\beta}}{k^{2}})\cdot\frac{2}{k^{2}},&kH\gg
1.\end{dcases}$ (2)

Hence, correspondingly, $E(k)$ behaves as:

$E(k)\sim\begin{dcases}\frac{3H^{2}\eta^{2}k}{4}\langle\omega^{2}(\textbf{k})\rangle+\frac{H^{2}}{4k}\langle
F^{2}(\textbf{k})\rangle,&kH\ll 1;\\\
\frac{\eta^{2}}{k}\langle\omega^{2}(\textbf{k})\rangle,&kH\gg 1.\end{dcases}$
(3)

Here $\eta$ is the shear viscosity and $\omega(\textbf{k})$ the vorticity
field in Fourier space. Furthermore, the scaling behavior of
$\langle\omega^{2}(\textbf{k})\rangle$ is governed by the competition between
$k$ and the wavenumber corresponding to the vortex size ($D_{v}^{-1}$). Under
the assumptions that the active vortices are uncorrelated and exponentially
distributed in size (Fig. 13), $\langle\omega^{2}(\textbf{k})\rangle$ scales
respectively as $k^{0}$ and $k^{-3}$ for $kD_{v}\ll 1$ and $kD_{v}\gg 1$ [21].
Altogether, in the limit $F\ll\eta k\omega$, the interplay among $k$,
$H^{-1}$, and $D_{v}^{-1}$, give rise to four different scaling regimes, see
Fig. 5b. We implement the empirical relation between $D_{v}$ and $H$ from the
experiments (Fig. 2c) and display the predicted scaling regimes for our system
in Fig. 5c.

The $k^{1}$\- and $k^{-4}$-scaling laws are revealed to be universal, and they
govern the kinetic energy distribution on the large and small scales,
respectively. On large scales ($k\ll H^{-1},D_{v}^{-1}$), the growth of long-
wavelength modes are suppressed by the presence of two no-slip boundaries, and
such size selection underpins the $k^{1}$-scaling. On the small scales ($k\gg
H^{-1},D_{v}^{-1}$), the $k^{-4}$-scaling marks the energy dissipation within
vortices, which is dominated by viscosity [11]. Meanwhile, $k^{-2}$ and
$k^{-1}$ represent two mutually exclusive scaling laws, which govern the
transitional regimes in the $k$-space for 2D and 3D samples respectively. The
$k^{-2}$-scaling governs the regime $D_{v}^{-1}\ll k\ll H^{-1}$ and becomes
evident in the 2D limit ($H^{-1}\gg k,D_{v}^{-1}$). On the other hand, the
$k^{-1}$-scaling governs $H^{-1}\ll k\ll D_{v}^{-1}$ and emerges in the 3D
limit ($H^{-1}\ll k,D_{v}^{-1}$). In all, our model supports the existence of
universal scaling for active turbulence in both the 2D and 3D limits.

Consistently, our experimental results agree with the predicted universal
scaling. In 3D samples, the measured scaling $(\gamma_{1},\
\gamma_{2})=(-1,-4)$ (Fig. 3c), and $k^{1}$-scaling are observed at small
wavenumbers (Fig. 12). These all align perfectly with the prediction, see the
right vertical dashed line in Fig. 5c. For samples in the 2D-3D crossover
regime, e.g., around $H=H_{1}$, the kinetic energy spectrum is predicted to
evolve from $k^{1}$-scaling to $k^{-4}$-scaling with no transitional scaling
(middle vertical line in Fig. 5c). It is exactly what we observe
$(\gamma_{1},\ \gamma_{2})=(1,-4)$ around $H=H_{1}$, see Fig. 3b and e.
Lastly, in the 2D limit, only the $k^{1}$\- and $k^{-2}$-scaling are observed
(Fig. 3a) because the predicted $k^{-4}$-scaling only exists on the scales
smaller than the bacterial size. On such scales, the bacterial suspensions
cannot be regarded as a continuum and thus the model becomes not applicable,
as is shown by the gray shaded area in Fig. 5c.

Figure 5: Scaling regimes predicted by the hydrodynamic model. (a) Scheme of
the hydrodynamic model, where boundary effects are approximated by two images.
(b) Asymptotic scaling behaviors results from the competition among $k$,
$H^{-1}$, and $D_{v}^{-1}$. (c) Predicted scaling regimes for the bacterial
turbulence. The green lines represent the $D_{v}-H$ relation. Each regime is
marked by its scaling exponent. Vertical dashed lines: (from left to right)
measured ranges corresponding to Fig. 3a-c respectively.

## Discussion and summary

In summary, we have observed the universal scaling behaviors of bacterial
turbulence in the 2D and 3D limits, which are independent of bacterial
activity, bacterial length, and system height. The universal scaling laws are
confirmed theoretically by considering the hydrodynamics of active fluids
confined between two no-slip walls. From 2D to 3D, the active turbulence
exhibits a two-step transition: spectral scaling governing the small length
scales changes at $H_{0}$ and a characteristic $k^{-1}$-scaling for 3D
turbulence emerges after $H_{1}$. The two critical heights are also found in
evolution of real-space structures. $H_{0}$ emerges from the interplay of the
confinement height $H$ and the cell length; while $H_{1}$ reflects the
competition between the system height and the vortex size. Lastly, a
characteristic $H^{0.5}$-dependence of vortex size is uncovered experimentally
for 3D bacterial turbulence and is derived using stability analysis in
confined active polar systems.

Our hydrodynamic model considers active fluids as a continuum force field and
use image systems to capture the hydrodynamic screening effect of the no-slip
boundaries as the leading-order approximation. Here we do not consider the
non-uniformity in the density distribution, which could be important in
presence of giant density fluctuations [40, 41] and motility-induced phase
separation (MIPS) [42]. Also, the inner structure of the vortex is largely
simplified in the sense that the continuum assumption ignores irregular shapes
and complex kinetics of realistic vortices at small scales (high-$k$).
Nevertheless, the model still captures our experimental findings excellently,
showing that the details above have only minor effects in determining the
scaling properties of active turbulence. One possible reason is that bacterial
turbulence is dominated by hydrodynamic interaction [2], which suppresses the
density fluctuation and MIPS [42]. In addition, our experiments confirm that
the exponential vortex size distribution, which applies for active nematic
turbulence [21, 18], is also a good approximation in bacterial turbulence.
Hence, we believe our theoretical model captures the essential elements of
scaling behaviors of active turbulence from 2D to 3D. The model can be applied
in hydrodynamic-dominated active systems where the vortex-vortex interaction
is negligible and MIPS is absent.

The scaling laws found in our experiments and theories show wide
applicability. The energy spectra and vortex structures reported by other
experiments that employ different bacteria species and/or are conducted at
higher concentrations [15, 31] are consistent with both our experimental
results and theoretical predictions. Even in dense suspensions of mammalian
sperms where the active agents are an order of magnitude larger than bacteria
and are of a distinct propulsion mechanism, the energy spectrum still displays
$k^{-4}$-scaling when $k\gg 1/H,1/D_{v}$, in line with our results [4].
Lastly, it is noteworthy that the existence of the four scaling regimes does
not rely on a specific $D_{v}-H$ relation but only requires that $D_{v}$
increases with $H$ slower than linearly in the 3D limit (so that $1/H$ and
$1/D_{v}(H)$ can intersect). As the $H^{0.5}$ scaling is also reported in
active nematic turbulence [43, 24, 26, 27], we expect these systems to exhibit
similar four scaling regimes from 2D to 3D.

Besides hydrodynamic-interaction-dominated bacterial turbulence confined
between two no-slip solid surfaces, the collective motion of active matter
systems can be dominated by other interactions, like short-ranged steric
interaction [44] and long-ranged electromagnetic interactions [45, 46]. The
confining boundary could be liquid-liquid interfaces [18] or liquid-air
interfaces [47], which provide distinct couplings within active fluids. How
other interactions and boundary types affect the universal scaling is a
question worth further exploring.

## Acknowledgments

Y.P. acknowledges the support from the National Natural Science Foundation of
China (NSFC) (No. 12074406 and No. T2221001). F.M. acknowledges supports from
the National Natural Science Foundation of China (NSFC) (No. 12047503 and No.
12275332), Chinese Academy of Sciences (No. XDPB15), Max Planck Society (Max
Planck Partner Group) and Wenzhou Institute (WIUCASICTP2022).

## Author Contributions

Y.P. and F.M. conceived the project. D.W. performed the experiments and
analyzed the data. F.M., R.G., Y.Y. and X.W. conducted the theoretical
modeling and calculations. D.W., Y.P. and F.M. wrote the paper with input from
M.L. and R.G.. All authors discussed the results and reviewed the manuscript.

##  Appendix A Methods

###  A.1 Fluidic chamber

We make fluidic chambers with glass slides and coverslips as shown in Fig. 6a.
After being filled with bacteria suspension, chambers are sealed with UV-cured
glue (NOA63, Norland Optics Inc.). To tune the chamber heights, we use
polystyrene beads as spacers for chamber height $H<10$ µm, and different tapes
for chambers of $H\gtrsim 7$ µm. Besides chambers of uniform height, sloped
chambers are made using tapes of different heights ($H^{\prime}\neq
H^{\prime\prime}$) and are calibrated per piece. They allow fine tune over $H$
and high-throughput experimentation. In such chambers, the bottom and top
plates are approximately parallel (incline angle $<0.4\degree$). Turbulent
statistics measured in those two types of chambers are same, see Fig. 6b.

Figure 6: Experimental samples.(a) Schematics of a chamber. (b) Velocity
correlation length measured in the sloped chambers ($H^{\prime}\neq
H^{\prime\prime}$) and uniform-height chambers($H^{\prime}=H^{\prime\prime}$).
Symbols and error bars represent mean and one standard deviation respectively.

###  A.2 Bacterial suspension

We employ a wild-type E.coli strain (BW25113) in our experiment. To extend the
swimming period in closed chambers, we insert a light-driven transmembrane
proton pump, proteorhodopsin (PR), into the strain [48]. A typical cell body
is a 3 µm-long rod with 0.8 µm in diameter. We first culture bacteria in
Terrific Broth (TB) culture medium overnight (12-16 h) at 37 $\degree$C. The
saturated suspension is diluted for 100 times with fresh TB medium and then is
further cultured at 30 $\degree$C for 6-8 h. We harvest bacteria and
centrifuge them at $800-1000$ g for 5 min. The bacteria are resuspended in
adjusted Berg’s motility buffer, containing 6.2 mM K2HPO4, 3.8 mM KH2PO4, 67
mM NaCl, 0.1 mM EDTA, 4.5$\times 10^{-3}$% (v/v) TWEEN20, and 0.5% (w/v)
glucose. The final body volume of bacteria is $\phi\approx 0.06$ (number
density $n\approx 3.2\times 10^{10}$ cells/ml).

###  A.3 Cephalexin treatment

Bacterial length is tuned by inducing cephalexin (Solarbio Inc.) of different
concentrations ($\rm{[Cephalexin]}$) at 3.5 h after the second stage of
culturing begins. We measure the length for N$\approx$1000 cell at each
$\rm{[Cephalexin]}$. The length polydispersity is unvaried, see Fig. 7.
Additionally, bacterial motility is not affected by the cephalexin treatment.
The mean swimming speed of single cells are measured in H=8 µm samples
(N$\approx$20 for each concentration) and they remain constant in the tested
antibiotic concentrations, in line Ref. [49, 36]. In the series of experiments
involving cephalexin, we maintain the cell body volume fraction $\phi\approx
0.06$.

Figure 7: Bacterial length distributions after cephalexin treatment.
Distribution of scaled cell length $L_{b}/\langle L_{b}\rangle$. Distributions
are normalized by their respective maxima.

###  A.4 Flow field extraction

Bright field microscopy are performed with an inverted microscope (Nikon
Ti2e). A $40\times$ objective is employed primarily for measurement and a
$10\times$ objective is used to obtain a large field of view, see Fig. 12.
Videos are recorded with an sCMOS camera at 30 fps (PCO.edge). A typical
measurement lasts for 5 s. For each sample of uniform height, $\sim$20
measurements are taken; while in the sloped chamber, $\sim$8 measurements are
performed at each height (location). Recorded videos are analyzed with
particle image velocimetry (MATLAB PIVLab2.37 [50]). We use windows of
$16\times 16$ pixels with 50% overlap for convolution.

###  A.5 Measuring correlation length for $H=$400 µm

For the thickest samples measured ($H=400~{}$µm, empty red marker at the right
extreme of Fig. 2a), flow patterns at $d_{c}=H/2$ are too obscured to extract.
The presented values are asymptotic estimations by extrapolating the trend
measured with $\rho=20n_{0}$ from $d=0-120$ µm, and $\rho=40n_{0}$ from
$d=0-70$ µm.

Figure 8: Variation of the correlation length in $z$. Measurements of $L_{vv}$
as a function of $d$ (N=3-5 samples) are plotted as dots. Symbols with error
bars represent the mean and standard deviation of the binned data.

###  A.6 Data pooling and averaging

For each confinement height $H$, the highest possible mean turbulence speed
$\langle v\rangle_{\rm max}$ is employed to benchmark the activity level of a
particular bacterial suspension, see Fig. 3a-c. $\langle v\rangle_{\rm max}$
is obtained from measurements over multiple samples, which are prepared
following the same protocol. The energy spectrum of a particular $H$ displayed
in Fig. 3d-f and Fig. 4 is averaged over $\langle v\rangle/\langle
v\rangle_{\rm max}>0.5$. Although $L_{vv}$ and $D_{v}$ are found to depend
negligibly on $\langle v\rangle$ after the onset of turbulence (Fig. 9a), same
as reported previously [29, 19, 20]. However, for consistency, the same
criterion on activity is applied for $L_{vv}$ and $D_{v}$. In measuring
$L_{vv}$ as a function of $d$, we set no criteria for bacterial activity as
$\langle v\rangle$ varies with $d$. The raw measurements of $L_{vv}(d)$ are
displayed in Fig. 8.

Figure 9: Independence of vortex structure on activity and bacterial density
(a) Measurements of $L_{vv}$ at different bacterial activity $\langle
v\rangle$. Different symbols represent different experiments. Measurements are
perform at $d_{c}=H/2$. (b) $L_{vv}$ as a function of bacterial volume
fraction $\phi$ ($H$=85µm, taken at $d_{c}=H/2$). (c) $D_{v}$ as a function of
$\phi$. Black arrows mark the bacterial density used for the presented data in
the main text.

##  Appendix B Independence of turbulence structure on bacterial density

Bacteria are known to accumulate near no-slip boundaries [35]. Such
accumulation may induce number density variation along the system’s $z$ axis.
If the characteristic lengths of turbulence (e.g., $L_{vv}$ and $D_{v}$)
depend on bacterial number density, then the near-wall accumulation will
induce uncertainties in the reported data, for example, in $L_{vv}$’s
dependence on $d$ (Fig. 2). However, we show that bacterial number density
matters negligibly when it is sufficient for turbulence to emerge. Both
$L_{vv}$ and $D_{v}$ remain almost constant after the bacterial volume ratio
$\phi$ exceeds $\sim$3%, see Fig. 9b and c respectively. The presented data
are measured at $d_{c}=H/2$, and we confirm that the trends are qualitatively
the same for other heights ($d$).

##  Appendix C Critical heights measured in different representations.

The critical heights $H_{0}$ and $H_{1}$ manifest both in the characteristic
lengths of the turbulence against $H$ ($L_{vv}$ and $D_{v}$) and in the
kinetic energy spectra. Here we show that the measured values match to each
other. $H_{0}$ and $H_{1}$ measured in different representations are displayed
in Fig. 10a and b respectively. The $x$-axis of Fig. 10a shows $H_{0}$
measured as where $\gamma_{2}$ drops below -2 in the energy spectra; while the
$y$-axis is $H_{0}$ measured as the $x$-intercept of the linear fitting to the
initial fast-increasing trend in $L_{vv}$. The $x$-axis of Fig. 10b reports
$H_{1}$ measured in energy spectra as where $\gamma_{1}$ drops below 1. The
$y$-axis reports $H_{1}$ as where $H=D_{v}$ (with linear interpolation). Data
gather around the identity lines ($y=x$, gray lines), showing that the two
approaches are equivalent.

Figure 10: Critical heights extracted from velocity correlation functions and
from energy spectra. $H_{0}$ (a) and $H_{1}$ (b) measured in different
representations. The solid lines are the identity line ($y=x$). The error bars
in the $x$ and $y$ axes represent the uncertainties in $H$. Figure 11:
Enhanced particle tracking velocimetry to measure 3D flows. Scheme of
extracting the $z$-position of beads in particle tracking velocimetry (PTV).
(a) The normalized intensity ($\tilde{I}$) and (b) tracking quality factor
($\tilde{Q}$) of a fixed bead. Upper images display the beads at different
relative heights $Z_{f}$. Insets: relative heights between the focal plane
(solid line) and the bead, corresponding to the upper images respectively. (c)
The composite metric $\tilde{Q}-2\tilde{I}$ corresponds linearly to the
$z$-position. Dots: single beads; error bar: mean$\pm$std; dashed line: linear
fit. (d) The mean ratio of the $z$-component of speed, $u_{z}$, to the
$xy$-component of speed, $u_{xy}$. Insets: representative histograms of
$u_{z}$ and $u_{xy}$ in samples of $H=30$ µm (left) and 80 µm (right).
Histograms are vertically shifted for clarity.

##  Appendix D Measuring vertical flow using enhanced particle tracking
velocimetry.

To confirm the velocity field has become isotropic at $H=H_{1}$, we measure
the $z$-component of speed $u_{z}$ in the following way. We first develop an
enhanced particle tracking velocimetry (PTV) that extracts a bead’s $z$-axis
position by its appearance. The scheme of the technique is displayed in Fig.
11a-c. Beads are tracked with the TrackMate plugin of ImageJ and the algorithm
primarily depends on fitting the beads’ interior with a Laplacian-of-Gaussian
(LoG) kernel. We exploit the mean pixel intensity $I$ of the beads’ interior
and the quality factor $Q$ which measures how well the LoG kernel fits. Fig.
11a-b display how the normalized metrics ($\tilde{I}$ and $\tilde{Q}$) varies
with the focal plane height ($Z_{f}$, defined as shown in the schematic
drawing). $\tilde{I}$ saturates at $\sim$1 when $Z_{f}\leq 0$ and drops for
$Z_{f}>0$, i.e. the bead center turning from white to black. On the other
hand, $\tilde{Q}$ peaks at $Z_{f}=0$, i.e. the bead that is in-focus has the
highest quality factor. These two trends are illustrated by the images of
beads at the top of Fig. 11a. Combine these two metrics as
$\tilde{Q}-2\tilde{I}$, we obtain a metric that translates linearly to
$Z_{f}$, see Fig. 11c.

Practically, we employ polystyrene beads of $a=6$ µm diameter for PTV. The
relationship between $\tilde{Q}-2\tilde{I}$ and $Z_{f}$ is calibrated by $\sim
50$ beads per sample after the bacterial suspension has ceased moving. Such
configuration allows us to track all three coordinates of the beads in a slice
of $\sim 2a=12$ µm thick. When a beads goes out of this slice, a track is
considered finished. We pool $2\times 10^{3}-1\times 10^{4}$ tracks per
sample. All tracks are broken down as 1 s long segments. This duration for
each segment is determined by the characteristic time scale of the temporal
velocity auto-correlation (1-3 s). The statistics of the speed (total arc
length/1 s) in the $z$-axis and that in the $xy$-plane are displayed in the
insets of Fig. 11d. Altogether, the mean or median of $u_{z}/u_{xy}$ converges
to $\sqrt{2}$ when $H=40-50$ µm, which matches with $H_{1}$ observed from the
vortex size (Fig. 2c) or kinetic energy spectra (Fig. 3h). The ratio supports
that bacterial turbulence turns isotropic ($\langle u_{x}\rangle\approx\langle
u_{y}\rangle\approx\langle u_{z}\rangle$).

##  Appendix E Energy spectrum at low-$k$

Our theory predicts that the scaling $E(k)\sim k$ dominates the low-$k$ end of
the spectrum. To observe examine such spectral behavior for $H>50$ µm, we
measure samples with 10$\times$ magnification. Under this low magnification,
we are able to extract the flow field in samples up to $H$=190 µm, see Fig.
12. The $E(k)\sim k$ scaling regime presents clearly while the spectra
continue shifting to the left (Fig. 3d-e).

Figure 12: Kinetic energy spectrum at the low-$k$ end. Flow fields acquired
with lower microscope magnification (10$\times$) are used to probe the
spectral scaling at the low wavenumbers.

##  Appendix F Exponential distribution of vortex size

To validate one of theoretical assumptions that vortices in bacterial
turbulence follow exponential size distributions, we measure vortex sizes with
the algorithm used in Ref.[21]. The algorithm first computes the Okubo-Weiss
field
$Q=-\partial_{x}u_{x}\cdot\partial_{y}u_{y}+\partial_{x}u_{y}\cdot\partial_{y}u_{x}$
from numerical data of the flow field. A region with $Q<0$ and containing a
flow circulation (vortex core) is considered as a vortex. Areas of such
regions measure the vortex size projected in the $xy$-plane. We implement this
algorithm in samples of varying $H$ and confirm that the vortex size follows
exponential distributions, see Fig. 13.

Figure 13: Exponential distribution of vortex size. Vortex size distribution
in samples of different heights. Cross-sectional sizes of the vortices in the
$xy$-plane are extracted with the algorithm described in Ref. [21]. Dashed
lines are the exponential fits.

## References

  * [1] Dombrowski, C., Cisneros, L., Chatkaew, S., Goldstein, R. E. & Kessler, J. O.  Self-concentration and large-scale coherence in bacterial dynamics.  _Phys. Rev. Lett._ 93, 098103 (2004).
  * [2] Peng, Y., Liu, Z. & Cheng, X.  Imaging the emergence of bacterial turbulence: Phase diagram and transition kinetics.  _Sci. Adv._ 7, eabd1240 (2021).
  * [3] Liu, S., Shankar, M. C., S.and Marchetti & Wu, Y.  Viscoelastic control of spatiotemporal order in bacterial active matter.  _Nature_ 590, 80–84 (2021).
  * [4] Creppy, A., Praud, O., Druart, X., Kohnke, P. L. & Plouraboué, F.  Turbulence of swarming sperm.  _Phys. Rev. E_ 92, 032722 (2015).
  * [5] Martínez-Prat, B., Ignés-Mullol, J., Casademunt, J. & Sagués, F.  Selection mechanism at the onset of active turbulence.  _Nat. Phys._ 15, 362–366 (2019).
  * [6] Wu, K.-T. _et al._ Transition from turbulent to coherent flows in confined three-dimensional active fluids.  _Science_ 355, eaal1979 (2017).
  * [7] Duclos, G. _et al._ Topological structure and dynamics of three-dimensional active nematics.  _Science_ 367, 1120–1124 (2020).
  * [8] Blanch-Mercader, C. _et al._ Turbulent dynamics of epithelial cell cultures.  _Phys. Rev. Lett._ 120, 208101 (2018).
  * [9] Lin, S.-Z., Zhang, W.-Y., Bi, D., Li, B. & Feng, X.-Q.  Energetics of mesoscale cell turbulence in two-dimensional monolayers.  _Commun. Phys._ 4, 21 (2021).
  * [10] Karani, H., Pradillo, G. E. & Vlahovska, P. M.  Tuning the random walk of active colloids: From individual run-and-tumble to dynamic clustering.  _Phys. Rev. Lett._ 123, 208002 (2019).
  * [11] Alert, R., Casademunt, J. & Joanny, J.-F.  Active turbulence.  _Annu. Rev. Condens. Matter Phys._ 13, 143–170 (2022).
  * [12] Kolmogorov, A.  The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds’ Numbers.  _Akademiia Nauk SSSR Doklady_ 30, 301–305 (1941).
  * [13] Boffetta, G. & Ecke, R. E.  Two-dimensional turbulence.  _Annu. Rev. Fluid Mech._ 44, 427–451 (2012).
  * [14] Kellay, H. & Goldburg, W. I.  Two-dimensional turbulence: a review of some recent experiments.  _Rep. Prog. Phys._ 65, 845–894 (2002).
  * [15] Wensink, H. H. _et al._ Meso-scale turbulence in living fluids.  _Proc. Natl. Acad. Sci. U.S.A._ 109, 14308–14313 (2012).
  * [16] Bratanov, V., Jenko, F. & Frey, E.  New class of turbulence in active fluids.  _Proc. Natl. Acad. Sci. U.S.A._ 112, 15048–15053 (2015).
  * [17] Alert, R., Joanny, J.-F. & Casademunt, J.  Universal scaling of active nematic turbulence.  _Nat. Phys._ 16, 682–688 (2020).
  * [18] Martínez-Prat, B. _et al._ Scaling regimes of active turbulence with external dissipation.  _Phys. Rev. X_ 11, 031065 (2021).
  * [19] Dunkel, J. _et al._ Fluid dynamics of bacterial turbulence.  _Phys. Rev. Lett._ 110, 228102 (2013).
  * [20] Heidenreich, S., Dunkel, J., Klapp, S. H. L. & Bär, M.  Hydrodynamic length-scale selection in microswimmer suspensions.  _Phys. Rev. E_ 94, 020601 (2016).
  * [21] Giomi, L.  Geometry and topology of turbulence in active nematics.  _Phys. Rev. X_ 5, 031003 (2015).
  * [22] Wioland, H., Woodhouse, F. G., Dunkel, J., Kessler, J. O. & Goldstein, R. E.  Confinement stabilizes a bacterial suspension into a spiral vortex.  _Phys. Rev. Lett._ 110, 268102 (2013).
  * [23] Lushi, E., Wioland, H. & Goldstein, R. E.  Fluid flows created by swimming bacteria drive self-organization in confined suspensions.  _Proc. Natl. Acad. Sci. U.S.A._ 111, 9733–9738 (2014).
  * [24] Wioland, H., Lushi, E. & Goldstein, R. E.  Directed collective motion of bacteria under channel confinement.  _New J. Phys._ 18, 075002 (2016).
  * [25] Martinez, V. A. _et al._ A combined rheometry and imaging study of viscosity reduction in bacterial suspensions.  _Proc. Natl. Acad. Sci. U.S.A._ 117, 2326–2331 (2020).
  * [26] Chandrakar, P. _et al._ Confinement controls the bend instability of three-dimensional active liquid crystals.  _Phys. Rev. Lett._ 125, 257801 (2020).
  * [27] Varghese, M., Baskaran, A., Hagan, M. F. & Baskaran, A.  Confinement-induced self-pumping in 3d active fluids.  _Phys. Rev. Lett._ 125, 268003 (2020).
  * [28] Nishiguchi, D. & Sano, M.  Mesoscopic turbulence and local order in janus particles self-propelling under an ac electric field.  _Phys. Rev. E_ 92, 052309 (2015).
  * [29] Sokolov, A. & Aranson, I. S.  Physical properties of collective motion in suspensions of bacteria.  _Phys. Rev. Lett._ 109, 248109 (2012).
  * [30] Guo, S., Samanta, D., Peng, Y., Xu, X. & Cheng, X.  Symmetric shear banding and swarming vortices in bacterial superfluids.  _Proc. Natl. Acad. Sci. U.S.A._ 115, 7212–7217 (2018).
  * [31] Liu, Z., Zeng, W., Ma, X. & Cheng, X.  Density fluctuations and energy spectra of 3d bacterial suspensions.  _Soft Matter_ 17, 10806–10817 (2021).
  * [32] Saintillan, D. & Shelley, M. J.  Instabilities and pattern formation in active particle suspensions: Kinetic theory and continuum simulations.  _Phys. Rev. Lett._ 100, 178103 (2008).
  * [33] Aditi Simha, R. & Ramaswamy, S.  Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles.  _Phys. Rev. Lett._ 89, 058101 (2002).
  * [34] Pope, S. B.  _Turbulent flows_ (Cambridge, 2000).
  * [35] Berke, A. P., Turner, L., Berg, H. C. & Lauga, E.  Hydrodynamic attraction of swimming microorganisms by surfaces.  _Phys. Rev. Lett._ 101, 038102 (2008).
  * [36] Li, H. _et al._ Data-driven quantitative modeling of bacterial active nematics.  _Proc. Natl. Acad. Sci. U.S.A._ 116, 777–785 (2019).
  * [37] Liron, N. & Mochon, S.  Stokes flow for a stokeslet between two parallel flat plates.  _J. Eng. Math._ 10, 287–303 (1976).
  * [38] Meng, F., Bennett, R. R., Uchida, N. & Golestanian, R.  Conditions for metachronal coordination in arrays of model cilia.  _Proc. Natl. Acad. Sci. U.S.A._ 118, e2102828118 (2021).
  * [39] Ishida, S., Yang, Y., Meng, F. & Matsunaga, D.  Field-controlling patterns of sheared ferrofluid droplets.  _Phys. Fluids_ 34, 063309 (2022).
  * [40] Narayan, V., Ramaswamy, S. & Menon, N.  Long-lived giant number fluctuations in a swarming granular nematic.  _Science_ 317, 105–108 (2007).
  * [41] Zhang, H. P., Be’er, A., Florin, E.-L. & Swinney, H. L.  Collective motion and density fluctuations in bacterial colonies.  _Proc. Natl. Acad. Sci. U.S.A._ 107, 13626–13630 (2010).
  * [42] Cates, M. E. & Tailleur, J.  Motility-induced phase separation.  _Annu. Rev. Condens. Matter Phys._ 6, 219–244 (2015).
  * [43] Brotto, T., Caussin, J.-B., Lauga, E. & Bartolo, D.  Hydrodynamics of confined active fluids.  _Phys. Rev. Lett._ 110, 038101 (2013).
  * [44] Patteson, A. E., Gopinath, A. & Arratia, P. E.  The propagation of active-passive interfaces in bacterial swarms.  _Nat. Commun._ 9, 5373 (2018).
  * [45] Meng, F., Matsunaga, D. & Golestanian, R.  Clustering of magnetic swimmers in a poiseuille flow.  _Phys. Rev. Lett._ 120, 188101 (2018).
  * [46] Meng, F., Matsunaga, D., Mahault, B. & Golestanian, R.  Magnetic microswimmers exhibit bose-einstein-like condensation.  _Phys. Rev. Lett._ 126, 078001 (2021).
  * [47] Peng, Y. _et al._ Diffusion of ellipsoids in bacterial suspensions.  _Phys. Rev. Lett._ 116, 068303 (2016).
  * [48] Walter, J. M., Greenfield, D., Bustamante, C. & Liphardt, J.  Light-powering Escherichia coli with proteorhodopsin.  _Proc. Natl. Acad. Sci. U.S.A._ 104, 2408–2412 (2007).
  * [49] Nishiguchi, D., Nagai, K. H., Chaté, H. & Sano, M.  Long-range nematic order and anomalous fluctuations in suspensions of swimming filamentous bacteria.  _Phys. Rev. E_ 95, 020601 (2017).
  * [50] Thielicke, W. & Stamhuis, E. J.  Pivlab – towards user-friendly, affordable and accurate digital particle image velocimetry in matlab.  _J. Open Res. Software_ 2, e30 (2014).