Defect turbulence in a dense suspension of polar, active swimmers

We study the effects of inertia in dense suspensions of polar swimmers. The hydrodynamic velocity field and the polar order parameter field describe the dynamics of the suspension. We show that a dimensionless parameter $R$ (ratio of the swimmer self-advection speed to the active stress invasion speed) controls the stability of an ordered swimmer suspension. For $R$ smaller than a threshold $R_1$, perturbations grow at a rate proportional to their wave number $q$. Beyond $R_1$, we show that the growth rate is $\mathcal{O}(q^2)$ until a second threshold $R=R_2$ is reached. The suspension is stable for $R>R_2$. We perform direct numerical simulations to investigate the steady state properties and observe defect turbulence for $R<R_2$. An investigation of the spatial organisation of defects unravels a hidden transition: for small $R\approx 0$ defects are uniformly distributed and cluster as $R\to R_1$. Beyond $R_1$, clustering saturates and defects are arranged in nearly string-like structures.


I. INTRODUCTION
The instability [1,2] of a viscous suspension of motile organisms leads to an abundance of ordered and chaotic states [3][4][5][6][7][8][9][10], of which active turbulence [11] has attracted much attention.The majority of studies have investigated this iconic phenomenon in nematic systems, in which the local structure and characteristic defects are associated with uniaxial fore-aft symmetric order.In a recent work [1], we established the stabilizing role of fluid inertia in the transition to active turbulence in polar active suspensions, in which the local order is vectorial.We discovered that the steady-state behavior of the suspension changed across two thresholds R 1 and R 2 in the dimensionless ratio R ≡ ρv 2 0 /σ of inertia and activity, where ρ is the mass density of the suspension, v 0 is the self-propulsion speed, and σ the mean active stress.We found defect turbulence for R < R 1 , a fluctuating but ordered phase-turbulent steady state for R 1 < R < R 2 , and a quiescent ordered state for R > R 2 .The analysis of [1] considered only the suspension velocity and polar order-parameter fields u and p, implicitly treating the concentration of active particles as "fast", as would happen, for example, if birth and death led to a nonzero relaxation rate at zero wavenumber.
Highly concentrated systems present another scenario where we can ignore fluctuations in the number den- * perlekar@tifrh.res.insity of suspended particles, albeit in a radically different manner.In such a limit, it is reasonable to make the approximation that concentration fluctuations are prohibited, which means that the solute velocity field or, equivalently, the polar order parameter p, is solenoidal: ∇ • p = 0. We will term such systems as simply "dense".See [12] for the implementation of such a condition in dry active matter.
In this paper, we study stability and turbulence in a dense suspension of polar active particles (DSPAP).Section II describes the equations of motion and in Section III, we conduct a linear stability analysis of an ordered suspension with respect to perturbations with wavenumber q.We show that for extensile suspensions, the same dimensionless parameter R as in [1] governs the linear stability of the suspension.For R < R 1 ≡ 1 + λ, where λ is the flow-alignment parameter, we find an inviscid instability where the most unstable pure bend modes grow at a rate ∝ q.For R 1 < R < R 2 , we find that the pure bend perturbations grow at a rate ∝ q 2 .Finally for R > R 2 , orientational order is linearly stable.These instability mechanisms are identical to those in our earlier report [1].Further, dense suspensions of contractile swimmers are linearly stable.Next, in Section IV, we study the steady-state properties of the suspension using numerical simulations.We show the presence of vortexdefect turbulence for R < R 2 .Correlation length grows as we increase R and appears to diverge at R ≈ R 2 .Interestingly we do not observe a phase-turbulent regime [1].Instead, a detailed analysis of the spatial organization of defects using correlation dimension d 2 reveals a novel morphological transition.We find that defects are uniformly distributed (d 2 ≈ 2) for R → 0 (Stokesian suspension) and start to cluster on increasing R. Maximum clustering (d 2 ≈ 1) is attained around R = R 1 , and we observe that defects organize into stringy patterns.No further changes in defect organization are observed on increasing R beyond R 1 .

