## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# New class of turbulence in active fluids

Edited by Alexander J. Smits, Princeton University, Princeton, NJ, and accepted by the Editorial Board October 29, 2015 (received for review May 12, 2015)

## Significance

It is widely appreciated that turbulence is one of the main challenges of modern theoretical physics. Whereas up to now, most work in this area has been dedicated to the study of Navier–Stokes flows, numerous examples exist of systems that exhibit similar types of spatiotemporal chaos but are described by more complex nonlinear equations. One such problem of quickly growing scientific interest is turbulence in active fluids. We find that such systems can exhibit power-law energy spectra with nonuniversal exponents as a result of nonlinear self-organization, defining a new class of turbulent flows.

## Abstract

Turbulence is a fundamental and ubiquitous phenomenon in nature, occurring from astrophysical to biophysical scales. At the same time, it is widely recognized as one of the key unsolved problems in modern physics, representing a paradigmatic example of nonlinear dynamics far from thermodynamic equilibrium. Whereas in the past, most theoretical work in this area has been devoted to Navier–Stokes flows, there is now a growing awareness of the need to extend the research focus to systems with more general patterns of energy injection and dissipation. These include various types of complex fluids and plasmas, as well as active systems consisting of self-propelled particles, like dense bacterial suspensions. Recently, a continuum model has been proposed for such “living fluids” that is based on the Navier–Stokes equations, but extends them to include some of the most general terms admitted by the symmetry of the problem [Wensink HH, et al. (2012) *Proc Natl Acad Sci USA* 109:14308–14313]. This introduces a cubic nonlinearity, related to the Toner–Tu theory of flocking, which can interact with the quadratic Navier–Stokes nonlinearity. We show that as a result of the subtle interaction between these two terms, the energy spectra at large spatial scales exhibit power laws that are not universal, but depend on both finite-size effects and physical parameters. Our combined numerical and analytical analysis reveals the origin of this effect and even provides a way to understand it quantitatively. Turbulence in active fluids, characterized by this kind of nonlinear self-organization, defines a new class of turbulent flows.

Despite several decades of intensive research, turbulence—the irregular motion of a fluid or plasma—still defies a thorough understanding. It is a paradigmatic example of nonlinear dynamics and self-organization far from thermodynamic equilibrium also closely linked to fundamental questions about irreversibility (1) and mixing (2). The classical example of a turbulent system is a Navier–Stokes flow, with a single quadratic nonlinearity, well-separated drive and dissipation ranges, and an extended intermediate range of purely conservative scale-to-scale energy transfer (3). However, many turbulent systems of scientific interest involve more general patterns of energy injection, transfer, and dissipation. A fascinating example of these kinds of generalized turbulent dynamics can be observed in dense bacterial suspensions (4). Although the motion of the individual swimmers in the background fluid takes place at Reynolds numbers well below unity, the coarse-grained dynamics of these self-propelled particles display spatiotemporal chaos, i.e., turbulence (5⇓–7). Nevertheless, the correlation functions of the velocity and vorticity fields display some essential differences compared with their counterparts in classical fluid turbulence (8, 9). Moreover, the collective motion of bacteria in such suspensions exhibits long-range correlations (10), appears to be driven by internal instabilities (11), and depends strongly also on physical parameters like large-scale friction (12). Such results challenge the orthodox understanding of turbulent motion and call for a detailed theoretical investigation. There also exist many other systems with similar characteristics, including flows generated by space-filling fractal square grids (13) , turbulent astrophysical (14) and laboratory (15, 16) plasmas, and chemical reaction–diffusion processes (17).

