Correlation lengths in hydrodynamic models of active nematics


1.125 0.8cm \titlespacing*


0pt4pt4pt \titlespacing*




2 Introduction

Active systems are continuously driven out of equilibrium by energy injected at the local scale, resulting in collective motion at the large scale. Examples include bacterial colonies, in vitro extracts of cytoskeletal filaments and associated motor proteins, and monolayers of vibrated granular matter 1. Much interest has recently focused on active nematics - systems composed of active units with head-tail symmetry that can order into states with nematic liquid crystalline order, with the activity giving rise to a rich variety of collective phenomena. These include spontaneous laminar flow 2, 3, 4, large density fluctuations 5, 6, 7, pattern formation 8, 9, spontaneous unbinding of topological defects 10, 11, 12, 13, and active turbulence 14, 15, 16, 17, 18.

Theoretical interest in active nematics has been fueled by the observation of spontaneously flowing and turbulent states in suspensions of microtubule-kinesin bundles confined at an oil-water interface 19, 20, 21. Other recent experimental realizations were obtained by immersing living swimming bacteria (specifically, E. coli) in a lyotropic liquid crystal 22 and by plating dense layers of fibroblasts on substrates 23. Microtubule bundles and E. coli exert active extensile force dipoles on their surroundings, while the fibroblasts exert contractile force dipoles. In both systems active stresses couple to orientational order and induce flows and defect unbinding, with qualitatively similar non-equilibrium dynamics, albeit on very different time scales.

All these experiments have indicated that the appearance of spatio-temporal chaos in active nematics is accompanied by the proliferation of topological defects, which may mediate the onset of turbulence even at negligible values of the Reynolds number 11, 14. In active fluids, distortions of the local orientational order produce local shear flows that enhance the orientational deformation, ultimately leading to the unbinding of defect pairs. In nematics, these consist of pairs of disclinations - specific distortions of the orientational order that are the signature of the underlying broken symmetry. It was additionally shown 11 that the dynamics of topological defects in active systems depends on the nature of the forcing, i.e., whether the force dipole is contractile or extensile. This result provides a useful criterion for inferring the nature of active forces in systems with nematic symmetry 11, 16, 24.

In spite of much theoretical work, discrepancies still exist in the literature over the nature of the characteristic length scales that control the spontaneous proliferation and annihilation of topological defects, and the resulting dynamics in the so-called turbulent state. In particular, the dependence of such length scales on the strength of the active forcing, , remains unclear. A recent numerical study by Giomi17 examined the statistics of the activity-driven turbulent phase in two-dimensional nematic films by measuring the distribution of vortex sizes for a selection of activities. This work provided evidence that the key physics is determined by a single active length scale, , proportional to . In contrast, in a closely related work, but on a different continuum model of a quasi- nematic, Thampi et al. performed a detailed study that measured several orientational and hydrodynamical correlation lengths, suggesting that the length scale of structure in the fluid instead scales as .

An important aim of this work is to provide a unified understanding of these previous, apparently conflicting reports. In order to do so, we consider two different but related models of active nematic liquid-crystal hydrodynamics that are commonly used in the literature. The dynamics of these models is compared using two independently developed numerical codes. By varying the key dimensionless parameters over several decades, we obtain data to support the conjecture that, in both models, the mean defect spacing in the regime of full developed active turbulence is set by the length scale , defined by the balance of active and elastic stresses. Here is the magnitude of the active stress (commonly referred to as the activity) and parametrizes the free energy penalty that results from spatial variations in the director field 3, 4, 1, 13, 17, 25. We show that this result holds for both extensile and contractile systems, in both the flow aligning and the flow tumbling regimes. This active length scale also controls the onset of spontaneous laminar flow in an active film, a phenomenon that has been referred to in the literature as the spontaneous flow instability 2. Our study provides the first explicit demonstration that distinct constitutive models produce the same emergent length scale, i.e., they both produce quantitatively consistent scaling relations. We also demonstrate a regime of less highly developed turbulence in which a weaker scaling appears consistent with our numerical data.