II. EQUATIONS OF MOTION
We use a hydrodynamic formulation to study the dense suspension of polar active particles (DSPAP).In the incompressible limit (uniform suspension density and active particle concentration), the dynamics of the hydrodynamic velocity field u(x, t) and the orientation order parameter p(x, t) is described by the following equations [1][2][3][13][14][15]: Here, ρ is the suspension density, µ is the fluid viscosity, v 0 is the self-advection speed of the swimmers, λ is the flow alignment parameter, P is the hydrodynamic pressure that enforces incompressibility of the velocity field, Π is the pressure-like term that enforces incompressibility of the order parameter field p imposed by the constant concentration approximation, S ≡ (∇u + ∇u T )/2 and Ω ≡ (∇u − ∇u T )/2 are the symmetric and antisymmetric parts of the velocity gradient tensor ∇u.Σ a ≡ −σ 0 pp is the leading order apolar intrinsic stress associated with the swimming activity where the forcedipole density σ 0 > 0 (< 0) for extensile (contractile) swimmers [2], Σ r ≡ −λ + hp − λ − ph is the reversible thermodynamic stress [15] with λ ± = (λ ± 1)/2, h = −δF /δp is the molecular field conjugate to p with the free-energy functional that favors an aligned order parameter state with unit magnitude, and Γ is the rotational mobility for the relaxation of the order parameter field to the uniform ordered state prescribed by the free energy dynamics.For simplicity, we choose a single Frank constant K, which penalizes the gradients in p [16].We are primarily interested in the interplay of the self-propulsion speed v 0 and the leading order apolar active stress, hence similar to our earlier study [1], we have ignored the contribution from higher-order polar gradient terms in Σ a .

III. LINEAR STABILITY ANALYSIS
We analyse the stability of the ordered state (u, p) = (0, x) to small monochromatic perturbations of the form (δu, δp) = ( u, p) e i(q•x−ωt) , where q is the perturbation wavevector and ω is the frequency.As dense contractile suspension are linearly stable, here, we focus our study on extensile suspension.In what follows, we discuss the results for pure-bend perturbations, the most unstable modes for extensile systems [1,4].For a detailed discussion of the linear stability analysis including the stability of twist-bend and splay-bend modes, we refer the reader to Appendix A. A large wavelength (small q) expansion up to O(q 2 ) yields the following dispersion relation for the pure-bend modes: where we have defined the dimensionless numbers R ≡ ρv 2 0 /2σ 0 , R 1 ≡ 1 + λ and β = ΓKρ/µ.In Fig. 1, we show a qualitative phase diagram highlighting the three different stability regimes.In regime A, R < R 1 and irrespective of the value of β, pure-bend modes are unstable with a growth rate 2 /4β, pure-bend modes grow at a rate I(ω) ∝ q 2 .Finally, in regime C, when R > R 2 the ordered state is stable.

IV. DIRECT NUMERICAL SIMULATIONS
We perform direct numerical simulations of (1) on a square domain of area L 2 discretized with N 2 equispaced collocation points.Similar to [1], we use a pseudo-spectral method for spatial integration of the velocity field and a fourth-order finite difference method For R ≥ 0.9, we only show a subdomain of size (40π) 2 .Inter-defect separation grows with increasing R.For R < R1 = 1.1, we find that the defects are uniformly distributed over the simulation domain, but form clusters for R > R1 (See also Fig. 5).
for the order parameter.For time marching we use a second-order Adams-Bashforth scheme [17], and the incompressibility condition in the order-parameter field is implemented using an operator-splitting method [18].We undertake high-resolution numerical studies at various values of R and fixed β = 10 −2 to characterize the turbulent states arising from the instabilities of the ordered state.We initialize our simulations with a perturbed ordered state u(x, 0) = 0+A 10 i=1 cos (q i x) y and p(x, 0) = x + B 10 i=1 cos (q i x) y, where q i = i(2π/L), and choose A = B = 10 −3 .For R < R 2 , the perturbations destabilise the flow and a defect turbulence state is achieved.Table I summarizes all our simulation parameters.In what follows, we discuss the statistical properties of defect turbulence with varying R.

A. Defect turbulence
In Fig. 2 we show the streamlines of the order parameter field p in the statistically steady state at different values of R. In both regimes A and B, the order parameter field shows vortices and saddles, with no evidence of global polar order.Asters and spirals are ruled out by (A,B) Scatter plot of the vortex cores for R = 0.3 and R = 4.0 respectively.For R = 0.3 and R ≪ R1 in general, defects are uniformly distributed over the simulation domain.For R ≳ R1, we observe clustering of defects on one-dimensional string like structures.(C) Plot of the cumulative radial distribution function p(r) for different values of R. At smaller distances, we observe p(r) ∼ r d 2 , with d2 ∼ 2 for R ≪ R1 and d2 ∼ 1 for R ≥ R1 p(r) ∼ r 1 .At distances larger than typical cluster size, p(r) ∼ r 2 shows the same scaling for all R. Note that we have rescaled each p(r) curve by its maximum value to highlight different scaling regimes.Inset: Variation of the correlation dimension d2 vs R obtained from least-square fit on p(r) for small r.
the incompressibility constraint.To further verify that the turbulent states lack global order, we compute the magnitude of the average polar order parameter |⟨p⟩| in the statistical steady state at different values of R, where ⟨. ..⟩ denotes spatio-temporal averaging.Note that for the disordered states, |⟨p⟩| = 0, whereas |⟨p⟩| = 1 for a perfectly aligned state.In Fig. 3 we plot |⟨p⟩| with increasing R, and, as expected, it is close to zero in both the unstable regimes.

