Bend instabilities and topological turbulence in shear-aligned living liquid crystal

Abstract

Flagellated microswimmers B. subtilis dispersed in a nematic phase of a lyotropic chromonic liquid crystal form a living liquid crystal (LLC). The combination of the passive and active components allows us to analyze how the active component transitions from the shear-imposed alignment into topological turbulence. The lateral extension of the experimental cell is 1000 times larger than the 10-micrometer thickness, to avoid the effect of lateral confinement on the dynamic patterns. The surface anchoring is azimuthally degenerate to avoid permanent anisotropy. After the shear cessation, active forces produced by the microswimmers trigger self-amplifying bend undulations. The amplitude A of undulations grows with time and then saturates, while the wavelength λ increases only slightly. Strong bend stresses at the extrema of undulations are released by nucleating ±1/2 disclination pairs, which multiply to produce topological turbulence. The nucleating pairs and the symmetry axes of the +1/2 disclinations are orientationally ordered. In highly active LLCs, this alignment causes a nonmonotonous time dependence of the disclinations' number, as a fast-moving +1/2 disclination of one pair annihilates a −1/2 disclination of the neighboring pair. The spectra of elastic and kinetic energies exhibit distinct wavevector dependencies caused by the energy cost of the deformed passive nematic background. The study demonstrates how applied shears and a passive viscoelastic background affect the dynamics of active matter. Combined with the previously determined viscoelastic properties of the lyotropic chromonic liquid crystal, the results complete a comprehensive description of the experimentally assessable type of active matter.

---

1. Introduction

Active nematics are the most studied examples of active matter.^1–7 Their experimental realizations are microtubules activated by molecular motors,^8–11 tissue monolayers,^12–18 bacterial colonies,^19–21 and aqueous suspensions of swimming bacteria.^22–27 A uniform active nematic formed by extensile active units, aligned by shear in microfluidic channels¹¹ or by pressure-driven flows,²⁸ is unstable against the formation of bend undulations in the orientation of neighboring units and chaotic flows called “topological turbulence”. These instabilities, predicted theoretically^1,29–31 are readily observed in aqueous dispersions of bacteria^32–38 and microtubules activated by molecular motors.^8,11,28

A different class of active matter has been recently introduced,³⁹ explored,^40–53 and reviewed^27,54,55 as a living liquid crystal (LLC), also known as a living nematic. In an LLC, an isotropic water solvent is substituted with a (passive) nematic (N), a water-based nontoxic lyotropic chromonic liquid crystal,^56–59 thanks to which the solvent acquires viscoelastic properties. In nature, many bacterial types move in viscoelastic fluids,^60,61 some of which, such as mucus^62,63 and extracellular matrix,^64,65 exhibit orientational order⁶⁶ or impart such an order onto the dispersed microswimmers.^67–69

The orientational order and anchoring of the N at the bacterium surface align bacteria with an elongated shape along the local director n̂ of N,^39,70–76 which is apolar, n̂ = −n̂. Although the individual bacteria are polar, their global dispersion in a uniform LLC is not polar: the number of bacteria swimming up and down the director is statistically the same. Apolar LLC behavior is also evidenced by orientational defects in the LLC, disclinations, which are of a semi-integer strength: bacterial orientations change by ±π when one circumnavigates the defect core, thanks to the identity n̂ = −n̂.³⁹ Polar flows of bacteria in an LLC, either circular,^40,47 or rectilinear,⁴⁶ can be created by imposing gradients of n̂ through patterned surface anchoring of the N at bounding plates.⁷⁷ These surface anchoring patterns of the passive N command the dynamics of bacteria, producing directional flows,^{40,46,47,53,72,78} and targeted delivery of a micro cargo.^46,73,74 However, when the activity exceeds some threshold, the LLC develops dynamic instabilities of orientation.^{39,42,43,45–47,52,53,79} Since LLCs combine orientational order, viscoelasticity, and activity in a tunable manner and since they hold promise for applications such as controlled microscale transport, understanding activity-triggered instabilities in them is of fundamental and applied interest.

So far, the LLC has been explored in pre-aligned cells with guiding rails of surface anchoring imposed onto n̂.^39–52 Such alignment causes permanent anisotropy and makes it difficult to separate the intrinsic properties of LLC from the effect of surface alignment. In this work, to avoid the externally imposed permanent anisotropy, we study the evolution dynamics of an initially shear-aligned LLC to a topologically turbulent LLC in a macroscopically isotropic environment, using flat cells with no pre-imposed in-plane surface anchoring. The initial uniform alignment of the LLC is produced by a temporal unidirectional shear. Once the shear stops, hydrodynamic interactions of bacteria create active flows that produce bend instability and nucleation of disclinations; this behavior is contrasted to the relaxation of the passive N to a uniform state. We measure the undulation wavelength, amplitude, and growth rate of the LLC instabilities as a function of bacterial concentration. The nucleation and evolution of disclinations are quantified by a newly developed Delaunay mesh approach. Experiments enable us to extract density–density, orientation–orientation, and velocity–velocity correlation functions, analyze density fluctuations, and calculate the spectra of elastic energy, kinetic energy, and enstrophy. These multifaceted data present a comprehensive experimental picture of the transition scenario from an oriented LLC to topological turbulence controlled by the balance of active stresses and underlying elasticity in the absence of permanent external anisotropy.

2. Methods

2.1 Living liquid crystal preparation

We use a Bacillus subtilis strain 1085, which is a rod-shaped bacterium with a head that is 3–10 µm long and ∼0.7 µm in diameter. The strain is stored at −80 °C. The bacterial colony is grown on Lysogeny broth (Miller composition from Teknova, Inc.) agar plates at 35 °C for 12–24 h, then transferred to a vial with 10 mL of aqueous terrific broth (TB) (Sigma T5574) and grown in a shaking incubator at 35 °C for 7–9 h. The vials are sealed to increase the resistance of bacteria to oxygen starvation. The bacteria concentration in the growth medium is monitored by measuring optical density. The dispersions are removed from the incubator at the end of the exponential growth stage, at a concentration of c₀ = 8 × 10¹⁴ cells per m³. The bacteria are extracted from the liquid medium by centrifugation.

The N-forming organic component disodium cromoglycate (DSCG) (Spectrum Chemicals, purity >98%) is dissolved in the TB at a concentration 13 wt%. The solution is a homogeneous nematic at the temperature (25.0 ± 0.1), at which all experiments are performed. The extracted bacteria are mixed with the DSCG solution to form the LLC.

The chromonic N is not toxic to bacteria.⁵⁹ A building unit of the chromonic N is a cylindrical aggregate of ∼2 nm in diameter and a length that ranges from 1 to 10³ nm.^56–58 The gap between the neighboring aggregates is typically (2–3) nm.⁸⁰ For a bacterium such as B. subtilis, with a rodlike body of a length (3–10) µm and a diameter of - 0.7 µm, with attached flagella of a length (5–15) µm, such an environment appears devoid of structure. However, anisotropic surface interactions at the bacterium surface orient the N aggregates along the bacterial axis, so that the swimming direction is collinear with the director n̂ = −n̂ of the passive N.^39,70–76 In the LLC, the bacteria interact through surface anchoring-mediated elastic forces of the N background and by hydrodynamics, as each bacterium produces two fluid jets directed outwards along its axis.⁴² The orientational elasticity and the activity of the LLC could be adjusted independently, by the concentration of chromonic molecules that controls the viscosity and elasticity of the N^81–83 and by the concentration of bacteria, which controls the levels of activity.³⁹

2.2 Shearing cell

A Linkam optical shearing system CSS450 with plate–plate geometry is used to shear-align the LLCs. The top and bottom quartz plates are parallel and allow one to observe the sample in a transmission mode. The windows of 2.5 mm in diameter are centered at 6.25 mm from the axis of rotation, Fig. 1. The shearing direction is along the x-axis in Fig. 1a. The gap distance h between the two plates is 10 µm, which is much smaller than the diameter of the sample, about 7.5 mm. One plate rotates with a controllable angular velocity in the range (0.001–10) s⁻¹, which produces shear rates in the range of 0.75 ≤ [small gamma, Greek, dot above] ≤ 7500 s⁻¹. Each plate is in thermal contact with a silver heater, which sets the temperature at (25.0 ± 0.1) °C. In contrast to the previous studies,^39,40,46,47 no alignment layers are used; the plates are cleaned with deionized water and acetone before each use and impose degenerate tangential alignment of the N director n̂. In the absence of shear, n̂ assumes any orientation in the plane of cell. To demonstrate this, we subjected the passive N samples of DSCG to shear at a rate of 10 s⁻¹ and then stopped the shear and allowed the material to relax for 15 min. The relaxation resulted in different orientations of the director in the plane of the sample; see ref. 84 for details.

Fig. 1 (a) Scheme of shear; (b) velocities of bacteria; (c) orientational parameters of bacteria.

2.3 Imaging

Olympus BX40 and BX50 polarizing optical microscopes equipped with a camera AM Scope MU 2000 (5 fps, field of view 730 µm by 486 µm) are used to record videos of the LLCs during and after the shear. We also use polychromatic polarization microscopy (PPM) developed by Shribak,^85,86 equipped with a spectral polarization state generator and an analyzer. An observation window of CSS450 allows one to determine the parameters of interest, such as velocities and orientations of the bacteria, Fig. 1(b) and (c).