In many experimental realizations active nematics are confined to quasi two-dimensional geometries, e.g., on the surface of lipid vesicles 20, in flattened water-in-oil droplets 19, in thin-films 26, or squeezed between parallel glass plates 22. It is worth noting that the presence of confining walls (with no-slip boundary conditions) in the last of these modifies hydrodynamic interactions between active particles 27; of the above experiments, free-standing films may then provide the closest experimental realisation of our system. Clearly, in any numerical study, it is important to define carefully the considered dimensionality. In what follows we denote by the number of dimensions in which the relevant fields (nematic order parameter tensor, fluid velocity, etc.) are allowed to vary; and separately by the number of dimensions in which the nematic director is allowed to develop non-zero components. We shall perform two different studies. In the first we take a strictly two-dimensional model of an active nematic sheet, in which the order parameter tensor is allowed to develop non-zero components only in the directions (); and physical quantities are likewise allowed to vary only in the plane (). In the second study we consider a three-dimensional nematic () but in which all quantities are nonetheless still assumed to be spatially homogeneous in the direction of the layer thickness (). In the latter case the director can in principle point out of the simulated plane, and indeed this effect has been reported in a previous numerical study of active nematics in a cylindrical capillary 28. In that work confinement in the plane of simulation plays an important role whereas our study deliberately focuses on the bulk dynamics where we do not observe out-of-plane motion. Finally it we note that accurately resolving turbulent active hydrodynamics is computationally demanding, especially when performing large sweeps of parameter space. This necessarily restricts us to a 2D study, as is the case with many other studies of active turbulence 3, 14, 17.

The paper is structured as follows. In Section 3 we define the equations of motion for both models, and outline the parameter ranges that we explore for each. In Section 4, we define the observable length scales that can be used to characterize the fluid structure, and discuss the physical reasoning behind their definition. The results of our study are presented in Section 5. In Section 6 we provide a comparison with other works and offer our conclusions.

3 Models

In spatial dimensions we consider an incompressible uniaxial active nematic liquid crystal with a director that can orient in dimensions, with in our two respective studies. The nematic orientational order is parametrized by a symmetric and traceless tensor field , where is the order parameter magnitude and the director is a headless unit vector that characterizes the direction of broken orientational symmetry. The nematic is embedded in an incompressible fluid of constant density, , and constant viscosity, . The fluid velocity field is denoted by . The associated pressure field is determined by the incompressibility condition .

The equations of motion for an active nematic are derived from the well-known hydrodynamic equations for a passive liquid-crystal 25


where is the material derivative and is a rotational viscosity. Here and denote the symmetric and antisymmetric parts of the rate of strain tensor , respectively, with and . For other tensors the transpose, symmetric, antisymmetric and traceless parts are denoted by the superscripts , , and , respectively. For example, .

The relaxation dynamics of the alignment tensor in Eq. 2 is governed by the molecular field, , in which the Landau-de Gennes free energy 25, , is the sum of contributions from a bulk free energy density


and the distortion free energy density


Here and determine the bulk and distortion energy density scales respectively. For simplicity we have adopted the one-elastic constant approximation in the distortion free energy (Eq. 4).

The tensor is defined for an arbitrary tensor as




Here is the Leslie-Ericksen flow aligning parameter, which specifies how the nematic director responds to a shear flow: corresponds to flow-aligning nematics and corresponds to the flow-tumbling regime. (See Appendix A.)

The stress tensor in Eq. 1 is the sum of passive liquid-crystal and active contributions, . The passive part of the stress tensor is given by


In an active nematic there is an additional active stress contribution that arises from the dipolar forces exerted by active particles on their environment. This active stress is , where describes contractile stresses and extensile stresses. In the passive limit , the equations just described reduce to those of a passive liquid-crystal.

So far, the model that we have presented encompasses both of the numerical studies performed. We now outline the specific choices for parameter values and dimensionality made in each of the two numerical studies separately, and discuss how these choices affects the form of the equations. The main difference between the two variations will be the presence of higher order coupling terms in Model II, both in the liquid-crystal stress and in the coupling between orientational order and velocity gradients.

Model I (MI): This model describes a dimensional nematic sheet with a dimensional nematic order parameter. In this case the symmetric second rank tensor has only two independent components and identically. In Model I the mean-field free energy (Eq. 3) has coefficients and , where is a dimensionless parameter that controls the continuous transition from an isotropic to a nematic state, with the transition occurring at . The second term in Eq. 5 is also identically zero; the third term is of a higher order in and can safely be neglected 29, so that . We also exclude the last term in Eq. 7 in Model I. We assume a constant density for which the Reynolds number , where the velocity scale is 5. We choose such that the system is deep in the nematic state, with . According to Eq. 6, the system will be in the flow-aligning regime if and in the flow-tumbling regime for . All results shown below for Model I correspond to (flow-tumbling regime). We also choose for extensile systems () and for contractile systems () to guarantee , a condition that is required in to observe the initial flow instability in the ordered state 13.

