# Multi-scale statistics of turbulence motorized by active matter

## Abstract

A number of micro-scale biological flows are characterized by spatio-temporal chaos. These include dense suspensions of swimming bacteria, microtubule bundles driven by motor proteins, and dividing and migrating confluent layers of cells. A characteristic common to all of these systems is that they are laden with active matter, which transforms free energy in the fluid into kinetic energy. Because of collective effects, the active matter induces multi-scale flow motions that bear strong visual resemblance to turbulence. In this study, multi-scale statistical tools are employed to analyze direct numerical simulations (DNS) of periodic two- (2D) and three-dimensional (3D) active flows and compare them to classic turbulent flows. Statistical descriptions of the flows and their variations with activity levels are provided in physical and spectral spaces. A scale-dependent intermittency analysis is performed using wavelets. The results demonstrate fundamental differences between active and high-Reynolds number turbulence; for instance, the intermittency is smaller and less energetic in active flows, and the work of the active stress is spectrally exerted near the integral scales and dissipated mostly locally by viscosity, with convection playing a minor role in momentum transport across scales.

## 1 Introduction

The multi-scale processes observed in the types of flows discussed here are induced by active matter laden in a fluid (Wensink et al., 2012; Sanchez et al., 2012; Dunkel et al., 2013). These are a special class of multi-phase flows, where the constituent particles are self-propelled. Examples of biological active matter are cells, motor proteins and bacteria. Synthetic active matter can be manufactured in the form of mechanically, chemically or optically propelled particles. However, a unifying characteristic of laden active matter is that it transforms free energy in the fluid into systematic motion (Simha & Ramaswamy, 2002). Although such energy conversion occurs at the particle scales, the collective interactions among many of these particles oftentimes translate into unstable flow motion across much larger scales.

Despite the low Reynolds numbers involved, flows induced by active matter have been referred to as active turbulence in analogy to the unsteady multi-scale dynamics found in high-Reynolds number flows (Wensink et al., 2012; Bratanov et al., 2015; Giomi, 2015). However, whether these flows can be identified as turbulent in the classical sense is debatable. The conservation equations of active flows are considerably non-linear and formally more complex than the incompressible Navier-Stokes equations, and the nature of these non-linearities is fundamentally different from those responsible for high-Reynolds number turbulence. Non-linear dynamical systems are known to display chaotic particle trajectories and mixing behavior even at low Reynolds numbers (Ottino, 1990).

In this study multi-scale tools are employed to examine basic flow statistics in 2D and 3D. The remainder of this paper is organized as follows. The formulation and computational set-up are described in §2. The analysis of the numerical results is presented in §3. Finally, concluding remarks are provided in §4.

## 2 Formulation and computational setup

The approach employed here is based on the numerical integration of the conservation equations of dense active nematohydrodynamics. This continuum formulation extends the description of passive liquid crystals of De Gennes & Prost (1995) and has proven successful in reproducing spatio-temporal dynamics observed in experiments of active flows (see Doostmohammadi et al. (2016b) and references therein).

### 2.1 Conservation equations for active nematohydrodynamics

In this formulation, the mesoscopic orientational order of active particles is represented by the nematic tensor , where is the magnitude of the nematic order, is the director and is the Kronecker delta. The conservation equation for the nematic tensor is given by

(1) |

where are velocity components, is time, are spatial coordinates, and is proportional to a rotational diffusivity. Additionally, is a standard co-rotation tensor defined as

(2) |

which accounts for the response of the orientational field to the extensional and rotational components of the velocity gradients, with and being the strain-rate and vorticity tensors, respectively. The relative importance of the vorticity and strain rate is controlled by the parameter , which characterizes the alignment of the nematics with the flow. The term involving in (1) describes, through the auxiliary tensor , the relaxation of to a minimum of a free energy

(3) |

In this formulation, represents the variational derivative, denotes the trace, is an elastic constant, and , and are material constants that determine the equilibrium state of the orientational order. The first three terms in (3) correspond to the Landau/De Gennes free energy, while the last term represents the Frank elastic energy with a one-constant approximation (De Gennes & Prost, 1995). Equation (3) is the free energy expansion in terms of the order parameter, where the terms allowed by symmetry have been retained. Note that the Frank term in (3) gives rise to a diffusive flux of in the form in (1). The description of the flow field is completed by the mass and momentum conservation equations, namely

(4) |