2.4 Analysis

All collected videos are processed using a homemade Matlab algorithm to enhance the contrast. Matlab and Mathematica-based codes have been developed to quantify (1) bacterial orientations, (2) the wavelength and amplitude of undulations, (3) positions of bacteria, (4) correlation functions, (5) energy spectra, and (6) density modulations. All codes are described in the SI.

To assure statistical reproducibility and robustness of the data, we used standardized protocols and settings (the shear rate, shear chamber geometry, absence of azimuthal anchoring, passive nematic composition, etc.); the only parameter changing was the concentration of microswimmers that varied from 0.5c₀ to 10c₀. The LLC behavior was also contrasted to the passive counterpart, with no bacteria. The LLC dynamics was recorded in multiple videos for each concentration. We identified 2000–3000 bacteria in each of the 370–650 frames of each video to extract the system parameters.

3. Results and discussion

An LLC is a dispersion of motile B. subtilis in an N formed by 13 wt% solution of DSCG. B. subtilis swims parallel to the local N director by coordinated rotation of about 10–20 helical filaments.^39,74 A swimming B. subtilis generates a hydrodynamic force dipole:⁸⁷ the surrounding fluid is pushed away from the bacterium along the axial direction and pulled towards the bacterium along two perpendicular directions. The development of LLC patterns upon cessation of shear depends on the concentration c of bacteria, varied from c = 0.5c₀ to 10c₀; c controls the activity level. We also explore an inactive N with c = 0, as described below.

3.1 Shear-induced patterns

3.1.1 Shear-induced patterns of passive N. Under shear, the passive N experiences a cascade of structural transformations, indicative of tumbling behavior.⁸⁴ At low shear rates, [small gamma, Greek, dot above] ≤ 1 s⁻¹, the director n̂ realigns perpendicularly to the shear plane formed by the shear direction (the x-axis) and the shear gradient (the z-axis), entering the so-called log-rolling regime.⁸⁴ As [small gamma, Greek, dot above] increases to (1–100) s⁻¹, disclination loops nucleate.^84,88 At [small gamma, Greek, dot above] > 500 s⁻¹, the texture develops stripes,⁸⁴ in which n̂ is predominantly along the shear direction, with small periodic clockwise and counterclockwise ≈15° tilts away from the shear plane.

3.1.2 Shear-induced patterns of LLC. In the presence of bacteria, c = (0.5–10)c₀, at a low shear rate [small gamma, Greek, dot above] ≤ 1 s⁻¹, the LLC director experiences a log-rolling regime, aligning along the y-axis, SI Fig. S1(a). Shear rates [small gamma, Greek, dot above] ∼ 1–100 s⁻¹ produce polydomain textures, SI Fig. S1(b). At [small gamma, Greek, dot above] ∼ 1000 s⁻¹, the LLC develops stripes, similar to those in the passive N,⁸⁴ in which the bacteria align along the flow direction, a behavior different from the one reported for isotropic B. subtilis dispersions in water under the shear.⁸⁹ In the latter case, shear combined with the chirality of the left-handed flagellum realigns the bacteria counterclockwise, so that the microswimmers acquire a “rheotactic” velocity along the y-axis.⁸⁹ We did not notice this effect in the LLC, presumably because the shear-aligned director prevents the realignment of bacteria.

3.2 Relaxation upon shear cessation

3.2.1 Relaxation of passive N upon shear cessation. Once the shear at [small gamma, Greek, dot above] = 1000 s⁻¹ stops at t = 0 s, the passive N develops stripes parallel to the shear x-direction, which at t ≈ 1 s are replaced by bands along the vorticity y-axis, SI Fig. S2(a) and Video S1. Similar bands have been observed after cessation of flow in liquid crystal polymers,^90–92 nematic suspensions of rodlike fd-viruses,^93,94 and cellulose cholesterics.^95,96 The bands represent bend undulations of the N director; their wavelength λ₀ increases with time, approximately linearly from ∼18 µm at t = 0.2 s to 36 µm at t = 2.4 s. At t > 3 s, the undulations transform into a coarsening polydomain texture, SI Fig. S2(a). The relaxation of active LLC upon shear cessation is very different, as described below.

3.2.2 Relaxation of LLC upon shear cessation. Low concentrations, c = (0.5–2)c₀, exhibit a peculiar scenario with long-time memory of the shear-aligned state, as the bacteria continue to swim up and down the x-axis for a long time, about 3 min after the cessation of flow, SI Fig. S3. The elongated bacteria, aligned by shear, anchor the N director n̂, stabilizing the uniform alignment of LLC. The concentrations c = (0.5–2)c₀ correspond to the bacterial separations ≈(20–40) µm, which are too large for the hydrodynamic force dipoles of neighboring bacteria to overlap and to cause collective behavior and instabilities.⁴⁷

High concentrations, c = 3c₀, 5c₀, and 10c₀ show dynamics illustrated by Videos S2–S5. In these LLC compositions, the separations of bacteria are short, 20 µm or less, permitting hydrodynamic interactions that cause a cascade of instabilities, as illustrated in Fig. 2 for 5c₀. Upon cessation of shear, at t = 1–2 s, the bacteria swim up and down the x-axis. At t > 2 s, the bacterial orientations and the background N director n̂ develop periodic bends, Fig. 2 and SI Fig. S2(b). Their period λ is in the range (50–110) µm, depending on c, and is noticeably larger than the wavelength λ₀ of the short-lived undulations in the inactive N. In what follows, we present details of the activity-induced instabilities in 3c₀, 5c₀, and 10c₀ compositions.

Fig. 2 Evolution of the LLC, c = 5c₀, after shear cessation: (a) color-coded orientation of bacterial bodies as established by observations under an optical microscope in unpolarized light; (b) scenarios of LLC transition from the uniformly aligned state to bend undulations and topological turbulence for bacterial concentrations 3c₀, 5c₀, and 10c₀, see Videos S3, S4 and S5, respectively.

3.3 Bend undulations of LLC upon shear cessation

The orientation of a rodlike bacterium is specified by its axis n̂_b = (cos θ, sin θ), where θ is the angle measured counterclockwise from the shear direction, Fig. 2(a). At the start of relaxation, θ ≈ 0, but the range of θ expands with time, until all reorientations become equally probable in topological turbulence, t = 150 s in Fig. 3. Immediately after the shear cessation, the distribution represents an ellipse elongated along the shear direction. As the bend undulations develop, t = 50 s, the ellipse elongates along the vorticity y-axis; this reshaping is caused by the tilted orientations of bacteria at the shoulders of the undulation wave. Nucleation of disclination at the ridges of the undulations elongates the ellipse even more, t = 100 s, as the bacterial orientation within the nucleating disclination pairs is along the y-axis.

Fig. 3 Histograms of bacterial orientation for 5c₀ dispersion at different times after shear cessation.

Fig. 4 illustrates the development of the wavelength λ(t), and the amplitude A(t) of undulations for c = 3c₀, 5c₀, and 10c₀. These parameters are determined by fitting the experimentally observed trajectories, which transition from sinusoidal to sawtooth profiles,⁹⁷ SI Fig. S4, and Section S6. In 5c₀ and 10c₀ dispersions, λ is practically time-independent; a higher activity brings about a smaller λ, as λ ≈ 75 µm in 5c₀ and ≈50 µm in 10c₀. This behavior is very different from the strong increase in the wavelength λ₀(t) of bend in the passive N, see the inset in Fig. 4(a). The behavior of weakly concentrated 3c₀ LLC is intermediate between strongly active LLC and a passive N, as λ slightly increases with time, from λ ≈ 80 µm at the onset of undulations to ∼100 µm at later stages.

Fig. 4 Development of (a) the wavelength λ(t) and (b) the amplitude A(t) of undulations for 3c₀, 5c₀, and 10c₀ LLCs. The inset in (a) shows the wavelength λ₀(t) of undulations in the passive N after the shear cessation. The inset in (b) shows a fitted A(t) for 10c₀.

The amplitude A, calculated as half the distance between the maximum and minimum of the bending wave, first increases with time and then saturates, Fig. 4(b). The overall growth and saturation could be approximated by the function , where A_s is the amplitude at saturation, τ is the characteristic time at which A starts to saturate, and t₀ is the offset time (t₀ = −4.3 s for 3c₀, 2.4 s for 5c₀, and 0.6 s for 10c₀), associated with a delay in the undulation onset. The saturation time decreases strongly with the increase of activity/concentration: τ ≈ (54.2 ± 2.5) s for 3c₀, (39.2 ± 1.1) s 5c₀, and (3.2 ± 0.3) s for 10c₀. The amplitudes A_s do not change much: A_s ≈ (20.6 ± 0.5) µm in 3c₀, (19.9 ± 0.3) µm in 5c₀, and (16.4 ± 0.3) µm in 10c₀. These values are close to the length scale ≈(15–20) µm at which the hydrodynamic force dipoles of neighboring microswimmers overlap.⁴⁷

