Controlled route to active turbulence: filling an activity spot with topological defects

Abstract

We present a controlled route to active turbulence in an active paranematic fluid, i.e., a suspension of rods that exhibit nematic order only under extensile activity. To this end, we introduce a spot of radius r with non-zero activity, embedded in an otherwise passive fluid. Due to the open boundary, defects can enter and leave the spot. As r increases starting from the nematic coherence length, we first observe paranematic order with a uniform director field, then transient +½ topological defects, followed by spiral or swirling pairs of +½ defects. Additional defects progressively enter until bulk active turbulence occurs. While positive and negative defect charges grow with the square of the spot radius r, the total charge only increases linearly in r. This hints to a length along the rim of the spot, comparable to the active length, so that activity can induce the director distortions needed for a defect to enter. In addition, the extensional active flow realizes active anchoring at the rim, which establishes a baseline charge of +1. Two dynamic regimes mark the progression toward bulk turbulence. The enstrophy rises sharply when the spot allows the stable circling motion of the two +½ defects, and the finite-time Lyapunov exponent, characterizing the chaotic flow pattern, jumps to a non-zero value, when a third +½ defect enters the spot noticeably. For large radii, both measures approach their bulk-turbulence values.

---

1 Introduction

Active matter refers to a class of materials in which individual subunits consume locally available energy, drive the system out of equilibrium, and create coherent structures on length scales larger than those of the individual units.^1–5 In particular, bacterial suspensions,^6,7 microtubule-kinesin solutions,^8–10 and cell monolayers^11,12 can exhibit active turbulence. Here, complex flows occur characterized by chaotically moving fluid structures with strong vorticity that show local nematic ordering and chaotically moving +½ topological defects. They are continuously created and annihilated in pairs. In our previous article we investigated the spatial control of active turbulence¹³ by considering the active species as self-propelling rods and using Doi's hydrodynamic equations for a semidilute rod solution,¹⁴ augmented with an active-stress term.¹⁵ These equations enable us to study activity-induced nematic order while the passive fluid remains isotropic. We call this the active paranematic state,^13,16,17 relevant to a number of experimental realizations.^18–20 In contrast to the widely used active nematic models,^1,21–23 which assume an ordered equilibrium nematic state, our approach captures systems that are isotropic in the absence of activity. However, we stress that our dynamic equations are essentially the same as in the other models but we leave out higher-order terms in the Landau–de Gennes free energy that generate nematic order already in the passive case.

Controlling the dynamics of active fluids for designing functional materials, organizing flow patterns, and generating regularized flows is a central goal of the field.^5,13,24,25 Several mechanisms have been proposed, including geometric confinement,^10,26–29 substrate friction,^30,31 optimal-control policies,^32–34 and directly patterning activity.^13,35,36 For example, experiments show that confining active nematics within circular domains of appropriate size transforms chaotic defect motion into regular patterns of two orbiting +½ defects, with frequent unbinding of defect pairs at the boundary.²⁸ Likewise, theory predicts the formation of a circulating pair of +½ defects under such confinement.³⁷ More recently, patterning activity has attracted attention as a means to control the flow of active fluids.^13,35,36 Advances in designing photo-activatable bacteria^38,39 and motor-proteins^40–43 now make it possible to realize locally varying activity in an experiment. For example, a lattice of inactivity spots reorganizes otherwise turbulent flow into a multi-lane flow state, where neighboring flow lanes with alternating flow directions are separated by a street of vortices.¹³ Activity patterning can also steer, trap, and segregate topological defects.^36,40,44 For example, introducing a circular active region within a passive nematics, which is confined in a circular domain, allows to regulate defect patterns.⁴⁵

In this article we study a circular activity spot that creates an active nematic fluid region embedded in an otherwise passive isotropic fluid and vary its radius [Fig. 1(a)]. In contrast to previous works, the boundary of the spot imposes neither geometric confinement nor surface anchoring, so that topological defects [Fig. 1(b)] are constantly able to enter and leave the spot. We demonstrate that the total topological charge inside the spot consists of a charge +1 due to active anchoring at the spot rim, which is determined by the flow-driven director alignment. Additional +½ defects can enter the spot with increasing radius since activity provides the necessary energy input for the required director deformations. For spot radi around the nematic coherence length λ_n, the director field in the spot is uniform. But then for spot radii 3–5 times larger than λ_n, a single spiral defect or a stable pair of swirling +½ defects occur. With further increasing radius progressively more topological defects enter and ultimately bulk active turbulence is reached, which we quantify by the enstrophy. Furthermore, beyond the regime of two swirling defects, the flow field and thereby the defect motion becomes chaotic, as the finite-time Lyapunov exponent (FTLE) shows.

Fig. 1 (a) Schematic of our system setup using an active paranematic fluid. We consider a circular spot of radius r where activity is negative, W = −3.25, such that the fluid is in an active paranematic state, while in the surrounding region the activity is zero, W = 0. (b) Possible topological defects in an active nematic. Left: The +½ defect, where the director rotates by an angle π when circling in an anticlockwise direction around the defect core, has a polar structure. Therefore, it is advected by the self-generated fluid flow. Right: The +½ defect is nonpolar and, therefore, does not move. The dark grey region indicates the defect core, the diameter of which is determined by the nematic coherence length λ_n.

This article is structured as follows. Section 2 explains the theoretical modeling including the active paranematic state and system parameters. Section 3 presents the results from our simulations for increasing spot radius. We analyze the total topological charge in the spot, active flow-induced anchoring at its rim, the gradual transition to bulk active turbulence, and the finite-time Lyapunov exponent. We close with a summary and conclusions.

2 Theoretical modeling

We first introduce the continuum model of an active paranematic fluid, formulate nondimensional equations and comment on our numerical method. Then, we shortly review some facts on the active paranematic order and on active turbulence from our previous article,¹³ and, finally, mention our system parameters as well as some details for topological defects.

2.1 Continuum model of an active paranematic fluid

In contrast to standard models of active nematics,^1,21–23 we consider a continuum model for a two-dimensional solution of self-propelling rod-shaped particles that do not exhibit nematic order in the passive case. So, we leave out higher-order terms in the Landau–de Gennes free energy that generate nematic order already in passive nematics. In our previous work, we used this model to study activity-induced nematic ordering.¹³ Thus, we rely on the Doi–Edwards modeling of a suspension of rod-shaped particles in the semidilute regime, where rods interact through excluded-volume encounters but at concentrations below the isotropic–nematic phase transition. The relevant continuum theory follows by coarse-graining the Smoluchowski equation for the full orientational and positional distribution function of rigid rods.¹⁴ In the momentum equation, we add an active stress tensor¹⁵ that generates local flows, which in turn align the rods to create nematic order.¹³ Thereby, a feedback loop is closed where spatial variations of the nematic order generate stresses that drive further flow. Since nematic order is caused by flow, we call our system an active paranematic fluid.