Model II (MII): This model considers a dimensional layer of nematic liquid crystal described by the full Landau free energy given in Eq. 3, thereby in principle allowing the director to explore all dimensions. However it still neglects all spatial variations in the direction of the layer thickness, so taking as noted above. In this case the free energy in Eq. 3 sets , and , yielding a first order isotropic-nematic transition at . In the following we choose which places us at the spinodal stability limit of the isotropic phase and well within the nematic state, with . According to Eq. 6, the system will be in the flow-aligning regime for and in the flow-tumbling regime for . In all simulations using Model II we have fixed , corresponding to a flow-aligning system. We consider only extensile systems with this model, i.e., values of . Finally, in MII we take the limit of zero Reynolds number by setting .

parameter description dimensions
Frank constant ( in MI)
energy density scale ( in MI, MII)
rotational viscosity ( in MII)
alignment param.
IN control param.
solvent viscosity ( in MI)
solvent density
box size ( in MII)
Table 1: Summary of the various model parameters and their dimensions. The choices for mass (or equivalently stress ), length and time in each model are also indicated.

The full list of nine parameters (for both models) is given in Table 1. We are free to choose units of mass , length and time , or equivalently of stress , length and time , and we have noted in Table 1 which quantities we chose to set equal to unity in each of the two studies. This leaves six dimensionless groupings that we summarize in Table 2, three of which are fixed throughout. Therefore even though we choose our units differently in the two different simulation studies, all results are presented and compared in a consistent adimensional way between the two models. We choose parameters that produce flow-tumbling behaviour in Model I and flow-aligning in Model II. Note that due to differences in parameter selections, the linear instability thresholds in the two models differ by a factor ( in Model I 13 and in Model II 30). Accordingly, the onset of the turbulent regime in each model is separated by a similar factor, requiring us to explore different ranges of dimensionless activity, as noted in Table 2. In particular, in Model II we explore the transition from small to large activities, whereas Model I focuses on larger activities still (i.e., deeper into the regime of fully developed active turbulence).

parameter description MI value MII value
varied parameters
dimensionless activity
ratio of micro- to macroscopic length scales
ratio of viscosities 0.567
fixed parameters
alignment param. 0.7
IN control param. 2 3
Reynolds number 0
Table 2: Summary of the dimensionless parameters and their values in both models. The velocity scale in our units.

3.1 Numerical details

In order to demonstrate the robustness of our results with respect to numerical implementation, we use two independent codes (one for each model). In each case we perform simulations in a square box of side with biperiodic boundary conditions. The dynamics in Model I is time-integrated numerically on a square grid of points using a fourth order Runge-Kutta method, with a timestep . Gradients of are computed using a finite difference scheme. In Model II is integrated numerically using a Euler time-stepping scheme of timestep in the range on a grid of points (dependent on the magnitude of activity) and gradients of are treated using a semi-implicit Fourier method. In Model II, the velocities are determined instantaneously from the force balance equation, which we solve in Fourier space using a stream function formulation. In Model I, we integrate the Navier-Stokes equation (Eq. 1) with the same scheme used for the order parameter equation to obtain the velocity at every time step. We have verified that our results are quantitatively unchanged upon decreasing the timestep or grid spacing. Both simulations were initialised with a uniform director field orientated within the plane, subject to a small sinusoidal perturbation of magnitude .

4 Characteristic length scales

Irrespective of the specific details of the model used, we expect the resulting dynamics of the active nematic to be controlled by the interplay of key length- and time scales that govern the basic physics.

4.1 Model length scales

An inspection of the hydrodynamic equations and model geometry reveals three underlying length scales. The first is simply the system size . The second arises from balancing the bulk and elastic-distortion free energy terms in Eqs. 3 and 4, to obtain the equilibrium nematic persistence length, which deep in the nematic state is given by 6


This is the length scale over which spatial correlations in the nematic field decay deep in the nematic phase, where it is proportional to the defect core radius. The third lengthscale arises by balancing the elastic stress associated with a deformation over a length with the active stress scale , to give the active length scale


To guarantee that any physics on these lengthscales is not contaminated by finite size effects, we focus on the regime in which and .