3.4 Nucleation of defect pairs in LLC

As the amplitude of undulations increases, while the wavelength remains practically constant, the gradients in the orientational field become stronger, reaching maxima at the extrema of the undulation wave. The increasing stress is relieved by the nucleation of pairs of ±1/2 disclinations in bacterial orientations. Below, we describe how these defects are identified.

3.4.1 Identification of defects. To ensure an accurate determination of bacterial orientations, we impose the following criteria on bacterial images: (1) each image is large, with an area above 30 pixels; and (2) the major axis is at least 8 pixels (3.2 µm) and the major axis/minor axis ratio is larger than 1.75. The typical number of bacterial images that satisfy these requirements in one frame ranges from 2600 for 3c₀ to 3000 for 10c₀.

Since the N director n̂ is tangential to the bacterial body,⁴² the long axis n̂^(j)_b = (cos θ^(j), sin θ^(j)) of a bacterium j making an angle θ^(j) with the shear direction x represents the local director n̂(r^(j)) of the LLC; here r^(j) is the center of mass of a pixelated bacterial image. The identification of disclinations illustrated in Fig. 5(a)–(c) involves the following steps:

Fig. 5 Identification and evolution of disclinations in an LLC: (a) a frame at the onset of topological turbulence contains n₊ red triangles (containing +1/2 disclinations), n₋ green ones (with −1/2 disclinations), and n_BP BPs composed of magenta (+1/2) and cyan (−1/2) triangles; (b) fragment of frame (a) where the orientation of +1/2 disclination is defined by the unit vector p along the radial director toward the defect core and the vector d connects the −1/2 to the nearest +1/2 disclinations; (c) fragment of frame (a) with BPs; (d) fragment of Fig. S5(g) where the pair selection connects the disclinations nucleated in different events. The evolution of disclinations for 3c₀ (red), 5c₀ (black), and 10c₀ (green) dispersions: (e) number densities of the +1/2 disclinations ρ₊ and of the BP pairs ρ_BP; (f) the orientational order of +1/2 disclinations 〈cos 2φ_p〉 and of unit vectors d̂ = d/|d|〈cos 2φ_d〉; (g) average of the mutual orientation of the paired disclinations 〈p·d̂〉; (h) a scheme showing how +1/2 and −1/2 disclinations from two different pairs annihilate as a result of self-propulsion of the +1/2 defect and how the p vector of +1/2 disclination loses its initial orientation along the y-axis through rotational diffusion.

(1) The image is fragmented into a set of triangles with vertices at each bacterium location r^(j), Fig. 5(a). We use the Delaunay triangulation algorithm, provided by Mathematica.

(2) A disclination of a strength m implies that circumnavigation around its core results in a director rotation θ = 2πm; counterclockwise circulation and director rotation are considered positive. Using each triangle as a loop, we calculate the director rotation between vertices k and k′ as θ_k′k = θ_k′ − θ_k + lπ, where θ_k, k = 1, 2, 3, are the angles of three bacteria at the vertices, and the integer l keeps θ_k′k in the range (−π/2, π/2). Then, a non-zero value of the total rotation Δθ = θ₂₁ + θ₃₂ + θ₁₃ around the triangle indicates that it contains a disclination of a strength m = Δθ/2π = ±1/2, as illustrated in Fig. 5(b). In Fig. 5(b), θ₁ = −0.20π, θ₂ = −0.04π and θ₃ = 0.37π, which result in Δθ = 0.16π + 0.41π + 0.43π = π and m = +1/2.

(3) The highlighted triangles in Fig. 5(a)–(c) are of the following types:

(a) Containing an isolated disclination, either +1/2 (red triangles) or −1/2 (green triangles), surrounded by defect-free triangles, Fig. 5(b). The orientation of an isolated +1/2 disclination is defined by the unit vector p along the radial director pointing toward the defect core. Since disclinations nucleate in pairs, we introduce a vector d that connects a −1/2 disclination to the nearest +1/2 disclination, Fig. 5(b); the pairs are selected by minimizing the total length of distances |d| in the frame.

(b) Forming bound pairs (BPs), in which a +1/2 triangle shares a common side with a −1/2 triangle; these are marked with magenta and cyan colors, respectively, in Fig. 5(c).

(c) Forming clusters that contain three or more connected triangles, Fig. 5(a). A cluster of 2n triangles has an equal number of ±1/2 defects and can be considered as a connected set of n BPs, whereas a cluster of 2n + 1 triangles can be considered a connected set of n BPs and a triangle with an isolated defect.

The method above allows one to identify the disclinations and analyze their evolution, as discussed below.

3.4.2 Evolution of defects. Fig. 5(a) and SI Fig. S5 present selected frames of defects' evolution. Each frame exhibits the numbers n_BP, n₊, and n₋ of BPs and isolated ±1/2 disclinations, respectively. A small random imbalance of n₊, and n₋ results from the motion of disclinations across the frame borders.

The evolution of disclinations is characterized by the number densities of the +1/2 disclinations and of the BP pairs , Fig. 5(e), and by the orientational order of isolated ±1/2 disclinations, Fig. 5(f) and (g), defined through the averages of the unit vectors p and d̂ = d/|d|, Fig. 5(b). Here is the frame area. We start the analysis with the 3c₀ and 5c₀ LLCs that exhibit similar evolution.

3.4.2.1 3c₀ and 5c₀ LLC dispersions. At the early stages after shearing cessation, t ≤ 10 s, only a few BPs form; there are no isolated ±1/2 disclinations, Fig. 5(e) and SI Fig. S5(a) and (d). At t > 10 s, ρ_BP increases, and at t > 15 s, separation of +1/2 and −1/2 disclinations of BPs produces isolated disclinations, Fig. 5(a), (b) and Fig. S5(b), (e). The densities ρ₊ and ρ_BP in 3c₀ and 5c₀ dispersions grow with time, Fig. 5(e). In the well-developed topological turbulence, t > 60 s, ρ_BP continues to increase, while ρ₊ saturates in the 3c₀ and exhibits a small decay in 5c₀ system, Fig. 5(e). This decay in the density of isolated disclinations can be attributed to the well-developed active flows that bring +1/2 and −1/2 disclinations together and stimulate their annihilation.

At the onset of topological turbulence, Fig. 5(a) and (b), the isolated +1/2 disclinations and the defect pairs align along the y-axis, as 〈cos 2φ_p〉 ≈ 0.4 and 〈cos 2φ_d〉 ≈ 0.6, Fig. 5(f). The background N director in the nucleating pair is parallel to the line connecting the two cores: 〈p·d̂〉 ≈ 0.6, Fig. 5(g). This director is perpendicular to the shear direction, which implies that the nucleating defects randomize the alignment introduced by the shear. As time progresses, the orientational order weakens: 〈p·d̂〉 and 〈cos 2φ_p〉 decrease to zero, while 〈cos 2φ_d〉 remains finite, ≈0.2, even at the latest times, t = 124 s, Fig. 5(f). The observations suggest that the orientational diffusion of p is much faster than the translational diffusion of +1/2 cores along the x-axis, Fig. 5(h).

3.4.2.2 10c₀ LLC dispersions. The evolution of the highly active 10c₀ LLC differs substantially from the 3c₀ and 5c₀ LLCs by a nonmonotonic time dependence of ρ_BP and ρ₊, Fig. 5(e). At t = 3 s, multiple defect pairs have already nucleated, preserving orientation of +1/2 disclinations along the y-axis, as 〈cos 2φ_p〉 ≈ 0.5 is larger than 〈cos 2φ_d〉 ≈ 0.3, and 〈p·d〉 ≈ 0.3, Fig. 5(f) and (g). Both ρ_BP and ρ₊ peak at t ≈ 5 s, then undergo a dramatic decrease, reaching minima at t ≈ 11 s, then grow at t > 15 s and eventually saturate in the well-developed topological turbulence, Fig. 5(e). The evolution of ρ_BP and ρ₊ correlates with the behavior of 〈p·d̂〉 in Fig. 5(g). The nonmonotonous dynamics of disclinations is illustrated in Fig. 5(h) and explained below.

In the strongly active LLC, the nucleating defect pairs align along the y-axis. The +1/2 part of the defect pair is mobile; it propels along p.⁹⁸ The −1/2 defects are less mobile because of their symmetry. In the 10c₀ dispersion, both the concentration of nucleating pairs and the speed of +1/2 cores, proportional to the activity, are high. The fast-moving +1/2 defect runs into the −1/2 defect of a nearby pair and annihilates with it, Fig. 5(h). This “alignment and activity-triggered annihilation” decreases the defect density within the time interval (6–20) s, Fig. 5(e). Annihilation implies that the separation between defects temporarily increases; fast orientational diffusion of +1/2 disclinations, Fig. 5(h), hinders further annihilations. The annihilation of +1/2 disclinations aligned along the y-axis happens within the t = (20–40) s interval: 〈cos 2φ_p〉 ≈ −0.15 is negative, which means that defects that did not annihilate orient predominantly along the shear x-direction rather than along the y-axis. Misalignment of defects and a lower probability of annihilation, Fig. 5(f), increases their numbers, Fig. 5(e). This alignment and activity-triggered annihilation is absent in 3c₀ and 5c₀ since the nucleating defects are separated by longer distances.