The continuum equations are formulated in terms of the velocity field of fluid flow, u, and the alignment or order-parameter tensor Q. The latter is symmetric and traceless and quantifies both the preferred orientation and the degree of alignment of rod-shaped particles. In two dimensions, the eigenvector n (|n| = 1) associated with the largest eigenvalue S is called the director, and it gives the mean orientation of the rods. The scalar order parameter S quantifies the degree of orientational ordering and ranges from 0 in a fully isotropic state to 1 for perfect alignment. With both quantities the alignment tensor takes the form , where 1 is the identity tensor.

In Doi's theory infinitely thin rods are assumed and the dynamics of the alignment tensor Q is described by¹⁴

where L = ∇ ⊗ u is the velocity gradient and A·B = A_ijB_ij. The first two terms on the LHS of eqn (1) give the material time derivative of Q. The following two terms account for the reorientation of the tensor Q due to velocity gradients. Decomposing L into its antisymmetric part, the vorticity tensor Ω = ½(L − L^T), and its symmetric part, the strain-rate tensor E = ½(L + L^T), shows that they contribute with equal weight to reorientation. This corresponds to an alignment parameter set to 1, which for other values describes rods with a finite aspect ratio.^46,47 The last term on the LHS of eqn (1) guarantees the zero trace of Q, since it cancels any nonzero trace of the third and fourth term. According to the first term on the RHS of eqn (1), the order parameter relaxes to zero since we do not include higher-order terms in Q. Thus, the passive rod solution does not exhibit any nematic order. The relaxation coefficient cD_r is the effective rotational diffusion constant in the semidilute regime. Here, D_r is the rotational diffusion coefficient of an isolated rod, and c < 1 is a reduction factor that accounts for the hindered rotational motion due to surrounding rods. The last term penalizes any non-uniform nematic order. The coefficient D₂, which has the dimension of a diffusivity, is mainly determined by steric interactions of the rods.

In Doi's theory the dynamic equation for the velocity field u corresponds to the momentum balance of an incompressible fluid with ∇·u = 0 and¹⁴

where p = P/ρ with P the pressure and ρ the constant density. The three terms in the square brackets on the RHS of eqn (2) define the stress tensor. It consists of the conventional shear stresses with the kinematic shear viscosity ν, some anisotropic contribution due to nematic ordering, where η_p is proportional to the volume fraction of the rods, and active stresses. The latter are generated by the self-propelled rods that act with force dipoles oriented along the rods on the fluid. Thus, in the continuum description the active stress tensor is proportional to the alignment tensor and its strength, named activity W, is proportional to the force dipole and the number density of the rods.¹⁵ The sign of W determines the type of stresses generated by the active particles: W < 0 for extensile (pushers) and W > 0 for contractile (pullers).

2.2 Nondimensional equations and numerical method

We nondimensionalize the equations by defining the characteristic length and time as and , with the final nondimensional form provided in the Appendix. The resulting dynamic eqn (12) and (13) are solved numerically using a pseudo-spectral method for spatial discretization,^48,49 with time integration carried out by the fourth-order exponential time differencing Runge–Kutta scheme (HochOst4),⁵⁰ and the incompressibility condition enforced via the projection method.⁵¹ Simulations are initialized with Q = 0 and u = 0, corresponding to an isotropic orientational distribution of the active rods and a quiescent flow field.

2.3 Onset of active paranematic order and active turbulence

We shortly review results from our earlier publication.¹³ Linear stability analysis applied to the isotropic resting fluid, Q = 0 and u = 0, shows that this state becomes unstable for W < W_c = −8c.¹³ This indicates that a system of pusher-type rods develops nematic order together with a nonzero velocity field. Since the nematic order arises purely from activity, we call it an active paranematic state.

The numerical solution of the full equations further shows that for W < W_c = −8c the system develops active turbulence accompanied by the continuous creation and annihilation of +½ defects, which move irregularly throughout the system. A complete analysis of active turbulence is provided in (ref. 13), where we use various measures including the mean enstrophy, which quantifies the intensity of vortices generated in the system.

From the linear stability analysis, the wavelength λ_max of the most unstable mode is determined as

The last equation is valid for D₂ ≪ 1, which applies in our case, and away from the instability at W = −8c. We use λ_max as the active length λ_a of the system, which characterizes the typical vortex size and the mean distance between topological defects as we showed in (ref. 13).

The instability of the isotropic resting fluid at W < W_c = −8c can also be understood by an alternative argument. Retaining only the linear terms in the nondimensional dynamic eqn (13) for the Q tensor and neglecting the Laplacian term on sufficiently large lengths, gives Q = E/4c. This means that extensional flow along the principal axis of E generates nematic order along this axis. Now, substituting the result for Q into the active stress tensor in the momentum balance equation, one realizes that activity contributes to an effective shear viscosity ν_eff/ν = 1 + W/8c. If W < −8c, thus exactly at the onset of the instability, the effective viscosity becomes negative and thereby generates spontaneous shear flow. This produces nematic ordering, which in turn feeds back into the active stress, thus generating stronger flows in a positive feedback loop.

2.4 System parameters and topological defects

As described in Fig. 1, we consider a circular spot of radius r with nonzero activity surrounded by an otherwise inactive fluid. We always choose an activity W = −3.25 inside the spot and W = 0 outside. Thus, inside the spot paranematic order is generated. As we will see, the rim of the activity spot acts as an open boundary, allowing topological defects to enter and leave the spot.

According to ref. 52, we have implemented the following procedure to identify topological defects. First, we locate the position of a defect from the intersection of the two zero-contour lines of the scalar fields Q_xx and Q_xy. Then, for each defect the topological charge is determined from the rotation of the nematic director, when proceeding along a closed contour line enclosing the defect core in the counterclockwise direction. In practice, we use a square contour of linear size s, chosen to be smaller than the mean distance between neighboring defects. A rotation of π corresponds to a +½ defect, while a rotation of −π corresponds to a −½ defect.⁵³

As Fig. 1(b) shows, a +½ defect has a polar (comet-like) structure that can be characterized by a polarization vector p = (∇·Q(x_i))/|∇·Q(x_i)|, where x_i is the core position.⁵⁴ This polarity drives a self-generated flow with finite velocity at the defect core. So, the defect is advected by the flow it creates and self-propels along p for extensile activity, while in contractile systems the direction is reversed. In contrast, the −½ defect is nonpolar. It generates a symmetric flow field with zero velocity at the core and, therefore, the defect does not self-propel.^54–56

The size of a defects core, where the scalar order parameter develops from zero to its bulk value, is given by the nematic coherence length , which is another characteristic length scale of the system. Generally, it measures the distance over which the nematic order parameter assumes its bulk value, when fixed at some boundary, for example, the center of the defect core.

For all simulations, we set η_p = 3 and c = 0.1, which places the system in the semidilute regime of rod-like solutions, and we further choose D₂ = 10⁻³. More details are given in our previous article.¹³ For this parameter set, the characteristic lengths evaluate to λ_n = 0.05 (nematic coherence length) and λ_a = 1.0 (active length). The activity spot is embedded in a square domain of linear size L with periodic boundary conditions. The domain size L is increased with the spot radius r and the spectral resolution (number of retained wave vectors n_k) is increased accordingly. We have verified that the results are independent of both L and n_k.