B. Defect correlations and clustering
In Fig. 4(A), we plot the isotropic correlation function where ⟨. ..⟩ denote ensemble and angle averaging, for different values of R. With increasing R, the correlations of the order parameter increase.We fit the functional form δ at small r ≪ L to extract the correlation length ξ and the exponent δ which decreases monotonically from δ ≈ 1.7 for R = 0.2 to δ → 1 for large R. The correlation function C(r) collapses onto a unique curve when plotted versus the scaled separation r/ξ (see Fig. 4(A,inset)).In Fig. 4(B), we show that the correlation length ξ increases with R, and from the intercept of the linear fit on the 1/R axis, it appears to diverge at R = R 2 .Note that investigating the spatial structure of the order parameter field for R → R 2 becomes numerically unfeasible as ξ → ∞ and finite-size effects become important.
To quantify the spatial distribution of the defect cores, we begin by identifying the defect coordinates.In Fig. 4(B,inset), we plot the average nearest-neighbor separation d min for different values of R. For R < R 1 , the correlation length ξ and d min are indistinguishable, indicating that a unique length scale describes the defect dynamics [19].In contrast, for R > R 1 , correlation length increases, whereas the average nearest neighbour separation d min saturates.Consequently, in Fig. 4(C), the defect number density N (ξ) (number of defects per unit area), also scales differently with the correlation length ξ for R < R 1 and R > R 1 .For R < R 1 , where a single length scale governs the dynamics, we observe N (ξ) ∼ ξ −2 , which indicates that the defects are distributed uniformly over the entire domain [20][21][22].In contrast, for R > R 1 , we find that N (ξ) ∼ ξ −1 , indicating clustering of defects.
Indeed, the scatter plots of the defect coordinates in Fig. 5(A,B) indicate uniformly distributed defects for R < R 1 , whereas they appear clustered for R > R 1 .We further quantify the clustering by evaluating the correlation dimension d 2 from the defect positions.The correlation dimension d 2 is evaluated from the r → 0 scaling behaviour of the probability p(r) of finding two defects within a distance r [23,24].For R ≪ 1, we find d 2 ∼ 2, and it decreases with increasing R until it saturates around d 2 ∼ 1 for R > R 1 (see Fig. 5(C)).Thus we conclude that in dense suspensions, the cross-over in the O(q) to O(q 2 ) instability around R = R 1 is marked by an intriguing defect clustering transition.
In order to gain further insight on the defect clustering transition, we compute the Shannon entropy density of the polar order parameter H p and the defect arrangement H D using a two-dimensional extension of the pattern matching method from information theory [25][26][27].To find H p , we apply the pattern matching method on the discretized orientation field [28].For the entropy density H D , we apply the pattern matching algorithm on a boolean field which is set to one at the defect locations and zero everywhere else (See Fig. 5).In Fig. 6, we plot H p and H D versus R. As the defect position field contains less information in comparison to the order parameter field, H D is smaller than H p for all R.Both entropy densities decrease as we increase R, and scale inversely with the correlation length roughly as H D , H p ≈ ξ −1. 4 .However, we could not capture a pronounced change in the trend of H p or H D that corresponds to the clustering transition around R ≃ R 1 .