The quantitative analysis of the LLC instabilities presented above establishes a connection between the number densities of disclinations and their orientational order and uncovers a new effect of the alignment and activity-triggered annihilation in highly active dispersions, Fig. 5(e)–(g).

3.5 Velocity analysis

The velocity v^(j)_n of an jth bacterium in a frame n is calculated as v^(j)_n = [(r^(j)_n − r^(j)_n−1) + (r^(j)_n+1 − r^(j)_n)]/2Δt, where r^(j)_n is the bacterium's location in the nth frame. A bacterium swims along n̂^(j)_b = n̂ = (cos θ, sin θ), where θ is the angle between the bacterium and the x-axis. In the absence of other bacteria, v^(j) is collinear with n̂^(j). However, the velocities v^(j) are often observed at a substantial angle to n̂^(j). The perpendicular v^(j)_⊥ component of v^(j) is attributed to the active flows produced by bacterial hydrodynamic interactions; these flows advect a bacterium along a direction different from n̂^(j). We thus represent v^(j) as a sum of a component parallel to n̂^(j), v^(j)_‖ = (n̂^(j)·v^(j))n̂^(j), and perpendicular to it, v^(j)_⊥ = v^(j) − v^(j)_‖, SI Section S9 and Fig. S6.

In a sheared dispersion, θ = 0 and v^(j)_⊥ = 0, Fig. 2(a). Once the shear stops, a fluctuating bend in the orientation of neighboring bacteria is enhanced by the mechanism proposed by Simha and Ramaswamy^29,99 and illustrated in Fig. 2(b). Each bacterium is a pusher that produces two antiparallel fluid jets along n̂^(j). At short separations of bacteria, ≤(15–20) µm, these jets overlap.⁴⁷ For bacteria on the two shoulders of the bend wave, the overlaps produce a flow orthogonal to the x-axis. This flow amplifies the bend amplitude, while leaving the wavelength λ practically intact; the latter is defined by the balance of active forces which favor a shorter λ and the N elasticity which favors a longer λ.³⁹

The individual velocity v^(j) = v^(j)_s + v̄_a is the sum of the swimming velocity v^(j)_s = v^(j)_sn̂^(j) and the active flow v̄_a; v̄_a can be determined by fitting the individual velocities if v̄_a is assumed to be the same for neighboring bacteria:

v^(j)_⊥ = v^(j)_s − (n̂^(j)·v^(j)_s)n̂^(j) = v̄_a − (n̂^(j)·v̄_a)n̂^(j). (1)

The evolution of color-coded averaged y-component of the active flow velocity v̄_ay, obtained by fitting v^(j)_⊥ with eqn (1) within a disc of a radius 20 µm, Fig. 6(a), shows a strong correlation with patterns of orientation in undulations and turbulence, Fig. 3. The active flows are not simple vector summations of the individual swimming velocities, Fig. 6(b). Since bacteria are force-free swimmers, their propulsion involves two forces of the same magnitude but opposite orientation. The propulsion force F^(j)_prop is parallel to v^(j)_‖ while the drag force F^(j)_drag is antiparallel to v^(j)_‖. The active flow v̄_a could result from the addition of F⁽ⁱ⁾_prop and F^(j)_prop, or F⁽ⁱ⁾_drag and F^(j)_prop, or F⁽ⁱ⁾_drag and F^(j)_drag. For example, scenarios 3 and 4 in Fig. 6(b) show two large vectors v⁽ⁱ⁾_‖ and v^(j)_‖ directed to the right while the active flow v̄_a is directed to the left.

Fig. 6 (a) The color-coded averaged horizontal component v̄_ay of the orthogonal velocity overlayed with the parallel velocities of bacteria v^(j)_‖ shown as black arrows at t = 10, 20, 60, and 120 s in 5c₀; (b) a few enlarged scenarios (taken at t = 20 s) demonstrating that the polarity of the active flow v̄_a is not a result of a simple vector summation of individual velocities v^(j)_‖ of bacteria. Scenarios 1 and 2 involve an active flow to the right, while 3 and 4 involve an active flow to the left.

3.6 Analysis of correlations

In this subsection, we analyze the density–density, orientation–orientation, and velocity–velocity correlation functions.

We start with the pair distribution function P₂(r, φ, θ_l, θ_k) of finding two bacteria k and l, oriented along the directions θ_l and θ_k, and separated by a distance r = r(cos φ, sin φ); here, the angles φ, θ_k and θ_l are calculated counterclockwise from the shear direction x, Fig. 1(c). Using the asymptotic behavior P₂(r → ∞, φ, θ_l, θ_k) → ρP₁(θ_l)P₁(θ_k), where ρ is the density of bacteria per unit area and P₁(θ) is the single bacterium's orientational function that obeys the condition , we introduce the dimensionless correlation function C(r, φ, θ_l, θ_k) = ρ⁻¹P₂(r, φ, θ_l, θ_k) − P₁(θ_l)P₁(θ_k) and represent it as a Fourier series

where the Fourier amplitudes are integral characteristics of C(r, φ, θ_l, θ_k):The indices n_l,k define the type of bacterium's characteristics that the correlation function quantifies. A bacterium is polar, with polarity defined by the swimming direction, Fig. 1(c). For n_k = 0, the kth bacterium is a point-like object; the corresponding correlation function quantifies the density correlations. The terms with n_k = ±1 describe polar correlations caused by the bacterium's polarity. The terms with n_k = ±2 describe the quadrupolar correlations induced by the elongated bacterial body; the body could be represented as an ellipse. The index n defines the azimuthal dependence of the corresponding contribution to the correlation function. The Fourier amplitudes with n = 0 describe the azimuthally averaged correlations; with n = 1 are not calculated since we do not distinguish between the head and the tail of a bacterium; with n = 2 describe quadrupolar correlations; these are different for two particles placed in a director field side-by-side and for two particles placed head-to-tail. The amplitudes exhibit the following properties: (a) since C(r, φ, θ_l, θ_k) is real, (b) since the system is globally apolar, when n_l + n_k + n is odd, (c) since the system shows a global mirror symmetry, x ↔ −x, then is real when n_l + n_k + n is even. We neglect the terms with any of |n_l|, |n_k|, and |n| being larger than 2, as they are small and do not represent the basic characteristics of the bacterial shapes and orientational correlations.

3.6.1 Density–density and orientation–orientation correlations. The density–density correlation function C_dd(r, φ) corresponds to n_l = n_k = 0, and n = 0; 2. It can be represented as the sum of isotropic and anisotropic terms:

C_dd(r, φ) = C⁽⁰⁾₀₀(r) + C⁽²⁾₀₀(r)cos 2φ, (4)

where the isotropic part describes the radial dependence of the azimuthally averaged density–density correlation function and the azimuthal amplitudeexhibits the azimuthal dependence of C_dd(r, φ), namely, C⁽²⁾₀₀(r) > 0 when C_dd(r, 0) > C_dd(r, π/2) and vice versa.

The isotropic part C⁽⁰⁾₀₀(r) of the density–density correlation in 5c₀ LLC is negative at short distances r < 5 µm and shows a pronounced maximum at r ≈ 7 µm at all times, Fig. 7(a). The negative values of C⁽⁰⁾₀₀(r) are associated with the “steric” repulsions of the bacteria, which are caused not only by the direct contacts but also by the repulsions through hydrodynamic and orientational elastic forces of the N background. The length 7 µm is a typical length of the bacterium's head. We thus associate the maximum of C⁽⁰⁾₀₀(r) with the tendency of bacteria to form chains.⁴² A swimming bacterium leaves a track in the N background that is the least-resistance direction for the bacteria that follow. The chaining effect is evidenced by the behavior of the anisotropic parts of the density–density correlations, Fig. 7(b). The corresponding function C⁽²⁾₀₀(r) is negative for r < 5 µm right after the shear cessation when the bacteria are still oriented along the x-axis. At t = (15–30) s, C⁽²⁾₀₀(r) shows a small positive maximum at r≈10 µm. At t > 30 s, when the defects nucleate and realign the N director by 90°, C⁽²⁾₀₀(r) exhibits a negative minimum at 10 µm < r < 20 µm, clearly seen in Fig. 7(b). The azimuthally dependent density–density correlation C⁽²⁾₀₀(r) shows an interesting far-field behavior, with two positive maxima at r = 40 µm and 80 µm, and two minima at 60 µm and 100 µm. These extrema of C⁽²⁾₀₀(r) correlate clearly with the wavelength λ ≈ 80 µm of the bend undulations in Fig. 4(a).

Fig. 7 The spatial density–density and the orientation–orientation correlation functions as a function of r at different times for 5c₀: (a) the azimuthal average C⁽⁰⁾_dd(r) and (b) the azimuthal amplitude C⁽²⁾_dd(r) of the density–density correlation function, eqn (4); (c) the azimuthal average C⁽⁰⁾_qq(r) and (d) the azimuthal amplitude C⁽²⁾_qq(r) of the quadrupole–quadrupole correlation function, eqn (7).