Alternatively, from a dynamical viewpoint one might consider the system to be controlled by two timescales: the passive structural relaxation time , which controls the relaxation of a distortion to the nematic order on a length scale , and the active time scale , which controls the relative rates of injection of active stresses and stress decay via viscous dissipation. The length scale that results when these timescales are equated is then


4.2 Emergent length scales

The length scales discussed above were motivated by simple dimensional analysis of the model parameters and flow geometry. In our numerical simulations, we find that (for a high enough level of activity) an initially homogeneous state gives way to a spatio-temporally complicated state with defects in the nematic director field, and associated local flows in the velocity field, as found earlier by several authors 11, 14 and shown in the snapshots of Figs. 2 and 3. An important aim of the present work is to elucidate how the length scales associated with these emergent structures depend on the underlying model length scales just discussed. We denote these emergent length scales by the common symbol , but in fact there are multiple possible scales that we might choose to characterize the spatio-temporal dynamics, as we now describe.

  • Mean defect separation : We define the mean defect separation


    where is the areal density of defects, calculated by adapting the defect tracking method of Ref. 31.

  • Director correlation length : The normalised director correlation function defined as


    This characterizes the probability that two director orientations a distance apart are the same (respecting the fact that are equivalent for a nematic). Here and throughout, the angular brackets indicate an average over space and time. We then choose to be the length at which , as in Fig. 1 (inset).

  • Velocity correlation length : Analogously, the velocity correlation function,


    defines the velocity correlation length according to .

  • Vorticity correlation length : Finally, we define the correlation function for the local vorticity, , as


    and define the vorticity correlation length by .

Fig.  1: Nematic correlation function defined in Eq. 12 obtained from Model II for an extensile system in the regime of spatio-temporally chaotic behavior for and activities in the range (red) to (blue). Inset: Unscaled data, demonstrating our definition . Main: the same data collapse onto a single curve when rescaled by the active length .

4.3 Scaling hypothesis

Simple dimensional analysis based on the model length scales discussed in Sec. 4.1 suggests that the length scales of Sec. 4.2 characterizing the emergent structures in the fluid (whether or ) should obey a simple scaling relation of the form


where is a general scaling function.

Previous simulation studies3, 32, 30, 17, 15 have shown that all characteristic length scales, denoted generically by , decrease with increasing activity . At low activity, typically just a few defects are seen in the simulation box, as in the snapshots in Figs. 2c and  3b. At higher activity one obtains a state of fully developed turbulence with a much higher density of defects (Figs. 2d and 3c). In this highly turbulent regime we expect the emergent length scale to become much smaller than, and therefore independent of, the system size . The above scaling form is then accordingly expected to reduce to


In our simulations all scaling law measurements are taken safely within this regime of fully developed turbulence, such that the emergent length scales are free of finite size effects. We also explicitly demonstrate that finite-size effects indeed return when is no longer small, as illustrated by the snapshot of Fig. 2c.

It is also worth noting that at extremely large activities the defect density could in principle become so large that the defect spacing approaches the microscopic length scale . In this regime we would expect to be unable to decrease further upon any additional increase in activity, and so to saturate. However our simulations do not reach this limit and the inequality is always respected.

4.4 Form of the scaling function

Having proposed the existence of a scaling function in Eqn. 16, we now consider possible specific forms for this functional dependence of on the model parameters. Conflicting scaling laws for have been proposed in the existing literature 15, 16, 17. While all of these studies agree that , there remains an apparent discrepancy over the scaling of with the activity.

Using Model II, Thampi et al. 15, 16 have proposed that , which would correspond to in Eqn. 16 having a square root dependence on its first argument. In contrast, using Model I, Giomi 17 suggested the relation , which would correspond to a linear dependence of on its first argument.

Fig.  2: Results from Model II for the nematic correlation length (empty symbols) and defect spacing (filled symbols) as functions of the dimensionless activity for an extensile nematic (). (a) Lengthscales vs for various values of the microscopic correlation length: (red circles), (green squares), and (blue triangles). The remaining parameter values are given in Table 2. At small activity we see saturation due to finite size effects. (b) The curves collapse when and are rescaled by . In both frames the black dashed lines show . In Fig. 2a we also mark the power law obtained by Thampi et al. as a purple dot-dashed line. (c,d) Representative snapshots of for (c) and (d) . We set in both snapshots. Defects of topological charge are identified by green dots (+) and red squares (-). For videos see supplementary material.