In the present work, we study—numerically as well as analytically—the spectral properties of a continuum model that has recently been suggested as a minimal phenomenological model to describe the collective dynamics of dense bacterial suspensions (4, 18, 19). A basic assumption of the model is that at high concentrations the dynamics of bacterial flow may be described as an incompressible fluid obeying the following equation of motion for the velocity field **1** amounts to a straightforward multidimensional generalization of the Kuramoto–Sivashinsky (KS) equation that has found application in describing magnetized plasmas (20, 21), chemical reaction–diffusion processes (22, 23), and flame front propagation (24, 25). It is widely regarded as a prototypical example of “phase turbulence.” (26) As a hallmark, if both kinetic parameters are positive (*k*, similar to other paradigmatic models of nonlinear dynamics, e.g., the Swift–Hohenberg model (27). For active systems this feature emulates energy input into the bacterial system through stress-induced instabilities (11). The growth of these linearly unstable modes is limited by nonlinear and dissipative terms. The main dissipation mechanism in Eq. **1** is mediated through the cubic nonlinearity on the right-hand side, **1** serves as a simple but generic test case to address some of the fundamental questions in the field of active turbulence. Via an appropriate choice of parameters, one can describe several different physical systems as explained in more detail in Table S1.

First and foremost, the similarities and differences between low and high Reynolds number turbulence remain to be elucidated. In particular, there is still a lack of understanding of the energy flow between different length scales. Here, we address the above questions by a systematic analysis of the turbulent features of Eq. **1**, combining numerical and analytical approaches, and we give a comprehensive picture of the spectral energy balance facilitating the understanding of the interactions among different spatial scales. Furthermore, extensive numerical simulations confirm the existence of a spectral power law at the largest scales of the system with its steepness depending on the parameters of the system (both of the linear and of the nonlinear terms in Eq. **1**). The form of the turbulent energy spectrum is an important quantity related to the frictional drag between the system and the surrounding walls (31, 32). In the present work, insight into the remarkable feature of a variable spectral exponent is gained by analyzing the role of the different terms in the equation for the spectral energy balance. As expected for a 2D incompressible fluid, there exists an inverse flow of energy from intermediate to large scales (33). Nevertheless, in contrast to classical, fully developed 2D Navier–Stokes turbulence, there is no inertial range characterized by a constant energy flux. Instead, we find that at large scales the nonlinear frequency corresponding to the Navier–Stokes energy flux is constant for the whole range characterized by spectral self-similarity. This differs fundamentally from the classical Navier–Stokes case, where this nonlinear frequency, the inverse of the nonlinear eddy turnover time, is a function of wave number and energy. In the model at hand, this energy flux is balanced by a linear dissipation/injection and a cubic dissipation term. For the latter, we derive an analytic approximation that compares very favorably with the numerical results and allows for an analytic closure predicting the type of dependence of the power law on the model parameters that is also confirmed numerically.

## Results

We have studied the 2D version of the continuum model defined by Eq. **1** both analytically and numerically. Our computational approach relies on a pseudospectral code where the linear terms are computed in Fourier space and the nonlinearities in real space. The details of this procedure and the necessary normalization are described in *SI Text*, section S1. All numerical results reported in this paper use a resolution of 1,024 effective Fourier modes in each direction, unless stated otherwise. A typical velocity is given by **1**, one reads off the wave number of the fastest growing mode, *β* and

### Spectral Analysis.

For the analysis of the flow of energy between different spatial scales mediated by the various terms in Eq. **1** we use a Fourier decomposition of *SI Text*, section S1). *SI Text*, section S1. The ensuing spectral energy balance equation reads*α*) either injects (**1** result from the nonlinear terms, i.e., advection term and cubic nonlinearity. They couple different wave numbers and provide a flow of energy in spectral space that (on average) balances the local injection or dissipation. The different terms in Eq. **2**, obtained from a numerical solution of Eq. **1**, are shown in Fig. 2. We have averaged over nearly *k*). Note that the advective nonlinearity (green curve) is positive for small *k* but negative for intermediate *k* and, thus, transports energy from small to large scales. This inverse energy flow is characteristic of 2D turbulent systems and is due to the constraint of enstrophy conservation (34). In the present context, it takes energy from the intermediate wave numbers where the *k*.

### Spectral Shell Decomposition.

To further assess the energy transfer among different length scales, we divide the spectral space into circular shells *SI Text*, section S2). Moreover, we introduce the projection operator

Applying **1** leads to an evolution equation for the energy *SI Text*, section S1).

The terms *I* and *J* (*SI Text*, section S2); i.e., summing over both indexes gives zero. This shows that (in an incompressible system) the Navier–Stokes nonlinearity neither injects nor dissipates energy but only redistributes it among the different shells *A*. In addition to verifying the antisymmetry, this also illustrates the direction of energy transfer in spectral space. There is a combination of forward and inverse energy flows. At intermediate wave numbers, there is mainly a local forward energy flux; see the areas next to the diagonal in Fig. 3*A*, where red above the diagonal and blue below it indicate a flow from smaller to larger *k*. Additionally, there is also a nonlocal inverse energy flow dominating at small wave numbers, represented by the smaller side branches in Fig. 3*A*. The green curve in Fig. 2 represents the cumulative effect of the 2D structures seen in Fig. 3*A*. The Navier–Stokes nonlinearity extracts energy from the intermediate wave numbers (negative contribution) and supplies it to both smaller (inverse cascade) and larger (forward cascade) wave numbers. The contribution of the cubic nonlinearity, on the other hand, is symmetric (*SI Text*, section S2). Because every second-rank tensor can be uniquely decomposed into a symmetric and an antisymmetric part, *B* clearly show that the entries of *B*, *Inset* resembles closely the blue line in Fig. 2. Moreover, the diagonal entries are negative, indicating the dissipative nature of the cubic nonlinearity. This feature together with the different physical interpretation of the cubic term represents the central result that the shell-to-shell decomposition yields. Both aspects are essential for the cubic interaction and should be captured by a successful closure approximation.