The quadrupole–quadrupole correlation function is defined as

where is the scalar order parameter of the bacterial array calculated with respect to the shear x-axis. It could be written in a form similar to the density–density correlation function, eqn (4) and (5)

C_qq(r, φ) = C⁽⁰⁾₂₂(r) + C⁽²⁾₂₂(r)cos 2φ, (7)

where C⁽⁰⁾₂₂(r) and C⁽²⁾₂₂(r) are, respectively, the azimuthal average and the azimuthal amplitude of the quadrupole–quadrupole correlation function C_qq(r, φ); C⁽²⁾₂₂(r) > 0 when C_qq(r, 0) > C_qq(r, π/2) and vice versa.

For small separations r<30 µm, C⁽⁰⁾₂₂(r), Fig. 7(c), and C⁽²⁾₂₂(r), Fig. 7(d), follow C⁽⁰⁾₀₀(r) and C⁽²⁾₀₀(r), respectively, Fig. 7(a) and (b), since the bacteria are parallel to each other. The far-field correlations within the time interval 15 s < t < 60 s are nonmonotonous with maxima at r ≈ 40 µm and ≈80 µm, Fig. 7(c), which correlates with the undulation wave in bacterial orientations. The radial modulations of C^(k)₂₂(r) in Fig. 7(c) and (d) are stronger than the similar modulations of C^(k)₀₀(r) in Fig. 7(a) and (b), which means that the modulations are indeed caused by the undulations with the period λ ≈ 80 µm. At later stages, t > 60 s, in topological turbulence with vortices, (i) the radial modulations of C⁽⁰⁾₂₂(r) indicate that the average distance between the vortices from the same pair corresponds to the periodicity of initial undulations, Fig. 7(c); (ii) the variation of anisotropic term preserves the structure typical for earlier stages of the undulations, Fig. 7(d), and corroborates the long-term memory of the pair orientation 〈cos 2φ_d〉, Fig. 7(f).

3.6.2 Velocity–velocity correlation. We calculate the velocity–velocity correlation function C_vv(R) = 〈v̄_⊥y(r)v̄_⊥y(r + R)〉 for the quantity v̄_⊥y, which is the y-component of the perpendicular velocity v̄_⊥. We represent C_vv(R) as two separate parts: one depends on the separation projected onto the shear direction, C_vv(R_x) = 〈v̄_⊥y(x)v̄_⊥y(x + R_x)〉, and the second depends on the separation along the vorticity direction, C_vv(R_y) = 〈v̄_⊥y(y)v̄_⊥y(y + R_y)〉. Fig. 8(a) shows that during the bend instability, t < 60 s, C_vv(R_x) reaches a negative minimum at R_x ≈ 40 µm, caused by an anticorrelation between the flow velocities at the neighboring extrema of the undulation wave, and a maximum at R_x ≈ 80 µm, which corresponds to the full wavelength λ ≈ 80 µm of undulations, Fig. 8(a). As the LLC transitions to turbulence at t > 60 s, C_vv(R_x) and C_vv(R_y) behave similarly, as expected for a system that lost the memory of its original orientation and undulation waves.

Fig. 8 The spatial velocity–velocity correlation function C_vv(R) = 〈v̄_⊥y(r)v̄_⊥y(r + R)〉 of 5c₀ calculated along the shear direction (a) C_vv(R_x), (b) vorticity axis C_vv(R_y) at different times t = 10, 20, 40, 60, and 120 s, and (c) snapshots of the v_x and v_y at different times t = 10, 60, and 120 s.

The correlations of perpendicular components of the velocities in Fig. 8(a) and (b) illustrate the structure of undulations well but fail to describe the turbulence. In Fig. 8(c), we show the snapshots of the y- and x-components of the total velocity field. In the undulation regime, the maps indicate active flows along the y-axis, which are initially weak, Fig. 8(c), and enhance as the turbulence sets in. The behavior of x-component is opposite, with strong flows caused by bacteria at the onset of undulations, followed by random motion at later stages. The snapshots also demonstrate that the typical scale at which the velocities reverse (which can be associated with the size of flow vortices) is about 20–35 µm, somewhat smaller than half-period λ/2 ≈ 40 µm of undulations.

3.7 Energy spectra

We quantify the energy spectra that describe spatial distributions of three quantities, the N elastic energy, kinetic energy, and enstrophy. These quantities are represented as the integrals over the area of each frame. The spectra are calculated by transforming the integrals over space coordinates into integrals over the scalar wavenumber q. The spectra are calculated for selected instants of time, to illustrate the evolution from uniform alignment to undulations and turbulence, Fig. 9 and 10.

Fig. 9 (a) The elastic energy spectrum E_els(q) (units N m⁻²) as a function of the wave vector q measured in µm⁻¹, (b) and (c) two cropped snapshots of the patterns in the 5c₀ LLC at 70 s and 115 s, respectively, taken from Video S6.

Fig. 10 (a) Kinetic energy spectrum K_s(q) (units (µm s)⁻²) and (b) enstrophy spectrum Ω_s(q) (units s⁻²) as a function of the wave vector q (µm⁻¹) at different times for 5c₀. The spectra K_s(q) and Ω_s(q) are averaged over 5 consecutive frames (or 1 s). (c) Vorticity curl v for 5c₀ at different times.

3.7.1 Elastic energy spectrum. E_els(q) is defined through the averaged area density of the elastic energy in real space , where q is the wavenumber, L_xL_y is the frame area, and Q_kl(r) represents the continuous extrapolation of the order parameter Q_kl(r^(j)) = 2n^(j)_kn^(j)_l − δ_kl defined by the bacterial orientations n̂^(j); the Frank elastic constant is estimated³⁹ as K = 10 pN. Using the discrete Fourier transform of Q^(l)_ij, , we calculate the elastic energy spectrum as the weighted sum over q′ inside an annulus of a mean radius q:where w is the weighting function, derived in SI Sections S9, S10, and Fig. S7.

Fig. 9(a) presents the elastic energy spectra E_els(q) for the selected instants of time. For small q, corresponding to r ≈ 250–150 µm, E_els(q) ∝ q^ν at all times, as the log–log plots are linear with a slope ν increasing from 1.0 at the onset of undulations to 3.0 at well-developed turbulence. The most distinguished feature of E_els(q) is a peak with an almost constant wavevector q_p = 2π/r corresponding to r ≈ (70–80) µm, which grows with time. This peak is accompanied by a second-order peak, q_2p = 2q_p, that appears at later times, t > 40 s, Fig. 9(a). The first peak corresponds to the undulations period, λ ≈ 80 µm, Fig. 4(a). At early stages, the undulations are of a single-harmonic sinusoidal shape, and the second-order peak is absent. With time, the undulation amplitude A increases, Fig. 4(b), resulting in a stronger peak at q_p; a somewhat higher λ yields a slight decrease of q_p, Fig. 9(a). The increased A reshapes the sinusoidal bend into a zigzag, causing the second-order peak q_2p = 2q_p. The development of topological turbulence at t = 70–115 s causes a further slight decrease of q_p and broadens the small-q side of the peak, so that the curve E_els(q) has a constant slope for q < q_p.

The snapshots of 5c₀ LLC at 70 s, Fig. 9(b), and 115 s, Fig. 9(c), which correspond to the saturated number of disclinations, Fig. 5(e), reveal that the turbulence patterns preserve the characteristic length scale of the preceding undulations and positions of the E_els(q) peaks. These patterns also exhibit some anisotropy, as the wave of misalignment persists along the shear (vertical) direction in Fig. 9(b), (c) and Video S6; this anisotropy correlates with the evolution of 〈cos 2φ_d〉 in Fig. 5(f).

3.7.2 Kinetic energy spectrum. K_s(q) is defined through the area-averaged velocity field . The kinetic energy spectrum is calculated similarly to the elastic energy (see SI Section S10) asThe kinetic energy spectrum K_s(q) exhibits a smooth monotonous increase with q within the entire range available in the experiment, Fig. 10(a). The set of quasi-linear parallel lines in the log–log plot establishes that K_s(q) ∝ q for all q and times t, with two exceptions discussed later. K_s(q) ∝ q would be a natural result if the two-dimensional kinetic spectrum defined through is homogeneous, i.e., K_s⁽²⁾(q) = K_s⁽²⁾. The last equation implies K_s(q) = 2πK_s⁽²⁾q. The linear dependency K_s(q) ∝ q implies no correlations between the velocities of neighboring bacteria. The two exceptions occur at the onset of disclination nucleation, t ≈ 20 s, and in topological turbulence, t ≈ (60–120) s. At t ≈ 20 s, bacteria produce alternating velocity lanes along the y-axis, with a period ∼80 µm, Fig. 6(a). These lanes yield a strong superlinear behavior in Fig. 10(a) for t ≈ 20 s, with the slope ν increasing from 1.0 to 2.2. Furthermore, at t ≈ (60–120) s, topological turbulence produces spotty patterns of the velocity field, Fig. 6(a), and results in a weak sub-linear behavior with ν decreasing below 1.0.