We explore spot radii in the interval r ∈ [0.1, 20]. So, on the lower side the radii are comparable to λ_n and they become much larger than λ_a. With this setup in place, we now examine how the defect dynamics in the activity spots varies with r.

3 Results

First, we demonstrate and quantify how topological defects in the director orientation fill the activity spot with increasing spot radius. Second, we show that the total topological defect charge in the spot is determined by two contributions: due to active anchoring and due to activity-induced director deformations that allow defects to enter the spot. Finally, we measure the chaosity in fluid flow by the finite-time Lyapunov exponent.

3.1 Route to active turbulence: filling in topological defects

To illustrate how topological defects fill up the activity spot with increasing spot radius, we show in Fig. 2 snapshots of the dynamic director field configuration in steady state together with the color coded scalar order parameter S for several spot radii r = 0.10–2.5. In addition, we also indicate the motile +1/2 defects (red dots with arrow) and the non-motile −1/2 defects (blue dots). For all snapshots we provide videos in the SI.

Fig. 2 Snapshots of the dynamic nematic director field in steady state and associated topological defects in activity spots with different radii ranging from r = 0.10 to 2.5. The background colour represents the scalar order parameter S, black lines show the nematic director field, and the rim of the spot is indicated by the gray dashed line. +½ and −½ defects are marked by red and blue dots, respectively, with a red arrow indicating the orientation of each +½ defect. (a) r = 0.10: uniform director field, no defects; (b) r = 0.13: temporary +½ defect enters the spot; (c) r = 0.15: two +½ defects form a stable spiral defect; (d) r = 0.20: two +½ defects swirl with constant distance; (e) r = 0.25: third temporary +½ defect enters; (f) r = 0.40: additional +½ defects occur and −½ defects appear; (g) r = 1.05: increased population of both +½ and −½ defects; (h) r = 2.5: bulk-turbulent behaviour. The snapshots (b)–(h) are taken from Videos S1–S7 in the SI. The maximum values S_max of the scalar order parameter for panels (a)–(h) are in order: S_max = 0.22, 0.31, 0.27, 0.43, 0.49, 0.56, 0.63, 0.64.

For a spot with radius r = 0.10 comparable to the nematic coherence length (λ_n ≈ 0.05) [see Fig. 2(a)], a uniform director field occurs within the spot without any defects since activity cannot distort the director field on such small scales. With the radius slightly increased to r = 0.13 [see Fig. 2(b)], a +½ defect repeatedly enters the spot at the rim, deforms the director field, and exits again, so that the director field relaxes back to its uniform orientation. This is nicely illustrated in Video S1 in the SI. Thus, the spot is still too small to permanently host defects.