A possible origin of this discrepancy is the differing dimensionality of the order parameter between the two studies: while both have , Refs. 15, 16 had , whereas Ref. 17 had . This motivates us to compare numerical results for both within a single study. However our results below will rule out differences in as a source of discrepancy. Another potential reason could be that the two studies in fact explored different parameter regimes given the high dimensionality of the parameter space in these models. Therefore in order to ascertain the generality of these scaling laws, we systematically explore wide ranges for the three relevant adimensional parameters for both models. Our results will show that both forms suggested by the earlier studies can indeed apply, each in a different regime: one in the regime of fully developed active turbulence, the other when the system size plays a non-trivial role.

5 Results

We now present the results of our simulations. We focus on the regime of fully developed turbulence, corresponding to activity large enough to avoid finite system-size effects () and yet small enough to avoid saturation of the defect spacing at the microscopic length (). We systematically explore the functional dependence of the emergent correlation lengths defined in Section 4.2 on the model parameters. Specifically in Model I we vary the activity, , and viscosity ratio, , keeping all other parameters fixed to the values in Table 2. In Model II we vary the activity and the nematic persistence length , with all other parameters fixed to the values in Table 2. We will show that in the region of fully developed active turbulence all of the emergent length scales defined above scale with the active length , in both models. We will additionally demonstrate that a weaker exponent might be obtained in the regime of less well developed turbulence, where the typical size of the emergent structures is an appreciable fraction of the box size.

5.1 Correlation lengths

Fig.  3: Results from Model I for the nematic correlation length (empty symbols) and defect spacing (filled symbols) as functions of the dimensionless activity for an extensile nematic (). (a) Length scales vs for various values of the viscosity ratio: (red circles), (green squares), (blue triangles), and (magenta diamonds). The values of the other parameters are given in Table 2. The black dashed lines denote a slope of . (b,c) Representative snapshots of the alignment tensor for in (b) the low activity regime () with low defect density and (c) the high activity regime () with high defect density. The color scale represents the magnitude of the order parameter and the black lines denote the local orientation of the director field. Topological defects with charge are shown as green dots (+) and red squares (-).

In this section, we present our results for the correlation lengths defined in Sec. 4.2. Our main focus will be on an extensile nematic, corresponding to . We shall briefly discuss the contractile case at the end of this section.

Extensile active matter

Orientational correlations. We begin by considering correlations in the nematic order parameter . Figs. 2a and 3a shows the director correlation length and the defect spacing as obtained from Model II and Model I, respectively. For sufficiently large activity , we find that in both models both lengths obey a clear scaling law (black dashed lines). Note that the defect spacing correlation length is consistently larger than by a factor . This is to be expected as correlations at the halfway point between two defects () should be similar to those at .

At smaller activities (i.e., for ) the data obtained with Model II show a saturation in the power law (leftmost data points in Fig. 2a). This can be attributed to that fact that the length scale of nematic structure now spans an appreciable fraction of the system size, as seen in the snapshots of Fig. 2c. It is possible that fitting a power law in this saturation regime could result in a less negative exponent than the found in the regime of fully developed turbulence. Indeed we find that the scaling suggested by Thampi et al. (purple dashed dotted line in Fig. 2a) matches our data reasonably well in this regime.

The data in Fig. 2a also suggests that both and scale linearly with . We verify this scaling explicitly in Fig. 2b by plotting and against activity. The data for various values of collapse neatly onto a single curve, demonstrating a clear linear relation between both correlation lengths and .

The data obtained with Model I shown in Fig. 3a focus on large activities and verify that in this regime the scaling of both and with holds regardless of the model used. (They do not probe the saturation with system size seen at lower activities in Model II.) Data obtained for different values of the viscosity ratio can be collapsed when plotted as shown in Fig. 3b, suggesting , although the range of variation of the viscosity ratio is not sufficient to provide convincing evidence of scaling.

Fig.  4: Velocity (, filled symbols) and vorticity (, empty symbols) correlation lengths, normalized by for (a) an extensile () and (b) a contractile () system. We explore several values of the viscosity ratio: (red circles), (green squares), (blue triangles), and (magenta diamonds). Frame (c) shows the defect spacing (filled symbols) and the director correlation length (empty symbols) for a contractile () active nematic as a function of for the same set of values of . All lengths scale as . The black dashed lines represent a slope of .