The K_s(q) behavior of LLC is dramatically different from a 2D dense suspension of B. subtilis in an isotropic fluid, in which case K(q) ∝ q^5/3 in turbulence and reaches a maximum at q that corresponds to the diameter of vortices.^6,22 This difference is caused by the elasticity of the passive N background.

3.7.3 Enstrophy spectrum. Ω_s(q) is defined through the enstrophy area density asThe overall behavior of Ω_s(q) in Fig. 10(b) is very similar to that of the kinetic energy spectrum K_s(q) in Fig. 10(a) and reflects the same underlying physical effects. The only difference is in the values of the slopes, as the enstrophy spectrum generally scales as Ω_s(q) ∝ q^1.5, while at the onset of nucleation, t ≈ 20 s, Ω_s(q) ∝ q³.

The snapshots of the vorticity curl v in Fig. 10(c) demonstrate that the 5c₀ LLC does not show a clear length scale. A bacterium realigns the N director parallel to its body. The realigned director field persists for some time after the bacterium leaves and serves as a “easy channel” for bacteria following the leader.^42,70 The population splits into chain-like “jets”, with bacteria swimming one after the other. These jets form pairs of elongated red and blue areas in Fig. 10(c). Because of chaining, the LLC is very different from an active nematic with an isotropic environment. At 5c₀, the nematic elasticity is still strong enough to withstand the randomizing forces of active flows.

3.8 Bacteria density modulations, clustering, and giant number fluctuation

One key feature of active matter is the presence of strong density fluctuations.^{6,7,24,100–102} Below, we analyze the spatial distribution and temporal evolution of the bacteria number density fluctuations.

The spatial modulation of bacteria number density is observed in the bend undulations regime. Shortly after the shear stops, t = 3 s, the bacteria align mostly along the shear direction; there is no clear modulation of their number density, SI Fig. S8(a). As the amplitude A of undulations increases, Fig. 4(b), the bacteria condense at the shoulders of the wave, avoiding the extrema, at which the radius of curvature is too small to accommodate an elongated bacterium, SI Fig. S8(b). The density is modulated with a wavelength that is two times shorter than the undulation period λ. Similar density modulations have been predicted in simulations of pusher suspensions by Saintillan et al.¹⁰³ To quantify the density fluctuations evolution, we employ the Voronoi tessellation, Fig. 11.

Fig. 11 (a) Voronoi tessellation for four instances of time, t = 6, 20, 60, and 130 s for a 5c₀ LLC. (b)–(d) Voronoi polygons normalized area probability distribution function (PDF). (b) PDF of Voronoi polygons normalized area for 3c₀ at t = 6, 20, 60, and 110 s; (c) PDF for 5c₀ at different times of t = 6, 20, 60, and 130 s; and (d) PDF for 10c₀ at different times of t = 4, 25, 45, and 75 s. The normalization factor σ₀ is defined as the area of the frame divided by the total number of bacteria in the frame.

The Voronoi tessellation involves the following steps. (1) The image of each bacterium within the frame is divided into five equal parts along the long axis.¹⁰⁴ (2) One draws a perpendicular through the middle point of each part of the image. (3) The five middle points are connected to the similar points belonging to other bacteria in the neighborhood, and a perpendicular is drawn through the middle of each connecting segment. (4) The perpendicular lines intersect at the points that become the vertices of polygons known as the Voronoi polygons. We used the “shoelace” method to calculate the area of each polygon. For example, if (x_i, y_i), i = 1, 2, 3 are the coordinates of vertices of a triangular polygon, then the area is calculated as σ = (x₁y₂ − x₂y₁ + x₂y₃ − x₃y₂ + x₃y₁ − x₁y₃)/2. The Voronoi polygon area corresponding to a single bacterium is found as the sum of five polygon areas corresponding to five parts of the bacterium.

Fig. 11 shows the Voronoi polygons of bacteria for 5c₀ LLC. At the early stage, t = 6 s, the Voronoi polygons elongate along the shear direction; they are of approximately the same area, Fig. 11(a). Later times, t = 20, 60, and 130 s, see the emergence of large areas free of bacteria (which produce large Voronoi polygons) and clusters with an elevated concentration of bacteria (small Voronoi polygons), Fig. 11(a).

The analysis of the Voronoi polygons' area distribution reveals bacterial clustering. In Fig. 11(b)–(d), the distribution of areas is shown for a normalized quantity σ/σ₀, where the normalization factor σ₀ is the total area of the frame divided by the number of bacteria in it. For 3c₀, the probability distribution function (PDF) broadens with time, showing the appearance of large empty regions, Fig. 11(b). The peak of the distribution decreases and shifts to smaller areas. For 5c₀, the shift of the distribution peak to the smaller areas is more significant, indicating stronger clustering, Fig. 11(c). For 10c₀, the probability distribution does not change much, Fig. 11(d), since the bacterial density is already very high, so there are very few bacteria-free areas even in the well-developed turbulence.

3.8.1 Giant number fluctuations. The density modulations in undulations and clustering in topological turbulence result in density fluctuations ΔN ∝ 〈N〉^α with α > 1/2, which exceeds the one expected for systems at equilibrium, ΔN ∝ 〈N〉^1/2, where 〈N〉 is the average number of particles and is the standard deviation. For 3c₀, α increases from 0.72 to 0.88 over time, Fig. 12. These values are similar to ones found in other systems. Nishiguchi et al. measured α ≈ 0.63 for a 2D system of B. subtilis in an isotropic fluid,²⁴ while Liu et al. measured α to increase from α ≈ 0.6 to α ≈ 0.83 over 120 s in an active turbulent 3D system of E. coli, also in an isotropic fluid.²⁵ Different models predict different limits for α.^31,105–107 Recent simulations of an active nematic with multiparticle collision dynamics, which resembles an LLC, show that α can grow from 1/2 at low activities to 1 at intermediate activities and then decreases to about 0.8 at the highest activities, which diminish the degree of the nematic order.¹⁰⁸

Fig. 12 (a) Standard deviation ΔN of the number of bacteria as a function of the average number 〈N〉 of bacteria at different times for 3c₀; (b) the fitting exponent α of ΔN ∝ 〈N〉^α as a function of time.

4. Conclusions

Flagellated microswimmers B. subtilis dispersed in a nematic lyotropic chromonic liquid crystal DSCG form a living liquid crystal (LLC). We explore the emergence and development of topological turbulence from a uniformly aligned state produced by shear of a thin (10 µm) layer of LLC. The experimental setup is of a large area to avoid the effect of lateral confinement and eliminates in-plane surface anchoring, which could impose a permanent anisotropy on the system. Upon shear cessation, bacterial activity destabilizes the uniform alignment, producing self-amplifying bend undulations, with an amplitude that initially grows linearly as a function of time, then saturates, while the period λ remains almost constant. The active stress that drives the growth of amplitude is released via nucleating ±1/2 defect pairs at the extrema of the bend. The defects multiply until the undulation pattern is lost and a regime of topological turbulence is established. Increasing the bacteria concentration results in a smaller λ, a higher growth rate of undulations, and an increased number of the ±1/2 defects.

To identify the disclinations and analyze their dynamics in the bacterial population, we propose an approach based on the Delaunay triangulation mesh. The disclinations of strength ±1/2 appear as bound pairs (BP) at the onset of nucleation and develop into isolated disclinations at later stages. The isolated +1/2 disclinations and ±1/2 pairs are initially well aligned along the y-axis orthogonal to the shear plane, Fig. 5(f). As time progresses, the number densities of disclinations grow and saturate. The orientational diffusion of p is much faster than the translational diffusion of +1/2 cores along the vertical x-axis, as shown in Fig. 5(h).

The number density of disclinations in the strongly active 10c₀ LLC shows a nonmonotonic time dependence, Fig. 5(e), caused by the “alignment and activity-triggered annihilation.” Namely, a fast-moving +1/2 disclination of one nucleated pair reaches a −1/2 disclination of a similarly aligned neighbouring pair and annihilates with it, thus reducing the total number of defects.

The velocities of individual bacteria reveal the presence of two components: the actual individual swimming velocity and the velocity of an underlying active flow caused by bacterial interactions. Immediately after the shear cessation, only the first component is present, with a maximum speed ∼20 µm s⁻¹ of bacteria swimming along the shear direction. Emerging bend undulations produce active flows orthogonal to the shear direction. As a result, the bacterial velocity makes a significant angle with the body axis. At the transition from bend to turbulence, most bacteria realign along the active flow. The active flows result from hydrodynamic interactions of the bacteria. Each bacterium produces a pair of axial forces, F⁽ⁱ⁾_prop = −F⁽ⁱ⁾_drag. The active flow could result from adding F⁽ⁱ⁾_prop of a bacterium “i” and F^(j)_prop of a bacterium “j” or adding F⁽ⁱ⁾_drag and F^(j)_prop, or F⁽ⁱ⁾_drag and F^(j)_drag, Fig. 6(b).