C. Energy Spectrum
The shell-averaged energy spectra for the velocity and the order parameter field is defined as , and where u m and p m are the Fourier coefficients of the velocity u and the order parameter p fields.In Fig. 7, we plot E u (qξ) and E p (qξ) for different values of R. Consistent with our correlation function plots, we find that the spectra collapses onto single curve for q < q σ , where q σ = 2π/ℓ σ with ℓ σ ≡ µ/ √ ρσ 0 .Order parameter spectra shows two distinct power law scaling regimes: where ℓ σ = µ/ √ ρσ 0 , and L determines the smallest wavenumber available in the simulation domain.The q −3 scaling is consistent with Porod's law prediction [29].The kinetic energy spectrum E u (qξ) is determined by the balance of the viscous and active stresses where f q = P • iq • pp q , P is the projection operator, and ⟨⟨•⟩⟩ denotes temporal averaging.
In Fig. 7 we show an excellent agreement between the energy spectrum obtained directly from the velocity field and using (6), confirming the dominant balance between viscous and active stresses.Further assuming f q to be Gaussian random variables and using Gaussian integration by parts [30] we get (7) We observe that the prediction (7) matches well with the energy spectrum for 1 < qξ < q σ ξ (see Fig. 7).

V. CONCLUSIONS
We study spatio-temporal properties of dense wet suspensions of polar active particles.Using a linear stability analysis, we show how inertia can stabilize the orientational order against small perturbations.Our study reveals that the dimensionless parameter R ≡ ρv 2 0 /2σ 0 characterizes the stability of the aligned state.Although the neutral stability curve for DSPAP is identical to the previously studied case where concentration fluctuations are rendered fast [1], our numerical studies reveal that their steady-state properties are very different.For R < R 2 we observe defect turbulence, the order parameter flow consists of topological defects (vortices and saddles) with no global polar or nematic order.We unravel a hidden defect-ordering transition by investigating the spatial organization of defect centers.For R = 0, defects are uniformly distributed and start to cluster with increasing R. The clustering saturates around R = R 1 where we observe that the defects organize onto nearly one-dimensional string-like structures.Finally, we show that the spectrum of the order-parameter field shows a Porod's scaling for qξ >> 1 and a balance of viscous and apolar active stress determines the kinetic energy spectrum of the suspension velocity.S. R. acknowledges research support from a J. C. Bose Fellowship of the SERB, India.be verified directly from (A3) and (A4).For example, the pure splay modes (ϕ = π/2) relax with rates ω s = −i µ ρ q 2 , −2iΓ , and twist-bend modes are always stable as R < 0.

FIG. 1 .
FIG. 1. R − β phase diagram highlighting different stability regimes for incompressible polar active suspensions.In both the unstable regimes A and B, we observe defect turbulence (see Section IV).Red squares mark the simulations on the R − β plane with β = 10 −2 , R1 = 1.1, and R2 = 28.See TableIfor the rest of parameters.

FIG. 2 .
FIG.2.Order parameter streamlines superimposed over pseudocolor plot of ∇ × p highlighting vortices at different values of R. For R ≥ 0.9, we only show a subdomain of size (40π) 2 .Inter-defect separation grows with increasing R.For R < R1 = 1.1, we find that the defects are uniformly distributed over the simulation domain, but form clusters for R > R1 (See also Fig.5).

2 FIG. 3 .
FIG. 3. Average order | ⟨p⟩ | in the steady state at different R. For a disordered state | ⟨p⟩ | = 0 and for perfect alignment, | ⟨p⟩ | = 1.We do not observe any order in regime A or regime B. Ordered state is stable to perturbations in the green shaded region (R > R2).The error bars are smaller than the markers.
FIG. 4. (A) Steady state correlation function C(r) for different values of R. (Inset) Correlation function collapse onto a single curve when distance is scaled by the correlation length ξ. (B) Plot of the inverse correlation length 1/ξ versus 1/R.From the intercept of the linear fit on the 1/R axis, we conclude that ξ diverges around R ≈ R2.Inset: Comparison between the correlation length ξ and average nearest-neighbour distance dmin.For small R, ξ and dmin are identical.For large R, dmin saturates and ξ diverges.(C) Defect density N (ξ) as a function of the correlation length ξ.We observe two distinct scaling regimes for R < R1(> R1).Vertical black line marks the correlation length ξ(R = R1) computed from the linear fit in (B).

4 FIG. 6 .
FIG. 6.Plot of the Shannon entropy density Hp of the polar order parameter and the Shannon entropy density of the defect arrangement HD.With increasing R, both Hp and HD decrease and scale roughly as ξ −1.4 (dashed black guiding line).

4 FIG. 7 .
FIG. 7. (A)The order parameter energy spectrum Ep(qξ) for different values of R in both regime A and regime B. We scale the spectra by their respective peak values.Ep(qξ) shows a Porod scaling[29] for 1 < qξ < qσξ.(B) Comparison of the kinetic energy spectrum Eu(qξ) (solid markers) with (6) (solid lines) and (7) (dashed lines) for two different values of R.