where is the density, is the pressure and is the dynamic viscosity. The additional terms in the momentum equation (4) involve the elastic stress

(5) |

which represents the passive conformational resistance of the nematics to changes in the orientational order (Edwards et al., 1991), and the active stress , which corresponds to a coarse-grained collective effect of the stresslets set up by the active particles (Simha & Ramaswamy, 2002). In particular, the divergence of the active stress is responsible for the injection of kinetic energy through the particle alignment characteristics represented by the nematic tensor, whose conservation equation (1) is non-linearly coupled with the velocity field. The absolute value of the coefficient is proportional to the activity, where positive and negative values of correspond, respectively, to extensile and contractile particles, the former being the case addressed in the computations shown below.

### 2.2 Remarks on the conservation equations

Some aspects with regard to the formulation given above are worth briefly pointing out. The first one is related to characteristic dimensionless parameters. In conditions addressed by these simulations (i.e., see further details in §2.3) and most experimental observations, the Reynolds number is small, , with and as the integral scales for velocity and length, respectively. In this limit, the convective transport of momentum has a diminishing influence, with the dominant balance in (4) being between viscous and active stresses. Similarly, the orientational order diffuses slower than momentum, as indicated by the relatively high Schmidt number , which results in characteristic structures of the velocity field that are larger than those of the orientational order, as discussed in §3. As a consequence, the Péclet number of the nematic order suggests that advection may play a more important role in transporting orientational order than it does for transporting momentum.

The second aspect worth stressing is that the conservation equations outlined above do not include any obvious external forcing term aimed at sustaining the dynamics. A flow laden with active matter is different from an inactive flow externally driven by boundary conditions, imposed shear or volumetric forces added to the momentum conservation equation. In real biological flows laden with swimming bacteria, adenosine tri-phosphate (ATP) molecules are consumed by the bacteria via hydrolysis reactions and transformed into motion in such a way that the driving occurs internally at the expense of depletion of chemical energy. The source of chemical energy, however, is not represented in the formulation above, in that it only concerns systems where the rate of ATP depletion is infinitesimally small compared to hydrodynamic processes.

The active motorization of the flow can be understood by multiplying the momentum equation in (4) by and performing a surface (in 2D) or volumetric (in 3D) periodic integration, which leads to the balance equation for the spatially averaged kinetic energy , where is the viscous dissipation. The work done by the active forces, given by the term , represents the main source of . The power deployed by the active work is dissipated by viscosity as shown below, thereby yielding a stationary state in which remains mostly constant.

### 2.3 Computational set-up

The formulation described above is integrated numerically in quasi-2D and 3D domains with periodic boundary conditions. In the quasi-2D framework the velocity varies in 2D while the nematic directors could develop out of plane components, as in previous studies (Thampi et al., 2013; Doostmohammadi et al., 2016b; Saw et al., 2017). The dimensional parameters in the simulations are , for 2D and for 3D, , (tumbling particles), , , , and (). The baseline activity coefficient is for 2D and 3D, with an additional 2D computation being performed with a smaller activity . All parameters here are expressed in simulation units. This generic set of parameters has been shown in earlier work to reproduce qualitatively the active turbulent state observed in bacteria and microtubule/motor-protein suspensions (Wensink et al., 2012; Sanchez et al., 2012; Thampi et al., 2013). A standard lattice-Boltzmann approach is used to integrate (4), while (1) is solved by employing a second-order finite-difference predictor-corrector algorithm. The resulting set of ordinary differential equations is integrated in time using an Euler method. The number of grid points for the 2D and 3D simulations is and , respectively, with a minimum grid spacing of for all cases, with being the characteristic size of the topological defects in the orientational order field.
The time increment used in the numerical integrations is , with being a characteristic time scale for the damping of the activity by viscosity. The initial conditions consist of zero velocity everywhere while the directors are set to random orientations. Data are collected after approximately time steps for 10 snapshots spanning a time period of , which provide and statistical samples leading to converged probability density functions (PDFs) in 2D and 3D, respectively. In the results, spatial coordinates are normalized with . Additionally, velocities are normalized with (for 2D, with ) and (for 3D, and 2D with ), where .

## 3 Analysis of numerical results

### 3.1 Flow structures