This changes for spots with radii r = 0.15 and r = 0.20 [see Fig. 2(c) and (d) as well as related Videos S2 and S3, where two +½ defects appear inside the spot. For r = 0.15, they approach each other and form a stable spiral defect with total charge +1. For the larger radius r = 0.20, the two +½ defects separate, orient anti-parallel, and circle around each other while maintaining a constant separation as Video S3 shows. We discuss the velocity fields for these two cases in the appendix in Section B.

As the radius is increased to r = 0.25 [see Fig. 2(e) and Video S4], a third +½ defect intermittently enters the spot, and the system alternates between two configurations: (i) episodes in which two defects remain in the spot and perform circulating motion similar to what is observed at r = 0.20; and (ii) episodes in which the third defect enters at the rim, approaches one of the existing defects, and engages in circular motion, while the other defect leaves the spot, which returns the system to the two-defect state. As we have an open boundary at the rim of the spot, where the nematic director is not anchored by some surface energy, defects can enter and leave the spot, which causes the number of defects to fluctuate in the spot. Fig. 3(a) shows the time course of defect counts for t > 50, after the system has reached its dynamic steady state. The red and blue curves record the respective numbers of +½ and −½ defects. The number of +½ defects fluctuates between two and three, while −½ defects are not observed.

Fig. 3 Time evolution of the number of topological defects moving in activity spots of three representative radii, after the system has reached its steady state at t > 50. The red and blue curves give the respective count of +½ and −½ defects.

For a larger spot radius of r = 0.40 [see Fig. 2(f) and Video S5], additional +½ defects enter, causing the instantaneous count of +½ defects to alternate between two and four, and the defect dynamics becomes more irregular. Furthermore, −½ defects also appear in the spot, but only for a short period [see Fig. 3(b)]. They either enter and leave through the rim or are created inside the spot together with a +½ defect and then annihilate immediately. Further increasing the spot radius to r = 1.05 [ see Fig. 2(g) and Video S6], raises the number of defects of both signs within the spot. Eventually, at r = 2.5 the spot is filled with many defects and displays bulk active turbulence, as we show in Section 3.1.1. Video S7 illustrates the defect motion. Time variations of the defect count of +½ and −½ defects [see Fig. 3(c)] rise and fall together indicating that they are created and annihilated in pairs within the spot as in bulk turbulence.^13,23 But defects still enter separately through the rim as observed for small spot sizes. As a consequence, the total positive charge of the spot grows continuously. We will discuss this in more detail in Section 3.1.2.

3.1.1 Enstrophy. To quantify the transition towards bulk turbulence with increasing spot radius, we use the enstrophy per unit area,where the integral goes over the area A of the activity spot. It serves as a measure of the total vorticity generated within the flow field. In Fig. 4 we plot Ω versus spot radius r. When the two +½ defects start to perform their swirling motion at r ≈ 0.15 (first vertical dashed line), the enstrophy jumps from small to large values. It shows some oscillatory behavior, when the two defects separate from each other and swirl with constant separation, which varies with r. The curve becomes smoother beyond the second dashed line, where the third defect starts to enter, and then continuously increases, when even more defects enter the spot (third dashed lines). Ultimately, at r ⪆ 2.5 the curve reaches a plateau. It corresponds to the enstrophy value of bulk turbulence, which we determined in ref. 13.

Fig. 4 Enstrophy per unit area, Ω, inside an activity spot plotted versus its radius r. The red dashed line shows the bulk value of Ω for active turbulence. The first, second, and third vertical dashed lines indicate, respectively, the first noticeable appearance of two, three, and more defects within the spot.

3.1.2 Topological charges of the spot. To gain more insight into the positive, negative, and total topological charges within the activity spots, we first consider the time averages of the positive charge, 〈C₊〉 = +½ × 〈N₊〉, and the negative charge, 〈C₋〉 = −½ × 〈N₋〉, where 〈N₊〉 and 〈N₋〉 denote the time-averaged numbers of +½ and −½ defects, respectively. In Fig. 5(a), these quantities are plotted versus spot radius r. When the two +½ defects swirl around each other for r values in the region between the first and second vertical line (0.15 ≤ r ≤ 0.22), the total positive charge is +1. When the third positive defect starts to enter at r = 0.22, 〈C₊〉 jumps slightly and then increases steadily. The total negative charge becomes noticeable around r = 0.35 and then also increases steadily but it always falls below 〈C₊〉. Finally, for large spots (r ≳ 2.5), both 〈C₊〉 and 〈C₋〉 scale with r² (dashed line), which indicates that most of the defects are generated as +½ defect pairs within the spot indicating bulk active turbulence.

Fig. 5 (a) Time-averaged positive and negative topological charges inside an activity spot, 〈C₊〉 = +½〈N₊〉 (orange) and 〈C₋〉 = −½〈N₋〉 (blue), plotted versus spot radius r. The black dashed line indicates the scaling |〈C_±〉| ∝ r². (b) Total topological charge C_T plotted versus r using three different methods: (i) by counting defects, 〈C_defects〉 (red circles), (ii) by counting full directors turns along the rim, 〈C_turns〉 (blue stars), and (iii) by integrating the diffusive charge density over the spot area, 〈C_diffusive〉 (green squares). The black dashed line indicates the scaling ∝ r.

In Fig. 5(b) we plot the total topological charge versus spot radius r, obtained by different methods. We first concentrate on the direct defect counting (as before), and determine 〈C_defects〉 = ½〈N₊〉 −½〈N₋〉, [red circles in Fig. 5(b)]. The results are in line with the previous comments. Below r ≈ 0.15, the total charge rises from nearly zero with increasing radius and becomes one for the two swirling +½ defects. Then, beyond r ≈ 0.22, when the third +½ defect enters, total charge grows continuously with r. Interestingly, for large spots (r ≳ 2.5) the total charge scales linear in r. Thus, the density of the total defect charge goes to zero as expected in bulk turbulence behavior. Nevertheless, there is still some surplus of positive charge, which enters through the rim of the activity spot.

We verified the total-charge calculation for the active spot by two additional methods. First, we determined 〈C_turns〉 by counting the number of full director turns, when going around the rim of the spot.⁵³ Second, we determined 〈C_diffusive〉 by integrating the diffusive topological charge density⁵⁷ over the whole spot area,

setting S = 1 in the Q tensor. As Fig. 5(b) demonstrates, the total charges calculated from all three methods collapse onto each other.

3.2 What determines the topological defect charge?

We show in this section that the total topological charge in the activity spot consists of two contributions. First, although we do not impose any explicit surface anchoring for the director, the active flow along the periphery of the activity spot generates an effective active anchoring which, together with the circular geometry, introduces a baseline topological charge of +1, as we show in Section 3.2.2. Second, additional positive charges can enter through the rim of the activity spot, as we already mentioned. We discuss this further in the next section.

3.2.1 Characteristic length for defects to enter. We already saw in Fig. 5(b) that for large spots the total charge C_T grows linearly in r, which implies the existence of a characteristic length. To extract it, in Fig. 6 we plot the circumference of the spot, 2πr, divided by C_T − 1, where we subtract the charge +1 determined by geometry. We skip the data points for r ≤ 0.22, since the ratio diverges. Between the two dashed lines, (2 + 1) +½ defects exist in the activity spot. For small spots, the ratio is comparatively large and strongly size-dependent. But around r = 0.7 it reaches a plateau value of about 3.8, in our reduced units. This is the length to accommodate two additional +½ defects along the rim. The characteristic length of 1.9 for each additional +½ defect is comparable to the active length λ_a = 1.0 in our system. So the active motion provides the energy to generate the director distortions along the rim, necessary to introduce additional defects. The inset in Fig. 6 shows a snapshot of these director distortions for a spot with r = 1.05 and a momentary charge of +3.

Fig. 6 Ratio of circumference 2πr to C_T − 1 plotted versus spot radius r. The red line indicates the length required to accommodate an additional full charge +1 at sufficiently large r. The two vertical dashed lines enclose the region where (2 + 1) +½ defects are in the spot. Inset: Snapshot of director configuration at the rim for spot radius r = 1.05 with momentary total charge +3.

The director field along the rim is very dynamic, since defects constantly enter and leave the spot. This is illustrated in Fig. 7. For each time t, we start at the angular position θ = 0 on the rim (the vertical axis), walk along the rim up to the angular position θ, and determine how much the director is rotated relative to the local rim normal. The resulting number of this integrated director turn is given by Q(θ,t). Its values are color-coded in the θ–t plane in Fig. 7 for the spot radius r = 0.4. The director turn at θ = 0 is always zero by definition while after a full circle around the rim at θ = 2π, it can only assume integer and half-integer numbers. They correspond to the number of +½ defects, which have entered the spot, besides the +1 charge due to the active anchoring. Indeed, in Fig. 7 the additional charge Q(2π,t) assumes values of 0, ½, and 1 and strongly fluctuates in time.

Fig. 7 In the plane time t versus angular position θ on the rim, the accumulated director turn Q(θ,t) starting from θ = 0 is shown color coded for an activity spot of radius r = 0.4. Defects of charge +½ entering or leaving the spot are indicated by red or pink dots, respectively. All defects of charge −½, whether entering or leaving, are shown by light blue dots.

We can also identify horizontal lines where the color changes abruptly along the vertical time axis. Defects enter or leave at both ends of these lines. For example, at the left end a +½ defect enters when the color jumps to a brighter appearance, indicating an extra positive half turn, and it leaves when the jump is towards a darker color. The reverse is valid for the not actively moving −½ defects, where we do not distinguish between entering and leaving, and at the right end of the lines.

3.2.2 Active anchoring at the rim of the activity spot. Now we address the active anchoring at the rim of the spot, which is responsible for the +1 baseline topological charge. We remind of Section 2, where we showed using the linearized Q-tensor dynamics that strain-rate tensor E and Q tensor are proportional to each other, E ∝ Q. Thus, the nematic director aligns along the extensional axis of the local flow given by one of the principal axis of E. In the full nonlinear equations, one also has to take into account vorticity.

We first determine the extensional axis of the flow field and describe its orientation by the angle ψ_E, measured with respect to ê_θ, as illustrated in Fig. 8(b). Fig. 8(a) then shows probability distribution functions of ψ_E for three representative spot radii, r = 0.25, 0.4, and 2.5. They are determined from all angles ψ_E occurring along the rim and in time. All three distributions (solid lines) feature two pronounced peaks at ±45°. To understand these angles, we note that for an incompressible fluid, E_rr = −E_θθ, where we introduced the diagonal components of E in radial and tangential (azimuthal) directions, respectively. In the ideal case, where there is no variation along the tangential direction and neglecting the curvature of the rim, one has E_θθ = 0 = −E_rr. So, the extensional axis is oriented along ψ_E = ±45°. In the general case, one has

tan 2ψ_E = −E_rθ/E_rr. (6)

Fig. 8 The nematic director field and the extensional axis of flow at the rim of the activity spot. (a) Probability distributions of the director angle ψ (dashed) and the angle ψ_E of the extensional flow axis for three spot radii r = 0.25 (red), 0.40 (green), and 2.5 (blue). (b) Definitions of ψ and ψ_E relative to the tangential direction ê_θ. Their difference is the Leslie angle θ_d. (c) Probability distributions for E_rr (orange), E_rθ (blue), and the vorticity ω (green) at the rim of a spot with radius r = 0.40.

In Fig. 8(c), we plot the probability distributions for E_rr and E_rθ for one example radius, r = 0.40. Indeed, E_rr is narrowly concentrated around zero, while the distribution for E_rθ is broader. As a result, −E_rθ/E_rr often attains very large positive or negative values, which corresponds to angles ψ_E close to ±45°.

According to the linearized Q-tensor dynamics, we expect the local director with its orientation angle ψ [see Fig. 8(b)] to be oriented along the extensional flow axis, ψ = ψ_E = ±45°. However, in all three probability distributions of ψ in Fig. 8(a) (dashed lines) the maximum is shifted away from ±45° towards smaller absolute angles. To understand this shift, we need to consider the full nonlinear Q-tensor dynamics. From nematodynamics it is known that the tilt angle of the director away from the extensional axis of the flow field, the Leslie angle θ_d, is determined by⁵⁸

sin 2θ_d  =  ω/E, (7)

where ω is the vorticity and quantifies the strain rate. Eqn (7) implies that for |ω| < E the director tilts away from the extensional axis by an angle θ_d; in the limiting case ω = 0 it aligns exactly with that axis. Conversely, for |ω| > E the director undergoes continuous rotation (tumbling). As an illustrative example, we consider the case r = 0.40. Taking the locations of the peaks in the distribution of ω in Fig. 8(c) and combining them with the corresponding peaks in the distribution of E_rθ, which determines E, we obtain from eqn (7) the respective Leslie angles θ_d ≈ 17.25° and 14.65°. These angles nicely agree with the differences in the peak positions of P(ψ) and P(ψ_E), namely 17.87° and 14.65°.

We also realize that the maxima of the distributions P(ψ) are less pronounced compared to the distributions P(ψ_E). They flatten progressively as the radius increases. This can be understood by the fact that for defects entering or leaving the spot, the flow vorticity ω at the respective locations on the rim exceeds the strain rate E, which causes the director to tumble. Thus, the distribution P(ψ) widens, which becomes more pronounced at larger radii due to an increased number of defect traversals.

To conclude, the flow field with its extensional axis and vorticity generated by active motion at the rim causes the nematic director to align with a specific angle relative to the tangential direction. Therefore, this active anchoring causes the baseline +1 topological charge. However, defects entering and leaving the spot cause the director angle to deviate from this ideal value.

We stress that there is a clear difference between the active anchoring in our work and the active anchoring mechanism for extensile activity discussed in ref. 57 and 59. In general, shear flow develops at the nematic–isotropic interface whenever the director n is neither parallel nor perpendicular to the interface.⁵⁷ In the case of ref. 57 and 59, flow-aligning terms are explicitly omitted or small, so that only flow tumbling is considered, where shear flow can only rotate and not align the director. As a consequence, at the nematic–isotropic interface, shear rotates the director until planar alignment, where shear flow then switches off in the stationary state. In contrast, our system is dominated by flow alignment, which not only creates nematic order but also at the nematic–isotropic interface selects a steady Leslie angle rather than strict alignment parallel to the interface. Correspondingly, the interfacial shear flow does not switch off but persists in the stationary state. The anchoring angle is not universal but depends on local flow parameters, namely, the vorticity ω and the strain-rate strength E.

3.3 Measure of chaosity in fluid flow: Lyapunov exponent

Chaotic flow is one of the main characteristics of active turbulence. In this subsection, we measure the degree of chaos in the velocity fields generated within activity spots and see how it changes with the spot radius. Of course, defects are advected by the flow field and thereby reflect the degree of chaos with their motion. But note that +½ defects are also mobile.

To measure the degree of chaos, we determine the finite-time Lyapunov exponent (FTLE).⁶⁰ It quantifies how the distance between two tracer particles, dispersed in a flow field and initially at time t₀ separated by an infinitesimally small distance δ₀, evolves over time as δ₀e^λ(t−t₀), where λ is the Lyapunov exponent. Given the dynamic velocity field in our simulations, we are able to determine tracer trajectories as a function of time, x(t;t₀,x₀), starting from their initial position x₀ at initial time t₀. The flow map completely describes the tracer or particle trajectories, i.e., it flows each starting point x₀ at time t₀ forward to its new location at time t. The FTLE for the time interval t₀ to t is given by⁶⁰

where σ_max(x₀)^t₀→t is the largest eigenvalue of the right Cauchy–Green strain tensor, . Here, the gradient of the flow map, , also called Jacobian matrix, maps an initial small separation vector δx₀ to its appearance at time t, . Thus, [σ_max(x₀)^t₀→t]^1/2 determines the largest possible growth rate among all directions of δx₀.

For each activity spot we seeded a large ensemble of hypothetical tracer particles within the spot at positions x₀ and times t₀ and let them advect under the simulated velocity field until they exited the spot. This defines the end time t of a trajectory. For each trajectory we determined the largest eigenvalue σ_max(x₀)^t₀→t using the TBarrier package for part of the calculations.⁶¹ The eigenvalue σ_max(x₀)^t₀→t is plotted versus the elapsed time Δt = t − t₀ in the inset of Fig. 9). At every Δt, we then averaged the values of log σ_max over all trajectories that stayed in the spot for this time. This results in a single curve of mean values 〈log σ_max〉 (Δt), which we show in Fig. 9 together with the standard deviation for a spot radius r = 0.40. After an initial transient, the curve increases linearly in time Δt and twice the slope indicated by the fitted dashed line in Fig. 9 gives the averaged FTLE for one spot radius.