Taken together, the data obtained with the two models tests the functional dependence of the two nematic correlation lengths with respect to activity and the nematic persistence length . Once free of the system size, we find that both obey . Consistent with this scaling, replotting in Fig. 1 (main) the full director correlation function as a function of the rescaled coordinate gives good data collapse. Additionally, the data obtained with Model I suggest a scaling , but a larger range of values would be needed to verify this. Next we demonstrate that the same scaling form is observed for correlations lengths associated with the velocity field .

Velocity and vorticity correlation lengths. Using data obtained with Model I, we explore the dependence of the velocity correlation length and the vorticity correlation length (as defined in Section 4.2) on activity. In light of the results of the previous section, we directly plot both these lengths against the rescaled activity (see Fig. 4a). As shown previously for the orientation correlation lengths, we again observe that both and scale as , with all data sets falling approximately on a single curve. We stress that this behavior is different from that reported in Ref. 15, where it was argued that does not depend on activity, while scale as .

Contractile active matter.

So far we have presented data for extensile systems, corresponding to . However many examples of contractile active matter are found in nature, e.g., suspensions of Chlamydomonas algae33, or cytoskeletal actomyosin networks 34. Therefore in order to further demonstrate the generality of our results, we now briefly consider the contractile case (). Since the linear instability of the homogeneous state requires requires 13, for contractile systems we use . Our data, shown in Figs. 4(b,c), support the idea that the defect spacing (), director correlation length (), velocity () and vorticity () correlation lengths are all controlled by a single active length scale . We caution, however, that the mapping of rod-like extensile () onto disc-like contractile () only holds at the linear instability level; the full non-linear dynamics may be subject to additional instabilities depending on the specific parameter values. While we do not expect that this would significantly change the scaling behaviour, we defer a full study of these effects to future work.

5.2 Kinetic energy and enstrophy

Fig.  5: Scaling of kinetic energy (frames (a) and (c)) and enstrophy (frames (b) and (d)) with activity for both Models I and II. The left figure displays the results obtained from Model 1 by varying the viscosity ratio as shown. The right figure displays the results obtained from Model II by varying the nematic correlation length . The inset of frame (d) shows the scaling collapse of the kinetic energy when plotting against activity. In frame (d), data is shown for two numerical resolutions: dashed lines for , and solid lines for .

The above scaling relations were obtained using the correlation functions defined in Section 4, which are normalised so that each function, e.g., , approaches unity as the separation distance . (See Fig. 1.) The normalization constants themselves, however, (i.e., the denominators in Eqs. 13 and 14) also provide useful information as they are directly proportional to the mean kinetic energy and enstrophy of the system, given by


where the angular brackets again denote an average over space and time. These quantities can be obtained experimentally, for instance by using particle image velocity (PIV) to quantify the flow fields of active liquids, as done by Dunkel et al. 35 in suspensions of extensile B. subtilis bacteria. We now use our earlier findings to motivate the expected scaling relation of these flow properties with activity, and then verify our predictions with further numerical data from both models.

Using simple dimensional analysis, the characteristic velocity of activity-induced shear flows associated with distortion of the local nematic order over a length scale can be obtained from the force balance condition (Eq. 1) as


Our results indicate that for sufficiently large values of activity the physics is controlled by a single active length scale , with . Using this in Eq. 19, we find and . The scaling of the vorticity can be estimated as , which gives an enstrophy .

This scaling is consistent with the findings of Ref. 17 in which the author examined the typical size of vortex structures in the regime of spatio-temporal chaotic dynamics using what we refer to here as Model I and found that both the vortex size and the defect spacing appear to scale with the active lengthscale . Further evidence for this scaling can be found in the experiments of Ref. 35, which found that where is the characteristic vortex size: assuming that , this implies that as we have argued above. Our proposed scaling is not, however, in agreement with the findings of Ref. 14, 16. In those studies it was found that and , a result that cannot seemingly be reconciled with the simple assumption that .