The calculated density–density, orientation–orientation, and velocity–velocity correlation functions provide additional information on LLC behavior. The positive maximum of density–density correlations at the typical length of the bacterium body illustrates the tendency of bacteria to form chains; a leading bacterium leaves an easy passage in the N director field for other bacteria to follow. The orientation–orientation correlations show maxima at r ≈ 40 µm and ≈80 µm, which correspond to the wave of undulations; these typical correlation distances are traceable also in the turbulence. The velocity–velocity correlation along the shear direction shows a negative minimum at R_x ≈ 40 µm and a maximum at R_x ≈ 80 µm in the bend instabilities regime, Fig. 8. The negative minimum is caused by antiparallel active flows that develop at the two neighboring bend extrema, while the maximum corresponds to the full wavelength λ ≈ 80 µm of undulations, Fig. 6.

The elastic energy spectra exhibit a peak with an almost constant wavevector that corresponds to the period of undulations. When the sinusoidal undulations transform into zigzag patterns, a second-order peak appears at a doubled value of the wavevector, Fig. 9(a). The kinetic energy and enstrophy spectra associated with the velocity fields exhibit a monotonous increase with the wavevector within the entire range available in the experiment. The kinetic energy spectrum grows linearly with the wavevector, except at the onset of nucleation and in the topological turbulence. In the latter case, one observes a sublinear behavior that is dramatically different from the active nematic formed by B. subtilis swimming in an isotropic environment. Also different are the patterns of vorticity, which in the LLCs do not show a clear length scale, in contrast to their isotropic counterparts. The LLC vorticity is affected by the chaining of swimming bacteria.

Besides the differences in the energy and enstrophy spectra, orientational order of disclinations, chaining of bacteria, there are other unique dynamic properties of the LLC that distinguish them from active nematics in which the self-propelled units are in an isotropic Newtonian environment. A long memory of the imposed macroscopic orientational order is one example; an intriguing, related effect is that in diluted dispersions, c = (0.5–2)c₀, the bacteria continue to swim up and down the pre-imposed shear direction for a very long time, about 3 min after the cessation of flow, Fig. S3. In diluted LLCs, the separation of bacteria is larger than the range of their hydrodynamic interactions, which explains the long-lasting orientational order.

To conclude, we analyzed the LLCs in which the active nematic part, the dispersion of swimming bacteria, is controlled by an orientationally ordered background of the passive nematic, representing a lyotropic chromonic liquid crystal DSCG. The 10 µm scale of the bacterial bodies allows us to image and analyze the relevant parameters with great detail by optical microscopy. The tangential anchoring of the passive nematic at the bacterial bodies allows us to connect the passive director field to the bacterial orientations and to calculate important parameters characterizing the cascade of orientational disorder, first through bend undulations and then through nucleation of disclinations and turbulence. Combined with the previously determined material properties of DSCG, such as elastic constants and viscosities,^81,82 resolved for splay, bend, and twist of the director, DSCG behavior in confinement,^80,109–111 the results complete a comprehensive description of the experimentally assessable type of active matter in which the self-propelled microswimmers move in a viscoelastic orientationally ordered medium.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00744e.

Data for this article, including binary videos for the experimental data and the codes used for the analysis are available at the GitHub repository at https://github.com/HendBaza/Liquid_Crystals_based_Active_Nematics-.git.

Acknowledgements

The authors thank L. Reichel and A. Buccini for their valuable discussions and insightful comments. The work was supported by NSF grants DMS-1729509 and DMR-2341830.

References