Fig. 9 Ensemble average 〈log σ_max〉 plotted versus elapsed time Δt = t − t₀ for spot radius r = 0.40. Blue dots show the mean value and the shaded band indicates one standard deviation. The black dashed line is a linear fit to the post-transient region, whose slope yields the FTLE . Inset: Log σ_max(Δt) plotted versus Δt for various example trajectories of the tracer particles with different lengths.

In Fig. 10 the values of the averaged FTLE are plotted as a function of the spot radius. The FTLE remains vanishingly small for spot radii below r = 0.24, where only two defects swirl around each other and the flow field is regular. Then, when a third defect enters the spot noticeably, the FTLE jumps to values between 0.03 and 0.04, indicating chaotic flow in agreement with the irregular motion of three defects. From there up to r ∼ 2.5, the FTLE highly fluctuates but gradually increases with the the spot radius. Ultimately, it reaches a plateau value at approximately 0.04 for all larger radii. This plateau value coincides with the FTLE for bulk turbulence. So again, once the spot exceeds r ≈ 2.5, bulk behavior is observed.

Fig. 10 Averaged FTLE as a function of spot radius. A sharp jump at r = 0.24 marks the onset of chaotic flow. The red dashed line indicates the bulk value of the FTLE.

3.4 Self-constraint in active flows