In order to appraise these conflicting scaling laws, we perform simulations with both Model I and Model II and measure the kinetic energy (, see Figs. 5a, b) and enstrophy (, Figs. 5c, d). The data from both models clearly obey our expected scaling laws and (black dashed lines). With Model I our choice of units means that increasing the rotational viscosity is equivalent to reducing the solvent viscosity . Our dimensional analysis in Eq. 19 suggests that increasing should increase the characteristic velocity. We indeed observe this trend in our data in Figs. 4b, although simulations over a larger range of would be required to determine the exact scaling. By the same analysis, we also expect that should be proportional to for fixed , since . Our data from Model II explores several values of , and plotting against activity indeed leads to a reasonable curve collapse (Fig. 5d inset). Consistent with the findings of Giomi17, we observe no appreciable dependence of on . This follows again from the scaling, .

6 Discussion

Using two distinct continuum models that have been studied extensively within the literature, we have performed a detailed numerical study of an active nematic to examine the scaling with activity of a number of structural and hydrodynamic correlation lengths, including the mean defect spacing. Our findings are consistent with the suggestion first put forward in Ref. 17 that in the regime of fully developed active turbulence defect proliferation, and the associated turbulent-like dynamics of the active nematic, are controlled by a single length scale . This is also the length scale that controls the onset of spontaneous flow instability of active films 2, 3, 36. Our numerical data from both models show that all measures of correlation length considered scale with this length scale, for both extensile and contractile systems.

Two caveats must, however, be applied. First, for extremely large activities (i.e., ) activity-induced deformations below the nematic persistence length are expected to be suppressed. Secondly, at low activities, structures can form that span the system size, and correlation lengths will correspondingly saturate, (i.e., ). We have explicitly demonstrated this system-size saturation in our simulations, a result that reconciles the apparently conflicting power law exponents previously reported in the literature.

Finally, to further support our findings, we have calculated the average kinetic energy and average enstrophy of the system, quantities that are readily obtainable from experiment. Our numerical results show that the scaling of these quantities with activity is consistent with a simple dimensional analysis based on the assumption that the physics is controlled by the single length scale .

Our results show that the key scaling relations hold for both strictly 2D and quasi-2D models. Encouragingly, this implies that such models capture the dynamics of active nematics in a generic way, i.e., independent of the specifics of the model. How our results would compare with the equivalent fully- simulation of an active nematic remains an interesting open question.

7 Acknowledgment

We thank Mark Bowick for introducing us to Ref. 31 (that describes a method for defect tracking) and Luca Giomi for writing the code used to study Model I. We thank them both and Mike Cates for invaluable discussions. MCM and PM were supported by the National Science Foundation through award DMR-1305184 and by the Syracuse Soft Matter Program. EJH thanks EPSRC for a Studentship. SMF’s and EJH’s research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement number 279365. The authors thank the KITP at the University of California, Santa Barbara, where they were supported through National Science Foundation Grant No. NSF PHY11-25925.

Appendix A Flow-aligning Parameter

We include here a comparison between the tumbling parameter used here and the Leslie-Erickson (LE) tumbling parameter (where corresponds to flow-aligning regime and corresponds to flow-tumbling regime). This comparison is presented in Appendix B of Ref. 3 for the case , but to our knowledge has not been displayed before for the case .

In dimensions the nematic tensor of a uniaxial nematic can be written as


where and are the sum and difference, respectively of the two principal values of the moment of inertia tensor of uniaxial nematogens. In our case we use that corresponds to needle-like molecules37.

We write the dynamical equation for the alignment tensor in dimension using the notation of Olmsted29,


where , , and are parameters that couple order and flow.

Substituting the expression given in Eq. 20 for the alignment tensor into the dynamical equation Eq. 21, and assuming to be constant, we obtain an equation for the director,


Comparing Eq. 22 to the Leslie-Erickson equation 29,


we identify the correspondence between the Olmsted coefficients and the Leslie-Erickson coefficients as


Using in Eq. 24, we obtain


Finally, for the case , we find