Instantaneous contours of velocity and vorticity are provided in figures 1(a),(b) and figure 2(a) for 2D and 3D domains, respectively. Specifically, the effects of decreasing the activity coefficient from to in the 2D simulations are an increase in the integral length (computed from the ensemble-averaged kinetic-energy spectrum) from to , a decrease in the Reynolds number from to to , and a decrease in the dissipation from to . The resulting low-activity flow has a less dense pattern of flow structures, as observed by comparing the main (high activity) and inset (low activity) frames in figure 1(a) and noticing that the normalizing length is independent of activity. It is noteworthy that, for the same activity coefficient , moving from two to three dimensions translates into an increase in the integral length from to , a decrease in the Reynolds number from to , and a decrease in the dissipation from to .

The spatial variations in the nematic tensor are central to the generation of vorticity. This is easily observed by taking the curl of (4), namely

(6) |

where is the vorticity and is the permutation tensor. In 2D, the vortex stretch term is exactly zero and the dominant mechanism of vorticity generation is the curl of the divergence of the active stresses. The structures of vorticity, which are shown in figure 1(b), are different from the classic round vortices observed in high-Reynolds number 2D isotropic flows. Instead, vortical structures attain here band-like shapes, which are closely related to thinner elongated regions referred to as walls, where the magnitude of the nematic order tensor becomes small, as shown by the solid contours of overlaid on the director field (largest eigenvector of the tensor) in figure 1(b). The walls are characterized by bend deformations in the in-plane director field in figure 1(c) (note that the out-of-plane components in these 2D simulations are zero), which separate nematically-aligned regions () across interstitial isotropic states (), and are typically much thinner than the hydrodynamic structures of velocity and vorticity. The resulting topological defects, depicted by circles and triangles in figure 1(b,c) (see inset), represent singular, disordered regions of strong vorticity generation that are created and annihilate in pairs while propagating in the flow in a complex manner (Doostmohammadi et al., 2016a) that is beyond the scope of the present study.

In 3D, the vortex-stretch term in (6) represents a much smaller contribution than the active stresses because of the low Reynolds numbers involved. As observed in figure 2(a), the resulting vortical structures in 3D are elongated as well but smaller than the velocity ones. The 3D mechanisms of vorticity generation are mostly unknown since the description of topological defects in active nematics is not well understood. Further insight into rotational and straining components of the 3D flow field can be gained by examining the velocity-gradient invariants and , which are expedient to classify flow structures. In contrast to typical scatter plots for high-Reynolds turbulence where most of the activity is in the upper-right and lower-left quadrants, figure 2(b) shows that for active flows straining is predominant. As a result, the marginal PDF of is heavily skewed toward negative values (skewness ). The straining is caused by the cumulative effect of the stresslets from the extensile active particles (i.e., see flow sketch in figure 2(c)).

### 3.2 PDF moments of flow variables

The PDFs of velocity , velocity gradient , vorticity and magnitude of the nematic order are shown in figure 3, and some of their moments are listed in Table 1. The PDFs of velocity and vorticity have nearly Gaussian flatness, while the largest skewness of the velocity gradient is reached in the 2D high-activity case and equals . In all cases, the vorticity flatness and the velocity-gradient skewness are smaller than typical values observed in high-Reynolds turbulence (i.e., and , respectively). In 2D, small activities favor sub-Gaussian flatness for vorticity and velocity along with decreasing skewness of the velocity gradient. Note, however, that for the same activity coefficient the 3D case leads to a smaller skewness of the velocity gradient, which suggests that the lesser spatial confinement plays a role in the development of fluctuations. Additionally, small values of the nematic-order magnitude, which correspond to isotropic behavior and vorticity generation, are statistically favored by large activities.

2D () | 2D () | 3D () | |
---|---|---|---|

Velocity () flatness | 2.84 | 3.14 | 2.78 |

Velocity derivative () skewness | -0.03 | -0.10 | -0.02 |

Vorticity () flatness | 2.85 | 3.21 | 3.12 |

### 3.3 Intermittency analysis

Albeit small, intermittency may not be entirely ruled out in these systems, as suggested by a narrow-band filtering analysis of the results based on a discrete db-4 wavelet decomposition of the velocity and vorticity fields. This is shown in figure 4, which provides the scale-dependent flatness of the PDFs of the direction-averaged velocity and vorticity wavelet coefficients, and , normalized with their corresponding standard deviations and , with being a scale index that ranges from (equivalent to a length scale ) to (equivalent to ) and is related to a representative wavenumber as . In these simulations, (for 2D cases) and (for 3D cases). The inverted caret denotes the wavelet transform, e.g. , where is a positive-integer direction index ( and in 2D and 3D, respectively), are scale-dependent wavelet grids where the wavelets are collocated, with , , , and are wavelet basis functions that are here taken to be tensor products of orthonormal one-dimensional db-4 wavelets with four vanishing moments. The reader is referred to Meneveau (1991) and Schneider & Vasilyev (2010) for general applications of wavelets in turbulent flows, and to Nguyen et al. (2012) for a scale-dependent flatness analysis of homogeneous isotropic turbulence similar to the one performed here.