In a recent inspiring work on defect–flow coupling in ref. 62, the authors show that self-motile +½ defects are tightly constrained to viscometric surfaces, where the strengths of vorticity and strain rate balance each other. Thus, we analyze how close +½ defects move to the viscometric contours,and also to the zero-vorticity contours ω(r,t) = 0. Here, Ω and E are the respective vorticity and strain-rate tensors, introduced before, and ‖·‖ denotes the Frobenius norm, ‖A‖² = A_ijA_ij. The authors in ref. 62 reported that +½ defects preferentially localize near and move along viscometric contours, [scr Q, script letter Q] = 0, with a substantially stronger localization compared to zero-vorticity contours ω = 0 that are commonly used to delineate regions of positive and negative vorticity from each other.

Motivated by these findings, we perform the same analysis for our activity spots. For each snapshot, we determine the zero-vorticity contours, ω = 0, and the viscometric contours, [scr Q, script letter Q] = 0, and then measure the closest distance of the +½ defects from both type of contours. The distributions of distances, resulting from all snapshots, are shown in Fig. 11, where the red curves correspond to distributions of distances d_ω from zero-vorticity contours and the blue curves show the distributions of distances d_{[scr Q, script letter Q]} from viscometric contours. All distances are rescaled by the radius r of each spot. Our data show a consistently stronger localisation of +½ defects to contours of [scr Q, script letter Q] = 0 compared to ω = 0 across all spot radii considered. This supports the picture that +½ defects preferentially reside on and move along viscometric contours.

Fig. 11 Distributions of distances of +½ defects to the viscometric contours [scr Q, script letter Q] = 0 (blue) and zero-vorticity contours ω = 0 (red), for spot radii r = 0.25, 0.4, and 2.5; distances are rescaled by r.

The condition [scr Q, script letter Q] = 0 can occur in two ways. Generically, it corresponds to regions with linear shear flow, where the vorticity and strain-rate strengths balance, ‖Ω‖ = ‖E‖. It also includes the degenerate case where both vorticity and strain vanish. In two dimensions, the scalar vorticity ω fully specifies the antisymmetric part of ∇ ⊗ u, so points satisfying both [scr Q, script letter Q] = 0 and ω = 0 necessarily satisfy ‖E‖ = 0 and hence ∇ ⊗ u = 0, corresponding to locally uniform flow. In our simulations of the activity spot, Videos S8–S10 and representative snapshots (see Fig. 14 in Appendix C) indicate that +½ defects predominantly occupy and travel along [scr Q, script letter Q] = 0 contours, often close to ω = 0 contours. This means the defects move near locations, where the two types of zero-contour lines intersect or are close to each other. We capture this behavior quantitatively using the joint distribution function of defect distances, P(d_{[scr Q, script letter Q]},d_ω) (see Fig. 15 in Appendix C). An enhanced probability around d_{[scr Q, script letter Q]} = d_ω = 0 at small distances demonstrates that defects more often are simultaneously close to both types of contour lines.

4 Conclusions

In this article we investigated a circular activity spot that induces active paranematic order embedded in an otherwise passive isotropic fluid. Increasing the spot radius r enables a controlled route to active turbulence. To demonstrate this, we numerically solved the hydrodynamic equations of active paranematic flow for activity spots with increasing radii relying on our previous work.¹³ As the spot radius r grows from values comparable to the nematic coherence length λ_n, the director field in the spot first is uniform and then single transient +½ defects appear. Further increasing the radius, the spot robustly hosts two +½ defects that first form a spiral defect in the center; for larger spots they separate and circulate around each other with a nearly constant separation. Beyond this two-defect regime, additional +½ defects progressively enter the spot and at larger radii also −½ defects can enter from the rim. Finally, bulk active turbulence is achieved for spots with radii roughly 2.5 times the active length λ_a (see Fig. 2 and supplemental videos).

To quantify the topological defects within the spot, we determined the time-averaged positive and negative charges. They both scale as r² at large r so that their areal densities are constant. Furthermore, the total charge of the spot, the sum of the positive and negative charges, grows linearly with r because +½ defects are persistently entering the spot through its rim. So, the areal charge density tends to zero for large r, as expected for bulk behavior [see Fig. 5(b)], and the continuous creation and annihilation of pairs of +½ defects dominates. Finally, in addition to defect counting, we verified the total topological charge inside the spot by two further approaches: counting full director rotations along the rim and integrating the diffusive charge density over the spot area.

The total topological charge in the spot is determined by two contributions. First, the active flow at the rim of the spot aligns the director along the local extensional axis of the flow plus some additional tilt due to vorticity, which is quantified by the Leslie angle. This establishes active anchoring and a baseline topological charge of +1 in the absence of any explicit anchoring energy. Second, dynamic director distortions sustained by activity at the rim enable additional defects to continually enter and leave the spot. Time series confirm the corresponding changes in defect counts (see Fig. 3), and the space-time dynamics of the director distortions, beyond active anchoring, is illustrated in Fig. 7. For large spots we identified the arc length along the rim, which is necessary so that activity can induce the director distortions needed for a defect to enter through the rim. It is comparable to the active length λ_a.

Two dynamic regimes stand out. Once the spot hosts two +½ defects with either a spiral or circulating pair configuration, the enstrophy rises sharply indicating the appearance of strong vortices. The enstrophy increases further as more defects populate the spot. When a third +½ defect participates noticeably, the defect trajectories decorrelate and the flow becomes chaotic. This is signaled by a jump of the finite-time Lyapunov exponent starting from a nearly zero value. A similar behavior in a defect-free active-fluid model has also been reported recently.⁶³ Finally, at large radii with r > 2.5 both the enstrophy and the finite-time Lyapunov exponent approach their bulk plateaus.

Thus, our work identifies a controlled route from a quiescent paranematic spot to bulk active turbulence and provides a quantitative link between rim alignment, influx of defects, which cannot occur for solid boundaries, and flow complexity. For example, it should be realizable in bacterial active fluids, where the individual constituents have a typical length of ca. 10 µm. Assuming the nematic coherence length λ_n to be of the same order,^64,65 we estimate for the active length λ_a a value of 200 µm, so 20 times larger than λ_n as in our work. This length is comparable to values reported in literature^66–68 and it sets the vortex size and the typical distance of defects.

An immediate extension of our work is to study spots of different shapes and topologies, where in the latter case active anchoring sets different baseline charges. In all these cases different defect configurations and flow states are expected. For example, in an annular ring with respective inner and outer radii, r_in and r_out, active anchoring occurs at the inner and outer boundaries so that the baseline charge should be zero. We expect the total topological charge in the annular ring to scale with Δr = r_out − r_in. In particular, for Δr close to but larger than the nematic coherence length, we envision a director field in the ring that is tilted against the local normal along the radial direction, as dictated by active anchoring. For opposite tilt angles at the respective inner and outer boundaries, a uniform circling flow should occur, while for equal tilt angles counter-rotating swirls are possible. Increasing Δr, defects should enter along the ring and they might perform some regular motion along the azimuthal direction, including some possible braiding motion. Another geometry are arrays of spots, where active flow and defects could spill to neighboring spots, similar to the geometry and bacterial vortices in ref. 27 but now without solid boundaries.