1. S. Ramaswamy, Annu. Rev. Condens. Matter Phys., 2010, 1, 323–345.
2. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
3. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246.
4. M. Bar, R. Grossmann, S. Heidenreich and F. Peruani, Annu. Rev. Condens. Matter Phys., 2020, 11, 441–466.
5. S. P. Thampi, Curr. Opin. Colloid Interface Sci., 2022, 61, 101613.
6. R. Alert, J. Casademunt and J. F. Joanny, Annu. Rev. Condens. Matter Phys., 2022, 13, 143–170.
7. S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti and V. Vitelli, Nat. Rev. Phys., 2022, 4, 380–398.
8. T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431–434.
9. J. Hardouin, R. Hughes, A. Doostmohammadi, J. Laurent, T. Lopez-Leon, J. M. Yeomans, J. Ignés-Mullol and F. Sagués, Commun. Phys., 2019, 2, 121.
10. A. Senoussi, S. Kashida, R. Voituriez, J. C. Galas, A. Maitra and A. Estevez-Torres, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 22464–22470.
11. P. Chandrakar, M. Varghese, S. Aghvami, A. Baskaran, Z. Dogic and G. Duclos, Phys. Rev. Lett., 2020, 125, 257801.
12. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216.
13. K. Kawaguchi, R. Kageyama and M. Sano, Nature, 2017, 545, 327–331.
14. S. Z. Lin, D. P. Bi, B. Li and X. Q. Feng, J. R.Soc., Interface, 2019, 16, 20190258.
15. T. Turiv, J. Krieger, G. Babakhanova, H. Yu, S. V. Shiyanovskii, Q.-H. Wei, M.-H. Kim and O. D. Lavrentovich, Sci. Adv., 2020, 6, eaaz6485.
16. S. Z. Lin, W. Y. Zhang, D. P. Bi, B. Li and X. Q. Feng, Commun. Phys., 2021, 4, 21.
17. L. Balasubramaniam, R. M. Mège and B. Ladoux, Curr. Opin. Genet. Dev., 2022, 73, 101897.
18. Y. M. Luo, M. Y. Gu, M. Park, X. Y. Fang, Y. Kwon, J. M. Urueña, J. R. de Alaniz, M. E. Helgeson, C. M. Marchetti and M. T. Valentine, J. R.Soc., Interface, 2023, 20, 20230160.
19. H. Li, X. Q. Shi, M. J. Huang, X. Chen, M. F. Xiao, C. L. Liu, H. Chaté and H. P. Zhang, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 777–785.
20. A. Be'er and G. Ariel, Mov. Ecol., 2019, 7, 9.
21. V. Yashunsky, D. J. G. Pearce, G. Ariel and A. Be'er, Soft Matter, 2024, 20, 4237–4245.
22. H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Lowen and J. M. Yeomans, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 14308–14313.
23. S. D. Ryan, A. Sokolov, L. Berlyand and I. S. Aranson, New J. Phys., 2013, 15, 105021.
24. D. Nishiguchi, K. H. Nagai, H. Chaté and M. Sano, Phys. Rev. E, 2017, 95, 020601.
25. Z. Y. Liu, W. Zeng, X. L. Ma and X. Cheng, Soft Matter, 2021, 17, 10806–10817.
26. Y. Peng, Z. Y. Liu and X. Cheng, Sci. Adv., 2021, 7, eabd1240.
27. I. S. Aranson, Rep. Prog. Phys., 2022, 85, 076601.
28. B. Martinez-Prat, J. Ignés-Mullol, J. Casademunt and F. Sagués, Nat. Phys., 2019, 15, 362–366.
29. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.
30. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E, 2007, 76, 031921.
31. D. Saintillan and M. J. Shelley, J. R.Soc., Interface, 2012, 9, 571–585.
32. A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2012, 109, 248109.
33. K. Beppu and J. V. I. Timonen, Commun. Phys., 2024, 7, 216.
34. C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein and J. O. Kessler, Phys. Rev. Lett., 2004, 93, 098103.
35. E. Lushi, H. Wioland and R. E. Goldstein, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 9733–9738.
36. H. P. Zhang, A. Be'er, E. L. Florin and H. L. Swinney, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 13626–13630.
37. H. R. Xu, J. Dauparas, D. Das, E. Lauga and Y. L. Wu, Nat. Commun., 2019, 10, 1792.
38. A. Sokolov, I. S. Aranson, J. O. Kessler and R. E. Goldstein, Phys. Rev. Lett., 2007, 98, 158102.
39. S. Zhou, A. Sokolov, O. D. Lavrentovich and I. S. Aranson, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 1265–1270.
40. C. H. Peng, T. Turiv, Y. B. Guo, Q. H. Wei and O. D. Lavrentovich, Science, 2016, 354, 882–885.
41. M. M. Genkin, A. Sokolov, O. D. Lavrentovich and I. S. Aranson, Phys. Rev. X, 2017, 7, 011029.
42. S. Zhou, O. Tovkach, D. Golovaty, A. Sokolov, I. S. Aranson and O. D. Lavrentovich, New J. Phys., 2017, 19, 055006.
43. M. M. Genkin, A. Sokolov and I. S. Aranson, New J. Phys., 2018, 20, 043027.
44. C. Conklin, J. Viñals and O. T. Valls, Soft Matter, 2018, 14, 4641–4648.
45. A. Sokolov, A. Mozaffari, R. Zhang, J. J. de Pablo and A. Snezhko, Phys. Rev. X, 2019, 9, 031014.
46. T. Turiv, R. Koizumi, K. Thijssen, M. M. Genkin, H. Yu, C. Peng, Q.-H. Wei, J. M. Yeomans, I. A. Aranson, A. Doostmohammadi and O. D. Lavrentovich, Nat. Phys., 2020, 16, 481–487.
47. R. Koizumi, T. Turiv, M. M. Genkin, R. J. Lastowski, H. Yu, I. Chaganava, Q.-H. Wei, I. S. Aranson and O. D. Lavrentovich, Phys. Rev. Res., 2020, 2, 033060.
48. M. A. Boulé, S. Rainville and T. Galstian, Mol. Cryst. Liq. Cryst., 2020, 712, 10–23.
49. N. Figueroa-Morales, M. M. Genkin, A. Sokolov and I. S. Aranson, Commun. Phys., 2022, 5, 301.
50. A. Vats, P. K. Yadav, V. Banerjee and S. Puri, Phys. Rev. E, 2023, 108, 024701.
51. A. Vats, V. Banerjee and S. Puri, Phys. Rev. E, 2024, 110, 034701.
52. J. Li, L. Ohm and S. E. Spagnolie, Phys. Rev. Fluids, 2025, 10, 083301.
53. Z. Y. Mou, Y. Li, Z. H. You and R. Zhang, Phys. Rev. E, 2025, 111, 065410.
54. O. D. Lavrentovich, Liq. Cryst. Rev., 2020, 8, 59–129.
55. W. S. Kim, J. H. Im, H. Kim, J. K. Choi, Y. Choi and Y. K. Kim, Adv. Mater., 2023, 35, 202204275.
56. J. Lydon, Liq. Cryst., 2011, 38, 1663–1681.
57. H. S. Park and O. D. Lavrentovich, in Liquid Crystals Beyond Displays: Chemistry, Physics, and Applications, ed. Q. Li, John Wiley & Sons, Hoboken, New Jersey, 2012, ch. 14, pp. 449–484.
58. P. J. Collings, J. N. Goldstein, E. J. Hamilton, B. R. Mercado, K. J. Nieser and M. H. Regan, Liq. Cryst. Rev., 2015, 3, 1–27.
59. C. J. Woolverton, E. Gustely, L. Li and O. D. Lavrentovich, Liq. Cryst., 2005, 32, 417–423.
60. L. Gaojin, L. Eric and A. M. Ardekani, J. Non-Newtonian Fluid Mech., 2021, 297, 104655.
61. S. E. Spagnolie and P. T. Underhill, Annu. Rev. Condens. Matter Phys., 2023, 14, 381–415.
62. C. Viney, Biorheology, 1999, 36, 319–323.
63. H. Chi, A. Gavrikov, L. Berlyand and I. S. Aranson, Commun. Phys., 2022, 5, 274.
64. I. I. Smalyukh, J. Butler, J. D. Shrout, M. R. Parsek and G. C. L. Wong, Phys. Rev. E, 2008, 78, 030701.
65. R. Hartmann, P. K. Singh, P. Pearce, R. Mok, B. Song, F. Díaz-Pascual, J. Dunkel and K. Drescher, Nat. Phys., 2019, 15, 251–256.
66. T. G. J. Chandler and S. E. Spagnolie, Phys. Rev. Fluids, 2024, 9, 110511.
67. D. Tampion and R. A. Gibbons, Nature, 1962, 194, 381.
68. S. Liu, S. Shankar, M. C. Marchetti and Y. L. Wu, Nature, 2021, 590, 80–84.
69. W. T. Liao and I. S. Aranson, PNAS Nexus, 2023, 2, pgad291.
70. P. C. Mushenheim, R. R. Trivedi, H. H. Tuson, D. B. Weibel and N. L. Abbott, Soft Matter, 2014, 10, 88–95.
71. P. C. Mushenheim, R. R. Trivedi, D. B. Weibel and N. L. Abbott, Biophys. J., 2014, 107, 255–265.
72. P. C. Mushenheim, R. R. Trivedi, S. S. Roy, M. S. Arnold, D. B. Weibel and N. L. Abbott, Soft Matter, 2015, 11, 6821–6831.
73. R. R. Trivedi, R. Maeda, N. L. Abbott, S. E. Spagnolie and D. B. Weibel, Soft Matter, 2015, 11, 8404–8408.
74. A. Sokolov, S. Zhou, O. D. Lavrentovich and I. S. Aranson, Phys. Rev. E, 2015, 91, 013009.
75. M. Goral, E. Clement, T. Darnige, T. Lopez-Leon and A. Lindner, Interface Focus, 2022, 12, 20220039.
76. A. G. Prabhune, A. S. García-Gordillo, I. S. Aranson, T. R. Powers and N. Figueroa-Morales, PRX Life, 2024, 2, 033004.
77. C. H. Peng, Y. B. Guo, T. Turiv, M. Jiang, Q. H. Wei and O. D. Lavrentovich, Adv. Mater., 2017, 29, 1606112.
78. S. Mandal, T. J. Mason, A. C. Croft and M. G. Mazza, Phys. Rev. Lett., 2025, 134, 128302.
79. B. Gautam and J. S. Lintuvuori, Phys. Rev. Lett., 2024, 132, 238301.
80. L. Tortora, H.-S. Park, S.-W. Kang, V. Savaryn, S.-H. Hong, K. Kaznatcheev, D. Finotello, S. Sprunt, S. Kumar and O. D. Lavrentovich, Soft Matter, 2010, 6, 4157–4167.
81. S. Zhou, K. Neupane, Y. A. Nastishin, A. R. Baldwin, S. V. Shiyanovskii, O. D. Lavrentovich and S. Sprunt, Soft Matter, 2014, 10, 6571–6581.
82. S. Zhou, A. J. Cervenka and O. D. Lavrentovich, Phys. Rev. E, 2014, 90, 042505.
83. J. J. Yu, L. F. Chen, G. Y. Li, Y. R. Li, Y. Z. Huang, M. Bake and Z. Tian, J. Mol. Liq., 2021, 344, 117756.
84. H. Baza, T. Turiv, B.-X. Li, R. Li, B. M. Yavit, M. Fukuto and O. D. Lavrentovich, Soft Matter, 2020, 16, 8565–8576.
85. M. Shribak, Sci. Rep., 2015, 5, 17340.
86. M. Rajabi, O. Lavrentovich and M. Shribak, Liq. Cryst., 2023, 50, 181–190.
87. J. L. Hu, M. C. Yang, G. Gompper and R. G. Winkler, Soft Matter, 2015, 11, 7867–7876.
88. Q. Zhang, R. Zhang, B. L. Ge, Z. Yaqoob, P. T. C. So and I. Bischofberger, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2108361118.
89. Marcos, H. C. Fu, T. R. Powers and R. Stocker, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4780–4785.
90. J. T. Gleeson, R. G. Larson, D. W. Mead, G. Kiss and P. E. Cladis, Liq. Cryst., 1992, 11, 341–364.
91. J. Vermant, P. Moldenaers, J. Mewis and S. J. Picken, J. Rheol., 1994, 38, 1571–1589.
92. P. Harrison and P. Navard, Rheol. Acta, 1999, 38, 569–593.
93. M. P. Lettinga, Z. Dogic, H. Wang and J. Vermant, Langmuir, 2005, 21, 8048–8057.
94. J. K. G. Dhont, M. P. Lettinga, Z. Dogic, T. A. J. Lenstra, H. Wang, S. Rathgeber, P. Carletto, L. Willner, H. Frielinghaus and P. Lindner, Faraday Discuss., 2003, 123, 157–172.
95. L. G. Wang and Y. Huang, Liq. Cryst., 2003, 30, 1129–1137.
96. J. P. Canejo, N. Monge, C. Echeverria, S. N. Fernandes and M. H. Godinho, Liq. Cryst. Rev., 2017, 5, 86–110.
97. T. Ishikawa and O. D. Lavrentovich, in Defects in Liquid Crystals: Computer Simulations, Theory and Experiments, ed. O. D. Lavrentovich, P. Pasini, C. Zannoni and S. Zumer, Kluwer Academic Publishers, Dordrecht, 2001, pp. 271–301.
98. L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, Phys. Rev. Lett., 2013, 110, 228101.
99. R. A. Simha and S. Ramaswamy, Phys. A, 2002, 306, 262–269.
100. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
101. V. Narayan, S. Ramaswamy and N. Menon, Science, 2007, 317, 105–108.
102. H. Chaté, F. Ginelli and R. Montagne, Phys. Rev. Lett., 2006, 96, 180602.
103. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2008, 100, 178103.
104. X.-Q. Shi and H. Chaté, arXiv, 2018, preprint, arXiv:1807.00294 DOI:10.48550/arXiv.1807.00294.
105. F. Ginelli, F. Peruani, M. Bär and H. Chaté, Phys. Rev. Lett., 2010, 104, 184502.
106. S. Dey, D. Das and R. Rajesh, Phys. Rev. Lett., 2012, 108, 238001.
107. S. Ngo, A. Peshkov, I. S. Aranson, E. Bertin, F. Ginelli and H. Chaté, Phys. Rev. Lett., 2014, 113, 038302.
108. T. Kozhukhov and T. N. Shendruk, Sci. Adv., 2022, 8, eabo5788.
109. L. Tortora and O. D. Lavrentovich, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 5163–5168.
110. J. Jeong, Z. S. Davidson, P. J. Collings, T. C. Lubensky and A. G. Yodh, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 1742–1747.
111. R. Koizumi, D. Golovaty, A. Alqarni, B. X. Li, P. J. Sternberg and O. D. Lavrentovich, Sci. Adv., 2023, 9, eadf3385.

---