The case was previously reported in Ref. 3.


  1. footnotetext:  Department of Physics, Durham University, Science Laboratories, South Road, Durham, DH1 3LE, UK
  2. footnotetext:  Physics Department and Syracuse Soft Matter Program, Syracuse University, Syracuse, NY 13244, USA
  3. footnotetext: † Electronic Supplementary Information (ESI) available: videos of the two regimes of activity driven turbulence in Model II. See DOI: 10.1039/C6SM00812G.
  4. footnotetext: ‡  These authors contributed equally to this work.
  5. Note that the typical activity-induced velocity scale is in the range (see Fig. 5b or Eq. 19), meaning that the effective Reynolds number is in the range . These values are still small enough to ensure that any turbulence is activity driven (rather than inertial) in nature.
  6. In Model I the isotropic-nematic transition is continuous and the equilibrium nematic correlation length given by diverges at the transition. Deep in the nematic state where we can approximate .


  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. R. Voituriez, J. F. Joanny and J. Prost, Europhys. Lett., 2005, 70, 404–410.
  3. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E, 2007, 76, 031921.
  4. L. Giomi, M. C. Marchetti and T. B. Liverpool, Phys. Rev. Lett., 2008, 101, 198101.
  5. S. Ramaswamy, R. A. Simha and J. Toner, EPL (Europhysics Letters), 2003, 62, 196.
  6. S. Mishra and S. Ramaswamy, Phys. Rev. Lett., 2006, 97, 090602.
  7. V. Narayan, S. Ramaswamy and N. Menon, Science, 2007, 317, 105.
  8. L. Giomi, L. Mahadevan, B. Chakraborty and M. F. Hagan, Phys. Rev. Lett., 2011, 106, 218101.
  9. L. Giomi, L. Mahadevan, B. Chakraborty and M. F. Hagan, Nonlinearity, 2012, 25, 2245.
  10. L. M. Pismen, Phys. Rev. E, 2013, 88, 050502.
  11. L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, Phys. Rev. Lett., 2013, 110, 228101.
  12. Shi Xia-qing and Ma Yu-qiang, Nat Commun, 2013, 4, 1218.
  13. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. Cristina Marchetti, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2014, 372, 0365.
  14. S. P. Thampi, R. Golestanian and J. M. Yeomans, Phys. Rev. Lett., 2013, 111, 118101.
  15. S. P. Thampi, R. Golestanian and J. M. Yeomans, EPL (Europhysics Letters), 2014, 105, 18001.
  16. S. P. Thampi, R. Golestanian and J. M. Yeomans, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2014, 372, 0366.
  17. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
  18. S. P. Thampi, R. Golestanian and J. M. Yeomans, Molecular Physics, 2015, 113, 2656–2665.
  19. T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431–434.
  20. F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic and A. R. Bausch, Science, 2014, 345, 1135.
  21. G. Henkin, S. J. DeCamp, D. T. N. Chen, T. Sanchez and Z. Dogic, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2014, 372, 0142.
  22. S. Zhou, A. Sokolov, O. D. Lavrentovich and I. S. Aranson, Proc. Nat. Acad. Sci. U.S.A., 2013, 111, 1265–1270.
  23. G. Duclos, S. Garcia, Y. HG and P. Silberzan, Soft Matter, 2014, 10, 2346–2353.
  24. T. Gao, R. Blackwell, M. A. Glaser, M. D. Betterton and M. J. Shelley, Phys. Rev. Lett., 2015, 114, 048101.
  25. P. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford: Oxford University Press., 2nd edn., 1993.
  26. A. Sokolov, I. S. Aranson, J. O. Kessler and R. E. Goldstein, Phys. Rev. Lett., 2007, 98, 1–4.
  27. H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Lowen and J. M. Yeomans, Proc. Natl. Acad. Sci., 2012, 109, 14308–14313.
  28. M. Ravnik and J. M. Yeomans, Phys. Rev. Lett., 2013, 110, 026001.
  29. P. Olmsted, Rheologica Acta, 2008, 47, 283–300.
  30. E. J. Hemingway, A. Maitra, S. Banerjee, M. C. Marchetti, S. Ramaswamy, S. M. Fielding and M. E. Cates, Physical Review Letters, 2015, 114, 098302.
  31. D. Huterer and T. Vachaspati, Phys.Rev. D, 2005, 72, 043004.
  32. S. M. Fielding, D. Marenduzzo and M. E. Cates, Physical Review E, 2011, 83, 041910.
  33. S. Rafaï, L. Jibuti and P. Peyla, Phys. Rev. Lett., 2010, 104, 1–4.
  34. P. M. Bendix, G. H. Koenderink, D. Cuvelier, Z. Dogic, B. N. Koeleman, W. M. Brieher, C. M. Field, L. Mahadevan and D. A. Weitz, Biophys. J., 2008, 94, 3126–3136.
  35. J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 228102.
  36. S. A. Edwards and J. M. Yeomans, EPL (Europhysics Letters), 2009, 85, 18008.
  37. H. Stark and T. C. Lubensky, Phys. Rev. E, 2003, 67, 061709.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description