Thus, our work paves the way to programmable and controllable flows in active fluids using currently available experimental means to spatially pattern activity in photosensitive materials. For example, genetically engineered photo-responsive bacteria allow the activity to be tuned by light,^38,39 and in microtubule-based active nematics activity patterning can be achieved either via fuel-based schemes using caged ATP⁴² or light-sensitive motor proteins.^40,41,43 In the latter case, motors are uniformly distributed but their binding state and hence activity is switched on and off by a prescribed light pattern, which thus produces an activity pattern fixed in the lab frame.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information: Movies S1–S10 showing director-field dynamics and ±½ defect motion for varying inactivity-spot radii, with the scalar order parameter (Movies S1–S7) or the velocity magnitude plus viscometric contours and zero-vorticity contours (Movies S8–S10). See DOI: https://doi.org/10.1039/d5sm01076d.

The data supporting the findings of this study are available upon request.

Appendices

A Nondimensional equations

We nondimensionalize our equations by introducing a characteristic length x₀ and time t₀ defined, respectively, aswhere D_r is the rotational diffusion coefficient of a single rod. Rescaling length and time in our dynamic equations, we arrive at the reduced velocity and the following dimensionless parameters:

In writing the nondimensional eqn (1) and (2), we drop the tilde for ease of notation and arrive at,

−∇·u = 0, (12)

B Velocity field in spots

In Fig. 12(a) we show the velocity field (streamlines) and vorticity (color map) inside a spot with r = 0.15, where we observe a spiral defect in the director field. A clear vortex flow in the center occurs surrounded by a region of weaker vorticity with opposite sign. The total vorticity in the spot is non-zero as shown in Fig. 13 since the flow spreads beyond the spot. The same is valid for the flow field in Fig. 12(b) for radius r = 0.20, where a pair of +½ defects circles around the center. However, now the central vortex is elongated and on both sides regions of strong opposite vorticity exist. Thus, both cases show a clear nonzero circulation in the center with a spontaneously selected swirl direction. Similar swirling flows have been reported in active nematic fluids, both under physical confinement^28,69 and in unconfined settings.⁷⁰

Fig. 12 Velocity fields (flow lines) and vorticity (color map) in activity spots with radii (a) r = 0.15 and (b) r = 0.20. The rim of the spot is indicated by a black dashed line.

Fig. 13 Absolute value of mean vorticity in the spot plotted versus spot radius r.

Fig. 13 shows that mean vorticity is strong in the region between the first and second vertical dashed lines, where only two defects circle around in the spot creating a circulating flow field. As soon as further defects enter, the mean vorticity drops sharply. With increasing radius positive and negative vorticities in the spot cancel each other and total vorticity stays close to zero.

C Topological-defect positions relative to flow contours

In Fig. 14, the +½ defects (red arrowed dots) and −½ defects (cyan dots) are overlaid with the viscometric contours [scr Q, script letter Q] = 0 (white curves) and the zero-vorticity contours ω = 0 (blue curves), on a background showing the velocity magnitude. The snapshots indicate that +½ defects most often lie close to both [scr Q, script letter Q] = 0 and ω = 0, that is near locations where the two contour families intersect or come close to each other. In two dimensions, [scr Q, script letter Q] = 0 together with ω = 0 implies ‖E‖ = 0 and hence ∇ ⊗ u = 0, so these intersections correspond to locally uniform flow. We quantify this by examining the joint distance statistics of d_{[scr Q, script letter Q]} and d_ω (Fig. 15). The strong concentration of probability near (d_{[scr Q, script letter Q]},d_ω) = (0,0) shows that +½ defects are simultaneously close to both contour sets, consistent with localisation near intersection segments of [scr Q, script letter Q] = 0 and ω = 0.

Fig. 14 Topological defects (red arrowed dots for +½ defects and cyan dots for −½ defects) overlaid with viscometric contours [scr Q, script letter Q] = 0 (white curves) and zero-vorticity contours ω = 0 (blue curves). The background color shows the velocity magnitude and black lines show the director field. Panels show activity spots with radii (a) r = 0.25, (b) r = 0.4, and (c) r = 2.5.

Fig. 15 Joint distribution P(d_{[scr Q, script letter Q]},d_ω) for distances d_{[scr Q, script letter Q]} and d_ω. Left: Full range of distances. Right: Zoom into the near-origin region indicated by the white box in the left panel.

Acknowledgements

We thank Tyler Shendruk, Julia Yeomans, Cristina Marchetti, Ricard Alert, and Markus Bär for helpful discussions. We acknowledge support from the Deutsche Forschungsgemeinschaft in the framework of the Collaborative Research Center SFB 910 and by a grant under STA352/11 as well as from TU Berlin.

Notes and references