### Cubic Damping Term.

To make progress beyond a numerical analysis, we seek an approximate solution for the stationary state of the energy spectrum. The analysis is complicated by the fact that the right-hand side of Eq. **2** involves third- and fourth-order velocity correlation functions, **1** as explained in *SI Text*, section S3), a natural way to approach the cubic damping term in Eq. **2** is via the quasi-normal approximation (35), also known as the Millionshchikov hypothesis (36). According to it, third-order correlations, e.g., *SI Text*, section S3. Hence, the approximation of the cubic damping term in Eq. **2** is directly proportional to the energy spectrum *B*, showing that the diagonal terms are the dominant ones in *α* (necessarily positive) below which the dissipation due to friction is insufficient and cannot balance the energy that accumulates at the large scales as a result of the inverse energy flow in 2D Navier–Stokes systems.

### Advective Nonlinearity.

In contrast to the cubic damping term, the advective nonlinearity in Eq. **2** produces an expression that involves third-order correlations, meaning that the quasi-normal approximation is not directly applicable. Formulating an evolution equation for the third-order correlation leads to the known hierarchy problem, which here, due to the presence of the cubic term, would be even more convoluted. Such a hierarchical scheme can, nevertheless, lead to a closed system of equations after applying the Millionshchikov hypothesis, but the resulting system of equations is highly complicated and tractable only numerically.

Because our goal here is to arrive at an analytical approximation for the energy spectrum at small wave numbers, we choose a more heuristic approach. As already discussed, the advective nonlinearity redistributes energy only among the different modes. This implies an energy flux in spectral space, defined as *k*. Because *k* in the cascade range of Navier–Stokes turbulence, the main contribution to the integral comes from the part of the integrand around *k*, which relates to the locality of the classical energy cascade. Eq. **8** together with the integral expression for *k* as seen in Fig. 3*A*. In addition, in classical turbulence models large spatial scales are more energetic than smaller ones, giving the former the potential to shear and distort the latter. For Eq. **1**, however, the energy spectrum *k* up to some maximum and then decreases again; see the red curve in Fig. 2. Hence, for the spectral region we are interested in, the larger scales are not able to shear the smaller ones. In summary, our investigation of the advective nonlinearity in this model shows that at small wave numbers there is a distinct constant frequency

### Variable Spectral Exponent.