While the 2D fields remain nearly Gaussian at all scales, with a slight increase for the vorticity flatness observed in the largest scale, figure 4(b) shows that the 3D fields contain significant intermittency in the small scales, as indicated by the strong increase in with (main panel) and by the increasingly longer tails in the PDFs of the wavelet coefficients as the length scale increases (inset). Nonetheless, the energetic content of these small scales and their associated intermittent motion is small in all cases. This can be understood by noticing the rapid decay of the Fourier kinetic-energy and enstrophy spectra, and , as increases, as observed in figure 5(a,b,d,e). Specifically, the kinetic-energy spectra decay with a slope that decreases from to as the activity is increased ten-fold. As a result, the enstrophy spectra reach a maximum at the integral scale and decay rapidly thereafter, indicating that the small-scale gradients bear vanishing energy. This detracts dynamical relevance from the increased intermittency observed in the 3D small scales and sets a fundamental difference between these flows and high-Reynolds number turbulence; in the latter, the kinetic-energy spectra decays at a slower rate and the small-scale vorticity intermittency is energetic.

It is also of interest to note that the characteristic wavenumber where the nematic-order fluctuation energy spectra (defined such that the area under the curve is ) reach a maximum is approximately one decade smaller than the wavenumber corresponding to the maximum enstrophy spectra in the 2D low-activity case, as shown in figure 5(c). The distance between peaks decreases with increasing activity and when moving to three dimensions, as observed in figure 5(f). The wavenumber of maximum decreases as the activity increases and its value is closer to than to , which spectrally illustrates the fine structure of the nematic-order field compared to the coarser velocity field.

### 3.4 Spectral energy-transfer analysis

As discussed in §2.2, a crucial role in the dynamics is played by the divergence of the active stress . Because of the small involved, it is anticipated that the spectral transfer of the active energy is locally dissipated by viscosity, since the convective inter-scale transfer is a mechanism of secondary importance. Although there may exist additional triadic interactions resulting from (5) that could transport energy across scales, the 3D results provided in figure 6 for active (), convective (), and viscous () wavelet-based spectral energy-transfer fluxes support the view that locality may dominate the transfer. Specifically, these fluxes describe the rate at which the spatially averaged, wavelet spectral kinetic-energy density, , is transferred across scales. In particular, and the analogous and are shown and compared to their Fourier-based counterparts in figure 5(d-f). In this formulation, is a discrete wavenumber shell, and the subindexed bracketed operator represents spatial averaging over scale-dependent wavelet grids .

Upon wavelet-transforming the momentum equation in (4), multiplying by and summing over , the spectral-energy equation is obtained. Here, the source term represents the sum of spectral energy-transfer fluxes created by each term on the right-hand side of the momentum equation in (4). For any force , the corresponding spectral flux is given by , with and indicating, respectively, inflow and outflow of energy at a given . Note that for , for , and for , which satisfy , and .

Figure 6(a) indicates that the active stress act as a kinetic-energy source at all scales on spatial average, with the maximum mean of occurring at scales of the same order as the integral length. Conversely, the viscous flux is a sink of kinetic energy and has a trend that is exactly opposite to , as shown in Figure 6(b), which indicates that the active energy is mostly dissipated locally in spectral space by viscosity. The spatial localization of the transfer, which is illustrated by the unbracketed versions of , , and , is represented by the variabilities of the PDFs shown in figures 5(d-f) and 6. Specifically, the PDFs in figure 6 reveal that the viscous transfer flux is spatially correlated with the active one (correlation coefficient at ), suggesting that upon deployment the active energy is dissipated mostly locally also in physical space.

The physical picture implied by figure 6 provides no evidence for an energy cascade in momentum where the sink and main source of mean kinetic energy could act in disparate ranges of scales interacting through a crossing long-range mechanism. This is in contrast to high-Reynolds number turbulence and its clear separation of scales between the large-scale forcing range and the small-scale molecular-dissipation range. These conclusions could however be different for the nematic-order energy, in that the transport description of the latter is highly non-linear and involves cross-triadic terms with the velocity as in (2). These aspects will be the subject of future research.