1. 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.
2. A. Zöttl and H. Stark, J. Phys.:Condens. Matter, 2016, 28, 253001.
3. S. Ramaswamy, Annu. Rev. Condens. Matter Phys., 2010, 1, 323–345.
4. A. W. Zantop and H. Stark, Soft Matter, 2022, 18, 6179–6191.
5. M. J. Bowick, N. Fakhri, M. C. Marchetti and S. Ramaswamy, Phys. Rev. X, 2022, 12, 010501.
6. C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein and J. O. Kessler, Phys. Rev. Lett., 2004, 93, 098103.
7. H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 14308–14313.
8. T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431–434.
9. P. Guillamat, J. Ignés-Mullol, S. Shankar, M. C. Marchetti and F. Sagués, Phys. Rev. E, 2016, 94, 060602.
10. K. Suzuki, M. Miyazaki, J. Takagi, T. Itabashi and S. Ishiwata, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 2922–2927.
11. T. B. Saw, W. Xi, B. Ladoux and C. T. Lim, Adv. Mater., 2018, 30, 1802579.
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. A. Partovifard, J. Grawitter and H. Stark, Soft Matter, 2024, 20, 1800–1814.
14. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, Oxford, 1986.
15. R. Aditi Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.
16. S. P. Thampi, A. Doostmohammadi, R. Golestanian and J. M. Yeomans, Europhys. Lett., 2015, 112, 28004.
17. S. Santhosh, M. R. Nejad, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, J. Stat. Phys., 2020, 180, 699–709.
18. D. Dell’Arciprete, M. L. Blow, A. T. Brown, F. D. Farrell, J. S. Lintuvuori, A. F. McVey, D. Marenduzzo and W. C. Poon, Nat. Commun., 2018, 9, 4190.
19. J. Hardoüin, 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.
20. C. Joshi, S. Ray, L. M. Lemma, M. Varghese, G. Sharp, Z. Dogic, A. Baskaran and M. F. Hagan, Phys. Rev. Lett., 2022, 129, 258001.
21. S. Edwards and J. Yeomans, Europhys. Lett., 2009, 85, 18008.
22. L. Giomi, L. Mahadevan, B. Chakraborty and M. Hagan, Phys. Rev. Lett., 2011, 106, 218101.
23. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246.
24. R. Di Leonardo, L. Angelani, D. Dell’Arciprete, G. Ruocco, V. Iebba, S. Schippa, M. P. Conte, F. Mecarini, F. De Angelis and E. Di Fabrizio, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 9541–9545.
25. S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti and V. Vitelli, Nat. Rev. Phys., 2022, 4, 380–398.
26. H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 268102.
27. H. Wioland, F. G. Woodhouse, J. Dunkel and R. E. Goldstein, Nat. Phys., 2016, 12, 341–345.
28. A. Opathalage, M. M. Norton, M. P. Juniper, B. Langeslay, S. A. Aghvami, S. Fraden and Z. Dogic, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 4788–4797.
29. F. L. Memarian, D. Hammar, M. M. H. Sabbir, M. Elias, K. A. Mitchell and L. S. Hirst, Phys. Rev. Lett., 2024, 132, 228301.
30. P. Guillamat, J. Ignés-Mullol and F. Sagués, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 5498–5502.
31. K. Thijssen, D. A. Khaladj, S. A. Aghvami, M. A. Gharbi, S. Fraden, J. M. Yeomans, L. S. Hirst and T. N. Shendruk, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2106038118.
32. M. M. Norton, P. Grover, M. F. Hagan and S. Fraden, Phys. Rev. Lett., 2020, 125, 178005.
33. S. Ghosh, A. Baskaran and M. F. Hagan, J. Chem. Phys., 2025, 162, 134902.
34. C. Floyd, A. R. Dinner and S. Vaikuntanathan, Soft Matter, 2025, 21, 4488–4497.
35. J. Rønning, M. C. Marchetti and L. Angheluta, R. Soc. Open Sci., 2023, 10, 221229.
36. S. Shankar, L. V. Scharrer, M. J. Bowick and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2400933121.
37. M. M. Norton, A. Baskaran, A. Opathalage, B. Langeslay, S. Fraden, A. Baskaran and M. F. Hagan, Phys. Rev. E, 2018, 97, 012702.
38. J. Arlt, V. A. Martinez, A. Dawson, T. Pilizota and W. C. Poon, Nat. Commun., 2018, 9, 768.
39. G. Frangipane, D. Dell'Arciprete, S. Petracchini, C. Maggi, F. Saglimbeni, S. Bianchi, G. Vizsnyiczai, M. L. Bernardini and R. Di Leonardo, eLife, 2018, 7, e36608.
40. R. Zhang, S. A. Redford, P. V. Ruijgrok, N. Kumar, A. Mozaffari, S. Zemsky, A. R. Dinner, V. Vitelli, Z. Bryant and M. L. Gardel, et al., Nat. Mater., 2021, 20, 875–882.
41. L. M. Lemma, M. Varghese, T. D. Ross, M. Thomson, A. Baskaran and Z. Dogic, PNAS Nexus, 2023, 2, pgad130.
42. I. Vélez-Cerón, J. Ignés-Mullol and F. Sagués, Soft Matter, 2024, 20, 9578–9585.
43. Z. Zarei, J. Berezney, A. Hensley, L. M. Lemma, N. Senbil, Z. Dogic and S. Fraden, Soft Matter, 2023, 19, 6691–6699.
44. K. Chaithanya, A. Ardaševa, O. J. Meacock, W. M. Durham, S. P. Thampi and A. Doostmohammadi, Commun. Phys., 2024, 7, 302.
45. A. Mozaffari, R. Zhang, N. Atzin and J. J. de Pablo, Phys. Rev. Lett., 2021, 126, 227801.
46. G. B. Jeffery, Proc. R. Soc. London, Ser. A, 1922, 102, 161–179.
47. A. N. Beris and B. J. Edwards, Thermodynamics of Flowing Systems: With Internal Microstructure, Oxford University Press, Oxford, 1994.
48. H. Isliker, A Tutorial on the Pseudo-Spectral Method, 2004, https://www.astro.auth.gr/vlahos/GravitoplasmaWS1/pseudo-spectral_2.pdf, University of Thessaloniki, accessed September 10 2025.
49. L. N. Trefethen, Spectral Methods in MATLAB, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2000, vol. 10, p. 181.
50. M. Hochbruck and A. Ostermann, SIAM J. Numer. Anal., 2005, 43, 1069–1090.
51. A. J. Chorin and J. E. Marsden, A Mathematical Introduction to Fluid Mechanics, Springer, New York, 2nd edn, 1990, vol. 4, p. 168.
52. A. U. Oza and J. Dunkel, New J. Phys., 2016, 18, 093006.
53. P. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, 1993.
54. L. Angheluta, Z. Chen, M. C. Marchetti and M. J. Bowick, New J. Phys., 2021, 23, 033009.
55. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. C. Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365.
56. S. Shankar, S. Ramaswamy, M. C. Marchetti and M. J. Bowick, Phys. Rev. Lett., 2018, 121, 108002.
57. M. L. Blow, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2014, 113, 248303.
58. J. Ericksen, Kolloid-Z., 1960, 173, 117–122.
59. M. L. Blow, M. Aqil, B. Liebchen and D. Marenduzzo, Soft Matter, 2017, 13, 6137–6144.
60. G. Haller, Phys. D, 2001, 149, 248–277.
61. A. P. Encinas Bartos, B. Kaszás and G. Haller, TBarrier, 2023 DOI:10.5281/zenodo.6779400.
62. L. C. Head, C. Doré, R. R. Keogh, L. Bonn, G. Negro, D. Marenduzzo, A. Doostmohammadi, K. Thijssen, T. López-León and T. N. Shendruk, Nat. Phys., 2024, 20, 492–500.
63. M. Hillebrand and R. Alert, Nat. Commun., 2025, 16, 11169.
64. R. Zhang, Y. Zhou, M. Rahimi and J. J. de Pablo, Nat. Commun., 2016, 7, 13483.
65. J. M. Barakat, K. J. Modica, L. Lu, S. Anujarerat, K. H. Choi and S. C. Takatori, ACS Appl. Nano Mater., 2024, 7, 12142–12152.
66. J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 228102.
67. K.-T. Wu, J. B. Hishamunda, D. T. N. Chen, S. J. DeCamp, Y.-W. Chang, A. Fernández-Nieves, S. Fraden and Z. Dogic, Science, 2017, 355, eaal1979.
68. H. Baza, F. Chen, T. Turiv, S. V. Shiyanovskii and O. D. Lavrentovich, Soft Matter, 2026, 22, 1297–1313.
69. F. G. Woodhouse and R. E. Goldstein, Phys. Rev. Lett., 2012, 109, 168105.
70. M. R. Nejad and J. M. Yeomans, Phys. Rev. X, 2023, 1, 023008.

---