In the statistically stationary state, time-averaging Eq. **2** yields zero on the left-hand side,**8** with constant **10** shows that at small wave numbers (*δ* of this power law is not universal but depends (directly and indirectly) on various system parameters. Qualitatively, a stronger dissipation, i.e., a positive *α* and a higher factor of **1** is presented for two different values of *α*. In both cases, the system exhibits clear power-law spectra over more than one order of magnitude in wave-number space, and it is evident that our model predicts the correct qualitative dependence of the spectral exponents. A quantitative test of our semianalytical result can be undertaken by carrying out numerical simulations for different values of *α*. We note in passing that such a parameter scan requires that there are always enough instabilities to drive the turbulence and that statistical homogeneity and isotropy are ensured. The linear growth rate of the most unstable mode equals *α* once **1** tries to destroy the statistical isotropy of the system. Thus, the energy injected by the *α* term must be considerably smaller than that injected by the *α*. The result from such a parameter scan of the numerical solution of Eq. **1** is displayed in Fig. 6, where every point is obtained by fitting a power law on the left end of the energy spectrum. The data from our investigation show a linear dependence of the slope *δ* on the parameter *α*, which agrees with the expression for *δ* provided by our model. Further numerical simulations indicate that the dependence of the slope on the strength of the cubic interaction *β* is qualitatively the same but quantitatively much weaker. This can be due to the factor *δ*. Stability analysis shows that for *β*.

## Summary and Conclusions

In the present work, we investigated the properties of a continuum model describing the turbulent motion of active fluids driven by internal instabilities. In addition to the convective nonlinearity of Navier–Stokes type, the model contains a cubic nonlinearity. Analytical and numerical considerations revealed that at large scales, the latter behaves like an Ekman damping with a frequency that is set by the system self-consistently. The system displays power-law energy spectra even in the absence of an inertial range, but the spectral exponents depend on the system parameters. These properties should be observable in laboratory experiments.

How do these findings fit into a broader perspective on turbulence? In Navier–Stokes turbulence, the dynamics are characterized by an inertial range that is dominated by a single nonlinear term and free of energy sources/sinks, displaying universal properties. Several turbulence models in the literature deviate from this standard picture in that they introduce multiscale forcing and/or damping with a power-law spectrum, thereby removing the inertial range (in a strict sense) (38⇓–40). It can be shown, however, that, in general, this modification really affects the system only at very small or very large scales, i.e., in the asymptotic limit (41). If and only if the power-law exponent is such that the linear forcing/damping rates scale exactly like the nonlinear energy transfer rates, does the forcing/damping affect the entire scale range, inducing nonuniversal behavior (16, 42, 43).

The physical system discussed in this paper is fundamentally different, however. Here, the existence of a second nonlinearity provides additional freedom, such that the system is able to self-organize into such a critical state (characterized by a scale-by-scale balance between linear forcing/damping rates and nonlinear transfer rates), without the need for external fine-tuning. The observed nonuniversal behavior is a natural consequence of this feature. These properties define a new class of turbulence.

## SI Text

## S1. Normalization and Numerical Methods

### Normalization.

For the sake of completeness, we discuss here briefly the normalization of the quantities used and describe the numerical implementation of the continuum model. Written in dimensional units, Eq. **1** reads*ρ* denotes density of the fluid and is considered constant. In contrast to ref. 4, we have omitted the term *γ* of the spectral modes given by *β* of

Eq. **1** possesses five different parameters, which makes a thorough investigation of its complete parameter space a very tedious task. To bring some structure into this highly dimensional parameter space, we present some limit cases in Table S1. It is important to note that for **1** do not reach a stationary state and the kinetic energy of the velocity field grows without a limit. Thus, this case is of no physical relevance.

### Numerical Methods.

Numerical solutions of Eq. **1** have been obtained via the pseudospectral approach (44), which is common in computational fluid dynamics. According to it spatial derivatives are computed in Fourier space where they reduce to simple multiplication. The nonlinear terms, on the other hand, are computed in real space because they correspond to a convolution in Fourier space, which is a computationally expensive procedure. The time evolution is computed numerically by the exponential time differencing scheme first developed in ref. 45 and later improved in ref. 46, where a fourth-order Runge–Kutta method has been used. The underdetermined system of Eq. **1** is completed by the incompressibility condition ^{^} symbol. It should be clear from the corresponding argument whether we mean the real-space function or the Fourier component. According to the definition given above, the Fourier coefficients are determined as**1** the

The solution of Eq. **1** in the part of parameter space we are interested in exhibits rapid fluctuations typical for turbulent systems. The same applies also for quantities like the energy of a given mode or the total energy of the system. In the turbulence literature the so-called “ensemble average” is defined, which corresponds to an average over many possible realizations of the flow due to different initial conditions. However, such an average procedure is computationally extremely demanding. In the case of a system that reaches a statistically stationary regime, we use a time average instead, denoted by *t* the integral becomes a sum. Under the assumption of ergodicity the definition above yields the same result as an ensemble average.

In ref. 4, where the continuous model we study here has been introduced, the main goal has been to compare its results with experimental measurements and with the findings based on more basic self-propelled rods (SPR) models. This determined the size of the real-space domain used in the simulations. In this article, however, we aim for a more fundamental examination of the turbulent dynamics produced by Eq. **1**. Because we are interested in the spectral range at small wave numbers, we investigate the convergence of the numerical results when the box size becomes larger. The latter is inversely proportional to the smallest nonzero wave number we have, meaning that a larger real-space domain allows for a better representation of the small-*k* part of the energy spectrum. Numerical simulations of Eq. **1** for different domain sizes yield the energy spectra displayed in Fig. S1. The box size and number of points representing numerically the real-space domain that have been used for obtaining the red curve are nearly the same as those used in ref. 4 and the slope of *k* points such a grouping can lead to notable ambiguity in the value of *k* assigned to the whole group and, thereby, influence the form of the curve. Therefore, the slope of a power law calculated by using only the first few points will be very sensitive to such numerical details. The blue curve, on the other hand, has been obtained with a much larger domain size and number of points (effectively *k* with a prominent power law of nearly

## S2. Symmetry of the Nonlinear Transfer Terms

Generally speaking, nonlinear terms tend to play a key role in determining the behavior of complex systems. This is also true here. In the main part of the paper, these terms were studied quantitatively by introducing specific projection operators in spectral space as defined in Eq. **4**. One can easily verify that *x* or the *y* axis. This ensures that every shell contains an appreciable number of modes and thereby reduces statistical fluctuations.

Applying **1** gives us an evolution equation for the filtered velocity field **S7**; i.e.,

Applying the velocity decomposition on the term stemming from the cubic nonlinearity in Eq. **1** leads to**S13** and **S14** we can read off the corresponding the shell-to-shell coupling terms; namely**6a** and **6b** in the main text, representing merely their real-space formulation.

Due to the symmetry of the scalar product *I* and *J*; i.e., *I* and *J*. Denoting for brevity

## S3. Derivation of Eq. 7 from the Millionshchikov Hypothesis

As discussed in the main text, the quasi-normal approximation rests on the Millionshchikov hypothesis (36) according to which the velocity correlations of odd order are still nonzero, but the even-order correlations can be expressed approximately as the sum of all possible products of second-order correlations. Thus, for the fourth-order correlation in Eq. **2** we can write*i*, *j*, and *ł*. A preliminary computation that proves to be useful in this regard yields**7**. For the second part of Eq. **7**, one needs to relate the scalar correlation function *θ*, we can write**2**.

In the turbulence literature it is often tested how close the statistics of the velocity field are to the Gauss distribution by means of the two-point velocity increments *r* between the two points; i.e., *P* of the velocity increments. (Strictly speaking, this gives us the probability density function.) A numerical computation of *P* (blue circles) for the case of *Upper Row* displaying *Lower Row* displaying *A* and *σ* being the fit parameters. As seen in Fig. S2, the distribution of velocity increments is very close to a Gaussian at large scales, which supports the use of the quasi-normal approximation. In addition, note that all of the distributions in Fig. S2 have nearly the same parameters *A* and *σ*. To connect the results shown here to the energy spectra in wave-number space let us remark that the separations displayed in Fig. S2 correspond to *k* range studied in this work, because, for comparison, the peak of the energy spectrum lies around *σ*. Here it is considerably lower, meaning that the width of the distribution is smaller. The reason for this lies in the fact that for

## Acknowledgments

V.B. thanks A. Bañón Navarro, M. Oberparleiter, and D. Groselj for helpful discussions. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/European Research Council Grant Agreement 277870. It was also supported by the German Excellence Initiative via the program “NanoSystems Initiative Munich” (NIM) and the Deutsche Forschungsgemeinschaft (DFG) in the context of SFB 863 “Forces in Biomolecular Systems” (Project B02).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: jenko{at}physics.ucla.edu.

Author contributions: F.J. designed research; V.B., F.J., and E.F. performed research; and V.B., F.J., and E.F. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission. A.J.S. is a guest editor invited by the Editorial Board.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1509304112/-/DCSupplemental.

## References

- ↵.
- Xu H, et al.

- ↵
- ↵.
- Kolmogorov AN

- ↵.
- Wensink HH, et al.

- ↵
- ↵.
- Zhou S,
- Sokolov A,
- Lavrentovich OD,
- Aranson IS

- ↵
- ↵
- ↵.
- Thampi SP,
- Golestanian R,
- Yeomans JM

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Kuramoto Y,
- Tsusuki T

- ↵
- ↵.
- Sivashinsky GI

*Acta Astron*4:1177–1206 - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Tran T, et al.

- ↵
- ↵
- ↵.
- Monin AS,
- Yaglom AM

- ↵.
- Millionshchikov MD

- ↵
- ↵Sain A, Manu, Pandit R (1998) Turbulence and multiscaling in the randomly forced Navier-Stokes equation..
*Phys Rev Lett*81:4377–4380 - ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Gottlieb D,
- Orszag SA

- ↵
- ↵
- ↵.
- Canuto C,
- Hussaini MY,
- Quarteroni A,
- Zang ThA

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Physics