## 4 Conclusions

The multi-scale statistical analysis of DNS presented in this study provides quantitative comparisons between active and classic turbulent flows beyond superficial visual similarities, demonstrating clear distinctions in the intermittency characteristics and mechanisms of inter-scale energy transfer. It is shown that increasing activities lead to increasingly packed and dissipating structures that have increasingly larger departures from Gaussian statistics. For the same activity, the 3D flow has a larger integral length and smaller kinetic energy compared to its 2D counterpart. A velocity-gradient invariant analysis of the 3D flow indicates that straining structures dominate the topology as a collective result of the embedded stresslets induced by each individual extensile active particle. A wavelet-based, scale-dependent flatness analysis shows the occurrence of intermittency in the small scales, particularly in the 3D vorticity field. However, the spectral energy content associated with the small-scale velocity gradients is small in all cases. The work of the active stress is spectrally deployed near the integral scales and mostly dissipated locally by viscosity in both physical and spectral spaces.

### Acknowledgments

This investigation was performed during the 2016 CTR Summer Program at Stanford University. The authors are grateful to Maxime Bassenne, Dr. Jeonglae Kim, and Profs. Marie Farge and Kai Schneider for useful discussions.

### References

- Bratanov, V., Jenko, F. & Frey, E. 2015 New class of turbulence in active fluids. PNAS 112, 15048–15053.
- Doostmohammadi, M. F., Thampi, S. P. & Yeomans, J. M. 2016a Defect-mediated morphologies in growing cell colonies. Phys. Rev. Lett. 117, 048102.
- Doostmohammadi, A., Adamer, M. F., Thampi, S. P. & Yeomans, J. M. 2016b Stabilization of active matter by flow-vortex lattices and defect ordering. Nat. Commun. 7, 10557.
- Dunkel, J., Heidenreich, S., Drescher, K., Wensink, H. H., Bär, M. & Goldstein, R. E. 2013 Fluid dynamics of bacterial turbulence. Phys. Rev. Lett. 110, 228102.
- De Gennes, P. G. & Prost, J. 1995 The Physics of Liquid Crystals. Oxford Univ. Press.
- Edwards, B., Beris, A. N. & Grmela, M. 1991 The dynamical behavior of liquid crystals: a continuum description through generalized brackets. Molecular Crystals and Liquid Crystals 201 (1), 51–86.
- Giomi, L. 2015 Geometry and topology of turbulence in active nematics. Phys. Rev. X 5, 031003.
- Meneveau, C. 1991 Analysis of turbulence in the orthonormal wavelet representation. J. Fluid Mech. 232, 469-520.
- Nguyen, R. V. Y., Farge, M. & Schneider, K. 2012 Scale-wise coherent vorticity extraction for conditional statistical modeling of homogeneous isotropic two-dimensional turbulence. Physica D 241, 186-201.
- Ottino, J. M. 1990 Mixing, chaotic advection, and turbulence. Annu. Rev. Fluid Mech. 22, 207–253.
- Sanchez, T., Chen, D. T. N., DeCamp, S. J., Heymann, M. & Dogic, Z. 2012 Spontaneous motion in hierarchically assembled active matter. Nature 491, 431–434.
- Saw T. B., Doostmohammadi A., Nier V., Kocgozlu L., Thampi S., Toyama Y., Marcq P., Lim C. T., Yeomans J. M. & Ladoux B. 2017 Topological defects in epithelia govern cell death and extrusion. Nature 544, 212–216.
- Schneider K. & Vasilyev O. V. 2010 Wavelet methods in computational fluid dynamics. Annu. Rev. Fluid Mech. 42, 473-503.
- Simha, A. & Ramaswamy, S. 2002 Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles. Phys. Rev. Lett. 89, 058101.
- Thampi, S. P., Golestanian, R. & Yeomans, J. M. 2013 Velocity correlations in an active nematic. Phys. Rev. Lett. 111, 118101.
- Wensink, H. H., Dunkel, J., Heidenreich, S., Drescher, K., Goldstein, R. E., Lowen, H. & Yeomans, J. M. 2012 Meso-scale turbulence in living fluids. PNAS 109, 14308–14313.