Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Active turbulence and spontaneous phase separation in inhomogeneous extensile active gels

Renato Assante , Dom Corbett *, Davide Marenduzzo and Alexander Morozov
SUPA, School of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK. E-mail: Dom.Corbett@ed.ac.uk; dmarendu@ph.ed.ac.uk

Received 2nd September 2022 , Accepted 30th November 2022

First published on 1st December 2022


Abstract

We report numerical results for the hydrodynamics of inhomogeneous lyotropic and extensile active nematic gels. By simulating the coupled Cahn–Hilliard, Navier–Stokes, and Beris–Edwards equation for the evolution of the composition, flow and orientational order of an active nematic, we ask whether composition variations are important to determine its emergent physics. As in active gels of uniform composition, we find that increasing either activity or nematic tendency (e.g., overall active matter concentration) triggers a transition between an isotropic passive phase and an active nematic one. We show that composition inhomogeneities are important in the latter phase, where we find three types of possible dynamical regimes. First, we observe regular patterns with defects and vortices: these exist close to the passive–active transition. Second, for larger activity, or deeper in the nematic phase, we find active turbulence, as in active gels of uniform composition, but with exceedingly large composition variation. In the third regime, which is uniquely associated with inhomogeneity and occurs for large nematic tendency and low activity, we observe spontaneous microphase separation into active and passive domains. The microphase separated regime is notable in view of the absence of an explicit demixing term in the underlying free energy which we use, and we provide a theoretical analysis based on the common tangent construction which explains its existence. We hope this regime can be probed experimentally in the future.


1 Introduction

Active gels1–3 are fascinating examples of non-equilibrium soft matter. Some well-known realisations of these systems include solutions of cytoskeletal filaments, such as actin or microtubules, interacting with molecular motors, such as myosin or kinesin.4,5 Other instances come from living materials, and encompass microbial suspensions of algae or bacteria.6,7 In an active gel, the constituent particles exert non-thermal forces on their environments. Such forces can be modelled, at the simplest level, as force dipoles, whose direction defines a nematic order parameter, which is a fundamental quantity to describe the emergent physics of these systems.8

The activity arising from the distribution of force dipoles leads to a phenomenology which is strikingly different from that of passive colloidal particles or polymer suspensions. For instance, these materials harbour a “spontaneous flow” instability, which sets in for sufficiently strong activity, and comprises a non-equilibrium transition between a quiescent suspension and a state where activity fuels continuous motion.9–13 For sufficiently large activity, the flow and orientation patterns of the spontaneously flowing state are seemingly chaotic, and the associated state is known as “active turbulence”.14–19 In active turbulence, active gels self-organise into a random arrangement of vortices. Experiments and theories suggest that in the nematic phase these vortices have a typical length scale, arising from the competition between activity and elasticity,14 while recent work points to important fundamental differences between active turbulence and its more widely studied passive counterpart.18,20 Active gels also possess strongly non-Newtonian rheological properties,21 such as marked activity-induced thinning or thickening,6,11,13,22,23 Darcy-like flow,24 or negative drag in microrheology.22,25

Existing theories and simulations of active gel hydrodynamics typically consider systems with constant composition. In contrast, inspection of active turbulent patterns found with microtubule–kinesin mixtures in the presence of polyethylene-glycol (which causes adsorption to the oil–water interface5,26) shows that the concentration of active material is significantly inhomogeneous. While a linear stability analysis shows that in extensile gels, such as a microtubule–kinesin mixture, the onset of spontaneous flow depends on orientational bend fluctuations and compositional fluctuations should be irrelevant,27 active turbulence is a highly non-linear phenomenon and the relevance of composition inhomogeneities to its physics remains unclear. Additionally, passive colloidal particles aggregate in active nematics,28 through a mechanism reminiscent of path coalescence29,30 or fluctuation-dominated phase ordering.31 This nonequilibrium aggregation shows that even a one-way coupling between composition and spontaneous flow (as tracers are affected by the spontaneous flow but do not modify it) can in principle give rise to composition inhomogeneities, and more in general to nontrivial physics outside the reach of a constant-composition approximation.

To understand the role of compositional inhomogeneities in the physics of spontaneous flow and active turbulence, here we study the hydrodynamic equations of motion of a mixture of an isotropic fluid and an active nematic gel by means of computer simulations. While the overall system is always incompressible, the active gel component can change its concentration, as is realistic for the microtubule–kinesin mixtures considered in ref. 26. With respect to previous work on Cahn–Hilliard models coupled to active nematics focussing on the crossover between wet and dry systems driven by friction with the substrate,32–34 here our focus is specifically on the qualitative role of compositional inhomogeneities on the emerging physics and patterns. As in active gels of uniform composition, we find that the system displays a transition between a passive isotropic phase and an active nematic phase. This transition can be triggered either by increasing activity or the nematic tendency of the system (more specifically the coupling between active matter concentration and orientational order). In the active nematic phase, though, compositional inhomogenities play a fundamental role and give rise to some unique dynamical behaviour. We find three dynamical regimes in this phase. Close to the transition boundary, our lyotropic system settles into regular flowing patterns, with approximately ordered spiral defect arrangements creating a rotational active flow consisting of long-lived and stable vortices. For larger activity, or deeper in the nematic phase, the flow becomes chaotic and we observe active turbulence. Both regular patterns and active turbulence are regimes which can be found in active gels of uniform composition. There are differences though, as in our case the thermodynamic coupling between orientational and compositional order parameters favours a concentration minimum, or relative void of active matter, at defect cores. Additionally, active turbulent patterns are characterised by a very broad distribution of local concentrations. The last regime we observe is unique to inhomogeneous systems, and is found even deeper in the nematic phase with respect to active turbulence, but for low activities. Here the active mixture spontaneously phase separates into low-concentration disordered and high-concentration nematic domains of irregular shape. We show by a semi-analytical theoretical analysis that this phase separation is due to the coupling between composition and nematic ordering and hence is driven thermodynamically. The corresponding coarsening is arrested in our case by the spontaneous active flow, and the size of the steady state domains decreases with activity. We conclude by discussing ways in which our study can be taken forward, both theoretically and experimentally.

2 Equations of motion

To describe the equilibrium properties of an inhomogeneous active nematic system in the passive phase (i.e., when the activity parameter defined below is switched off), we employ a Landau–de Gennes free energy [scr F, script letter F], whose density we call f. The latter consists of a contribution which characterises the orientational order of the active liquid crystal, measured by a nematic tensor Qαβ, and of another contribution determining the physics of compositional fluctuations, depending on the compositional order parameter ϕ (which measures the local concentration of active material). The liquid crystalline free energy density we use is a standard one to describe passive nematic liquid crystals,35 and is explicitly given by
 
image file: d2sm01188c-t1.tif(1)
where the first term is a bulk contribution, describing the isotropic–nematic transition, while the second term is an elastic distortion term. In the equation above, A0 is a constant, γ(ϕ) controls the magnitude of order (so that it may be viewed as an effective temperature or concentration for thermotropic and lyotropic liquid crystals respectively), while K is an elastic constant – note we are using the (standard in this field) one-constant approximation.35 Here and in what follows Greek indices denote Cartesian components and summation over repeated indices is implied. The coupling between concentration and ordering arises through γ(ϕ), which equals
 
γ(ϕ) = γ0 + ϕ(r,t)Δ,(2)
where γ0 and Δ are appropriate constants. In our simulations described below, we fix γ0 and vary Δ (see Section 3 for full parameter list).

The free energy density used to describe compositional fluctuations is instead given by a simple function, used to describe binary fluid in the mixed (non-phase-separating) regime, and its form is simply given by

 
image file: d2sm01188c-t2.tif(3)
where a is a constant related to the compressibility of the active gel component. Note we do not include a surface-tension-like square gradient term, image file: d2sm01188c-t3.tif, required for stabilisation in conventional binary mixture models, as density variations are already penalised thermodynamically by the term proportional to the elastic constant K in eqn (1). Note also that the total free energy density f = fLC + fϕ has two contributions which depend on ϕ: besides the mixing free energy fϕ, eqn (3), there is also a ϕ dependence in the liquid crystalline free energy fLC, eqn (1), through the γ(ϕ) term.

The fluid velocity, u, obeys the continuity equation and the Navier–Stokes equation,

 
αuα = 0,(4)
 
ρ(∂t + uββ)uα = −∂αp0 + ηβ2uα + ∂βΠαβζβ(ϕQαβ),(5)
where ρ is the fluid density, η is an isotropic viscosity, and
 
image file: d2sm01188c-t4.tif(6)
 
image file: d2sm01188c-t5.tif(7)
is the stress tensor, where ζ is the activity,8 and measures the strength of active force dipoles. With the sign convention chosen here ζ > 0 means extensile rods and ζ < 0 means contractile ones.8 The molecular field H which provides the driving motion is given by
 
image file: d2sm01188c-t6.tif(8)
where Tr denotes the tensorial trace.

The equation of motion for Q is taken to be36

 
tQ + u·∇Q = ΓH + S,(9)
where Γ is a collective rotational diffusion constant. The first term on the left-hand side of eqn (9) is the material derivative describing the usual time dependence of a quantity advected by a fluid with velocity u. This is generalized for rod-like molecules by a second term
 
image file: d2sm01188c-t7.tif(10)
where
 
image file: d2sm01188c-t8.tif(11)
are the symmetric part and the anti-symmetric part respectively of the velocity gradient tensor ∂βuα. The constant ξ depends on the molecular details of a given liquid crystal. The first term on the right-hand side of eqn (9) describes the relaxation of the order parameter towards the minimum of the free energy.

Finally, the active material concentration, ϕ, obeys a Cahn–Hilliard-like equation,

 
tϕ + u·∇ϕ = M2μ,(12)
where image file: d2sm01188c-t9.tif is the chemical potential of the active mixture, and M is a mobility, which for simplicity we consider to be constant.

3 Numerical method

To study the dynamics of eqn (5)–(12), here we perform direct numerical simulations. Note that in our simulations we assume that the fields are two-dimensional and that the Q tensor describes nematic order in a 2D plane (i.e., we assume that there is no out-of-plane nematic order). We employ an in-house MPI-parallel code developed within the Dedalus spectral framework.37 Simulations are performed on a periodic rectangular domain [0,H] × [0,H] with H = 200 (here and below, all quantities are given in simulation units). All fields are represented by de-aliased, double periodic Fourier series with 512 × 512 Fourier modes. Our time-iteration scheme employs a 4th-order semi-implicit backward differentiation formula (BDF) scheme38 with the timestep dt = 0.1. We have confirmed that our spatial and temporal accuracy is sufficient to obtain numerically converged results. To obtain statistically converged averages, simulations are performed for 30[thin space (1/6-em)]000 time units.

In what follows, we fix ρ = 2, η = 5/3, ξ = 0.7, A0 = 0.1, K = 0.01, γ0 = 2, Γ = 1, and a = 0.003, and M = 4. These parameters are chosen as they are in line with those used in previous hybrid lattice Boltzmann simulations of both constant-composition active nematics11,12 and self-motile active gel droplets.39

4 Results

To discuss our results, we first present numerical simulations, and then a semi-analytical discussion of phase separation in the passive limit which sheds light on the phase diagram which we find.

4.1 Simulation results

To address the role of compositional inhomogeneities in the hydrodynamics of extensile active gels, we first report numerical simulations (for methodology details, see Section 3). We characterise the dynamical behaviour of the system as a function of two key parameters: activity, quantified by ζ, and tendency to acquire nematic order, quantified by γ(ϕ0) = γ0 + ϕ0Δ. In practice, we change ζ and Δ in our simulations, keeping other parameters fixed (see Section 3 for a full list).

For each set of parameter values, we compute: (i) the average largest eigenvalue of Q, denoted by 〈q〉, which we use to quantify the global magnitude of order; (ii) the average fluid velocity, 〈|u|〉; and (iii) the average variance of the compositional order parameter. In each case, these averages are computed first spatially, over a configuration, and then over time. We acquire data when the system is in a statistical steady state, removing any initial transient. These quantities allow us to build a phase diagram.

The magnitude of order, 〈q〉, is shown in Fig. 1 and allows us to map the boundary between the isotropic and nematic phase. The nematic phase can be reached by increasing either Δ (thereby γ(ϕ0)) or ζ. The former is a thermodynamic route, the latter is a non-equilibrium one. The non-equilibrium transition between isotropic and nematic phases arises because activity-induced flows create shear, which in turn generates nematic order. Analogous transitions in systems with constant composition or polar order have been studied in ref. 32, 40 and 41. A linear stability calculation analogous to that of40 shows that an isotropic state of constant ϕ is unstable, for sufficiently large ζ, to an active nematic state with non-zero flow and order. The critical threshold in the (Δ, ζ) plane is described by ζc(Δ), where

 
image file: d2sm01188c-t10.tif(13)
corresponding to a straight line on the (Δ, ζ) plane. This predicted boundary is shown as a yellow line in Fig. 1. It agrees well with our numerics for small Δ, sufficiently far from the passive isotropic–nematic transition (at ζ = 0 and Δ = 0.7). As the latter is a first-order discontinuous transition, which requires the inclusion of non-linear terms in the equation of motion to be accurately described, it is unsurprising that there is a quantitative discrepancy with our linearised calculation close to this point.


image file: d2sm01188c-f1.tif
Fig. 1 Heatmap of the magnitude of nematic order as a function of ζ and Δ. The plot can be used to define regions in parameter space where the system is in the isotropic passive (bottom left) or active nematic (top right) phase. The prediction for the onset of spontaneous flow (passive–active transition) from linear stability analysis (see text) is shown as a yellow line. Note that the isotropic–nematic transition in the passive limit (ζ = 0) occurs at γ(ϕ0) = 2.7 (corresponding to Δ = 0.7; cf. Section 4.2).

Throughout the nematic phase, except at ζ = 0, we find a non-zero flow in steady state. In other words, the emerging nematic structures are always active in the parameter range which we have explored. To characterise their behaviour, we show in Fig. 2 the dynamical regimes we observe, before discussing the phase diagram quantitatively.


image file: d2sm01188c-f2.tif
Fig. 2 From top to bottom rows represent: snapshots of concentration field, {ϕ(x,y)|x,y∈ [0,200]}; director field, [n with combining circumflex] (scaled by local magnitude of order, q, and coarse grained for visibility), overlaid on concentration (for x ∈ [87,113], y ∈ [168,194]); vorticity (ωz); largest eigenvalue of Q tensor (q); |Qxy|, qualitatively corresponding to a Schlieren pattern as could be obtained experimentally with crossed polarisers. Parameters are as in section except: (a–e) ζ = 0.16, Δ = 0.3; (f–j) ζ = 0.2, Δ = 1; (k–o) ζ = 0.02, Δ = 1. Colour scales are shared between (f–j) and (k–o).

Fig. 2(a–e) shows an example of regular active patterns. These are found close to the isotropic–nematic transition for moderate values of Δ. The example shown features regular compositional modulations in steady state (Fig. 2a), where local minima are collocated with spiral defects in the director field of topological charge 1. The latter are most easily visible through the 4-brush pattern in the Schlieren-like plot of |Qxy| in Fig. 2e,35 and also correspond to deep minima in the local nematic order parameter, as shown in Fig. 2d. There are also steady vortices associated with the pattern, because spirals continuously rotate as in previous models of defects in constant-composition active nematics.42 Such vortices can be identified as maxima and minima in the vorticity plot in Fig. 2c. This collocation with vortices provides support for the interpretation that these thermodynamically unstable structures are stabilised in steady state by the active flow. It is the thermodynamic coupling (proportional to Δ) between order parameter and concentration in our free energy that drives the local concentration depletion at the centres of spirals, giving rise to a correlation between concentration and order, where larger concentration is associated with larger order, and vice versa. Thus the elastic energy cost associated with defect formation is decreased through proximal depletion, which explains the collocation of defects and voids. Our simulations show that the regular defect patterns we find close to the transition can either be stationary or self-motile (see Movie S2, ESI for the dynamics correspondent to the snapshot in Fig. 2a–e)

Fig. 2(f–j) shows an example of a different dynamical regime, obtained deeper in the nematic phase, for sufficiently large ζ and Δ (see also Movie S3, ESI). Here, the patterns are never stable but display a chaotic dynamics and diffuse around in the system. In line with constant-composition active nematic literature,15,18,20 we refer to these spatiotemporally varying patterns as active turbulence. Our simulations show that active turbulent patterns are accompanied by the appearance of defects – mostly of half-integer topological charge, corresponding to two-brush patterns in |Qxy| (Fig. 2j) – as in active gels of uniform composition. We also find this regime is characterised by large concentration variations (Fig. 2f). As in the case of regular patterns, here too the concentration field tends to decrease close to defect cores (Fig. 2f and g). As the active chaotic spontaneous flow moves defects around, streaks of concentration voids form (Fig. 2f and g) which follow defect trajectories.

While the regimes in Fig. 2(a–e) and (f–j) are qualitatively reminiscent of those found in the literature for constant-composition active extensile gels,11,32,40 we also find a third dynamical regime which is unique to inhomogeneous systems. This regime is found for low activity and sufficiently large Δ, and consists in spontaneous phase separation into active and passive domains (Fig. 2(k–o), see also Movie S1, ESI). Active domains are nematically ordered. Importantly, this is an example of microphase separation, or arrested phase separation, as coarsening does not proceed indefinitely and there are multiple domains in steady state (Fig. 2k, l and Fig. S1, ESI). In other words, the late-time domain size – computed, for instance, via the inverse first moment of the structure factor – does not scale with system size and decreases with increasing ζ (see ESI, Fig. S1 and S2). The presence of this spontaneous microphase separation regime is at first sight surprising, as the ϕ-dependent part of the free energy does not promote demixing by itself. As we show in Section 4.2, however, phase separation is driven by the coupling between nematic order and local concentration of active matter.

To quantitatively delineate the phase diagram of the system, and the boundaries between the three different dynamical regimes shown in Fig. 2, we proceed as follows. We use the plot in Fig. 1 to identify the phases as isotropic (〈q〉 ≃ 0) or nematic (〈q〉 ≠ 0). We then classify cases where the kinetic energy reaches a plateau as regular patterns. To discriminate between active turbulence and spontaneous phase separation, we look at probability distribution functions for ϕ (calculated over space, and averaged over time, see Fig. 3 for examples). If this distribution is bimodal (i.e., it has two maxima), then we classify the pattern as phase separated (Fig. 3, yellow curve). A similar identification would arise from analysing the average concentration variance and the average flow magnitude (Fig. 4). Flow and concentration variance are highest in the active pattern and spontaneous phase separation regimes respectively. The phase diagram for our inhomogeneous extensile mixture is shown in Fig. 5.


image file: d2sm01188c-f3.tif
Fig. 3 Probability distribution functions showing frequency versus local concentration at a grid point, for the active pattern (blue curve), the active turbulence (maroon curve) and the spontaneous phase separation (yellow curve) regimes. It can be seen that the distribution corresponding to the active turbulent regime is unimodal, while that associated with the spontaneous phase separated one is bimodal, with peaks at ϕ ≈ 0.19, 1.41, cf. binodal points in Fig. 6.

image file: d2sm01188c-f4.tif
Fig. 4 Combined heatmap of the averaged flow magnitude, |u| (blue), and variance of ϕ (yellow), as functions of ζ and Δ. The plot can be used to identify the regions of parameter space where phase separation is most prominent (upper left) or where the steady active flow is strongest (right).

image file: d2sm01188c-f5.tif
Fig. 5 Phase diagram of the system in the (ζ, Δ) plane. Each square corresponds to one simulation. Regions (A–D) correspond respectively to isotropic passive, active patterns, active turbulence, and spontaneous microphase separation regimes.

The concentration distribution plots show that concentration variations are a generic feature throughout the active nematic phase. First we note that, without such concentration variations, the spontaneous microphase separation regime could not arise. The existence of this regime in a realisable system depends on the value of Δ, which will be determined by microscopic parameters. Assuming an Onsager-like theory for composition-dependent orientational order, we expect this could be achievable with rod-like active particles of sufficiently large aspect ratio. Note that we actually observe a bimodal distribution even in the regular pattern regime, however here the magnitude of the peaks differs by more than one order of magnitude. This corresponds to a detectable depletion of active material at the cores of long-lived defects rather than actual phase separation. A second observation is that the width of the concentration distributions is remarkably large in the active turbulent regime, as the sampled values of ϕ can vary from ∼0 to ∼2. This range is similar to the spread observed in the spontaneous phase separation regime, although the distribution remains peaked at ϕϕ0 = 1 for active turbulent patterns. This observation is in line with experimental microscopy of quasi-2D extensile mixtures of microtubules and molecular motors, which shows a highly inhomogeneous concentration of active material.26

4.2 Coupling-induced demixing in the passive limit

To shed more light on the physics underlying the occurrence of our spontaneous microphase separated regime, it is instructive to consider the passive limit of our nematic mixture, where the behaviour is solely determined by thermodynamic considerations. In this case, the free energy is minimised in equilibrium, and a full phase diagram can be computed semi-analytically using the common tangent construction (Fig. 6).
image file: d2sm01188c-f6.tif
Fig. 6 Coupling-induced demixing derived semi-analytically from the free energy for a uniformly aligned passive nematic phase given in eqn (14). Panel (a) shows the free energy, fhom(qmin(ϕ), ϕ), which, for a given value of ϕ, is minimized by qmin(ϕ) (see Section 4.2). The second derivative with respect to ϕ is shown with a dashed curve, and is negative in the spinodal region, highlighted in yellow. Here the system is linearly unstable to demixing into regions of high (ϕ+) and low (ϕ) concentration, the values of which are given by the respective binodal points, shown in grey. These points are found by the common tangent construction illustrated by the green line, and share the same chemical potential and pressure (Π, dashed in blue). In this plot Δ is fixed at 0.7, and the presence of the lower spinodal point at ϕ = 1 compares well with our simulation results, as illustrated in Fig. 5. Panel (b) shows the ϕΔ plane, which is a phase diagram derived from the loci of the spinodal and binodal points illustrated in panel (a), considered as we vary ϕ and Δ. The shaded spinodal and binodal regions are bounded by these loci, and show the parameter combinations for which the system is unstable to phase decomposition. Note that all ϕ-axis values signify the spatial average over the system, which is dropped from the notation for readability.

Recall that our model employs both conserved (composition ϕ) and non-conserved (nematic tensor Qαβ) order parameters, which are mutually coupled through the free energy given in eqn (1). The resulting mixture is an example of a lyotropic liquid crystal, such as that considered in ref. 43, albeit with an important distinction. Specifically, in the limit of no coupling (Δ = 0 in eqn (2)), the system we consider exhibits no phase separation and mixes freely. In what follows, therefore, we shall see that passive phase separation is driven only by the coupling.

To find the equilibrium state, we consider a uniformly aligned (homogeneous) nematic phase, so that the elastic term contribution in eqn (1) vanishes and the total free energy can be written in terms of the magnitude of nematic order, q, and of ϕ as follows,

 
image file: d2sm01188c-t11.tif(14)
This homogeneous state has everywhere the same value of the conserved composition. Thus, for each point (x,y) in the system, image file: d2sm01188c-t12.tif. The magnitude of order that minimizes the total free energy for this given homogeneous value of ϕ is
 
image file: d2sm01188c-t13.tif(15)
If the homogeneous state is stable against phase decomposition into regions of high and low concentration, then there is no configuration for which 〈f〉, the spatial average of the free energy density, is less than the value of fhom(qmin(ϕ),ϕ). However, as the common tangent construction illustrated in Fig. 6a shows, for each homogeneous configuration lying between the binodal points, there exists a corresponding inhomogeneous configuration with the same average composition, 〈ϕ〉, but a lower average free energy density, 〈fdemixed〉, the value of which is given by the tangent line itself. Moreover, for sufficiently large Δ and certain values of ϕ, we find that fhom(qmin(ϕ),ϕ) is a non-convex function. In this spinodal region, we therefore expect that the nematic mixture is linearly unstable to phase separation. By comparison, in the limit of vanishing activity, our simulation found phase decomposition in the parameter range 〈ϕ〉 = 1, Δ > 0.7 (see, for instance, the left border of Fig. 5). This simulated demixing is highlighted in Fig. 6b, and signifies good agreement with our semi-analytic thermodynamic prediction.

In greater detail, the spinodal region – where the system is linearly unstable to phase separation, for any perturbation however small – is precisely where fhom(qmin(ϕ),ϕ) is concave down. Thus it is bounded by the inflexion points, where the second derivative of fhom with respect to ϕ, fϕϕ(qmin(ϕ),ϕ), crosses the ϕ axis in Fig. 6a. The lower spinodal point is simply ϕc = (γcγ0)/Δ, whereas the upper spinodal point, ϕs, was found numerically as the only real root of a seventh order polynomial.

This effective free energy also gives, by a common tangent construction, the binodal points, namely the values of the coexisting concentrations of the phase separated systems. Specifically, as coexisting phases must share a common chemical potential (μ = fϕ(qmin(ϕ),ϕ)) and a common pressure (Π = fhom(qmin(ϕ),ϕ) − μϕ), we have a system of two equations in two unknowns. The two unknowns in question are the two concentrations, ϕ and ϕ+, for which a straight line tangentially touches the curve fhom(qmin(ϕ),ϕ) exactly twice, and were found numerically. For values of 〈ϕ〉 between the binodal points, this tangent represents the average free energy density of the phase separated system, which is visibly lower than the free energy density of the homogeneous system. Thus 〈fdemixed〉 < fhom(qmin(ϕ),ϕ) for ϕ < 〈ϕ〉 < ϕ+.

The loci of these spinodal and binodal points as we vary ϕ and Δ are the curves plotted in the ϕΔ plane in Fig. 6b. Here we can see that, in line with our simulation results shown in Fig. 5, there is a spinodal point at Δ = 0.7 for 〈ϕ〉 = 1.

We stress again that the analysis in this section refers to the ζ = 0 passive limit of our mixture. In the active case, macroscopic phase separation is generically arrested by the activity-induced spontaneous flow. This leads to microphase separation at steady state, similarly to what has been observed in sheared passive binary mixtures44–46 or in active model H.39

Quantitatively, our theory predicts a spinodal point at Δ = 0.7 for 〈ϕ〉 = 1, which compares well with our simulation results (Fig. 5). The predicted binodal points for Δ = 1, lie at ϕ = 0.18 and ϕ = 1.79 – these values should be compared with the histogram peaks shown in Fig. 3, viz. ϕ = 0.19 and ϕ = 1.41 – the lower value of the latter peak is likely due to the active flow which drives the system away from thermodynamic equilibrium (assumed in Fig. 6). Note our simulations only find phase separation in the spinodal region in Fig. 6b, as the initial perturbation of the homogeneous state are smaller than the scale required for nucleation.

5 Discussion and conclusions

In summary, we have used computer simulations to study the hydrodynamics of an inhomogeneous active nematic gel. With respect to conventional models for active gels, which only consider the velocity field and Q tensor, our theory also allows for the time evolution of the active matter concentration ϕ. Previous work has shown by a linear stability analysis that compositional fluctuations are irrelevant for the physics of the “generic instability” of active gels,27 which stands for the transition between the passive (quiescent) and the active (spontaneously flowing) phase. It has however remained unclear what their role is deep in the active phase, where nonlinearities are important; shedding light on this issue has been the focus of our current work.

Our main result is the quantitative characterisation of the phase diagram of inhomogeneous active nematic (Fig. 5). We have found that there are three regimes with distinct emergent behaviour in the active phase. First, close to the transition between the passive isotropic and active nematic phase, there are regular patterns typically composed of self-assembled rotating spirals. Second, deeper in the active nematic phase there is an active turbulent regime featuring chaotic dynamics of vortices and half-integer nematic defects. Third, for low activity and large enough nematic tendency (Δ in our phase diagram in Fig. 5), we find spontaneous phase separation into active and passive domains. This latter phase separation is arrested by the active flow, so that domains do not coarsen past a typical size, which decreases with increasing activity.

The regular spiral/vortex patterns we find are reminiscent of those observed with polar active gels in the ordered phase.42,47 While polar nematics can only exhibit defects with integer topological charge, defects in active (apolar) nematics normally have half-integer topological charge, so it is non-trivial that close to the passive–active transitions we observe spirals, whose topological charge is +1. Notwithstanding this qualitative resemblance, the patterns we observe also have a non-trivial spatiotemporal dynamics (Movie S2, ESI). The core of our spirals are also associated with notable concentration minima, or voids, which arise because of the coupling between nematic order and composition in the free energy of the system. The mechanism responsible for this coupling is the same that drives inert colloidal particles or isotropic droplets (with no anchoring on their surface) to the defect cores or disclinations in passive liquid crystals.28,49

The chaotic, active turbulent regime we find for sufficiently large ζ is an analog of the regime of the same name in active gels of uniform composition.18,20,50,51 An important feature of this regime in our simulations, though, is that there are very large compositional fluctuations (Fig. 3a). These are qualitatively in line with experimental observations of active turbulence in microtubule–motor mixtures, which show substantial inhomogeneities in microtubule concentration over a sample.5,26 Whether such concentration variations lead to a fundamental change in the scaling properties of active turbulence is an open question which we believe deserves further investigation, for instance by a quantitative detailed analysis of the scaling of velocity–velocity correlations.16,17,52

Regarding the spontaneous microphase separated regime, this is notable especially because there is no term in the free energy density fϕ which explicitly favours demixing. In other words, in the absence of fLC the system would remain uniform. Phase separation arises due to the coupling between composition and order, measured by the parameter Δ. In this sense, phase separation is not driven by activity but rather thermodynamically, and indeed it can be explained with a theoretical discussion of the free energy in the passive limit (Fig. 6). In simulations we observe a microphase separated pattern rather than macroscopic phase separation, as the active flow arrests coarsening, and controls the size of the steady-state domains observed at late times, similarly to the case of active model H.39 While experimental realisation of active nematics have shown plenty of instances of active turbulence,5,26 the spontaneous microphase separation regime appears to not have been found in the lab yet. Our model suggests that the most promising avenue to realise this regime experimentally is to control the composition-order coupling Δ. The latter may be estimated by monitoring how the isotropic–nematic transition point depends on the concentration of nematogenic particles (for instance, microtubules) in the passive limit of no activity.

Looking ahead, we can suggest a few directions in which our work can be carried forward. First, it would be of interest to understand from a more fundamental point of view the universal properties of the dynamical regimes we have identified. For instance, one could quantify the dependence of vortex correlation length and pattern size on physical parameters, and the scaling of the power spectra of the kinetic energy. This would allow to clarify the important theoretical question of whether or not inhomogeneous active turbulence is in the same universality class as turbulence in active gels of uniform composition. Second, from the experimental point of view, it would be desirable to compare more quantitatively concentration distributions in active turbulence with those predicted by our simulations. Third, with regards to computer simulations, it would be highly interesting to explore the phase behaviour and dynamics of inhomogeneous active nematics in 3D, comparing and contrasting it with the one found here in 2D.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank S. Samatas for useful discussions. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission. This work was partially funded by EPSRC (grant number EP/V048198/1). We would also like to thank EPSRC for the computational time made available on the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk).

References

  1. S. Ramaswamy, Annu. Rev. Condens. Matter Phys., 2010, 1, 323–345 CrossRef.
  2. J. Prost, F. Jülicher and J.-F. Joanny, Nat. Phys., 2015, 11, 111–117 Search PubMed.
  3. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 1–13 CAS.
  4. F. C. MacKintosh and C. F. Schmidt, Curr. Opin. Cell Biol., 2010, 22, 29–35 CrossRef PubMed.
  5. T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431–434 CrossRef CAS PubMed.
  6. Y. Hatwalne, S. Ramaswamy, M. Rao and R. A. Simha, Phys. Rev. Lett., 2004, 92, 118101 CrossRef PubMed.
  7. H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 268102 CrossRef PubMed.
  8. 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 CrossRef CAS.
  9. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101 CrossRef PubMed.
  10. R. Voituriez, J.-F. Joanny and J. Prost, EPL, 2005, 70, 404 CrossRef CAS.
  11. D. Marenduzzo, E. Orlandini and J. Yeomans, Phys. Rev. Lett., 2007, 98, 118102 CrossRef CAS PubMed.
  12. D. Marenduzzo, E. Orlandini, M. Cates and J. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921 CrossRef CAS PubMed.
  13. M. Cates, S. Fielding, D. Marenduzzo, E. Orlandini and J. Yeomans, Phys. Rev. Lett., 2008, 101, 068102 CrossRef CAS PubMed.
  14. L. Giomi, Phys. Rev. X, 2015, 5, 031003 Search PubMed.
  15. A. Doostmohammadi, T. N. Shendruk, K. Thijssen and J. M. Yeomans, Nat. Commun., 2017, 8, 1–7 CrossRef PubMed.
  16. J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo and A. Morozov, Phys. Rev. Lett., 2017, 119, 028005 CrossRef PubMed.
  17. D. Bárdfalvy, H. Nordanger, C. Nardini, A. Morozov and J. Stenhammar, Soft Matter, 2019, 15, 7747–7756 RSC.
  18. L. N. Carenza, L. Biferale and G. Gonnella, EPL, 2020, 132, 44003 CrossRef CAS.
  19. T. Kozhukhov and T. N. Shendruk, Sci. Adv., 2022, 8, eabo5788 CrossRef PubMed.
  20. R. Alert, J.-F. Joanny and J. Casademunt, Nat. Phys., 2020, 16, 682–688 Search PubMed.
  21. D. Saintillan, Annu. Rev. Fluid Mech., 2018, 50, 563–592 CrossRef.
  22. G. Foffano, J. Lintuvuori, A. Morozov, K. Stratford, M. Cates and D. Marenduzzo, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 1–11 CrossRef PubMed.
  23. V. A. Martinez, E. Clément, J. Arlt, C. Douarche, A. Dawson, J. Schwarz-Linek, A. K. Creppy, V. Škultéty, A. N. Morozov and H. Auradou, et al. , Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 2326–2331 CrossRef CAS PubMed.
  24. F. Mackay, J. Toner, A. Morozov and D. Marenduzzo, Phys. Rev. Lett., 2020, 124, 187801 CrossRef CAS PubMed.
  25. G. Foffano, J. Lintuvuori, K. Stratford, M. Cates and D. Marenduzzo, Phys. Rev. Lett., 2012, 109, 028103 CrossRef CAS PubMed.
  26. B. Martnez-Prat, J. Ignés-Mullol, J. Casademunt and F. Sagués, Nat. Phys., 2019, 15, 362–366 Search PubMed.
  27. A. Baskaran and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15567–15572 CrossRef CAS PubMed.
  28. G. Foffano, J. S. Lintuvuori, K. Stratford, M. Cates and D. Marenduzzo, Soft Matter, 2019, 15, 6896–6902 RSC.
  29. A. Maritan and J. R. Banavar, Phys. Rev. Lett., 1994, 72, 1451 CrossRef PubMed.
  30. M. Wilkinson and B. Mehlig, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 040101 CrossRef CAS PubMed.
  31. D. Das and M. Barma, Phys. Rev. Lett., 2000, 85, 1602 CrossRef CAS PubMed.
  32. S. P. Thampi, A. Doostmohammadi, R. Golestanian and J. M. Yeomans, EPL, 2015, 112, 28004 CrossRef.
  33. A. Doostmohammadi, M. F. Adamer, S. P. Thampi and J. M. Yeomans, Nat. Commun., 2016, 7, 1–9 Search PubMed.
  34. K. Thijssen, D. A. Khaladj, S. A. Aghvami, M. A. Gharbi, S. Fraden, J. M. Yeomans, L. S. Hirst and T. N. Shendruk, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2106038118 CrossRef CAS PubMed.
  35. P.-G. De Gennes and J. Prost, The physics of liquid crystals, Oxford University Press, 1993 Search PubMed.
  36. A. N. Beris, B. J. Edwardset al., Thermodynamics of flowing systems: with internal microstructure, Oxford University Press, 1994 Search PubMed.
  37. K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet and B. P. Brown, Phys. Rev. Res., 2020, 2, 023068 CrossRef CAS.
  38. D. Wang and S. J. Ruuth, J. Comput. Math., 2008, 26, 838–855 Search PubMed.
  39. A. Tiribocchi, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 115, 188302 CrossRef PubMed.
  40. S. Santhosh, M. R. Nejad, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, J. Stat. Phys., 2020, 1–11 Search PubMed.
  41. M. G. Giordano, F. Bonelli, L. N. Carenza, G. Gonnella and G. Negro, EPL, 2021, 133, 58004 CrossRef CAS.
  42. K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, Phys. Rev. Lett., 2004, 92, 078101 CrossRef CAS PubMed.
  43. A. Matsuyama, R. Evans and M. Cates, Eur. Phys. J. E: Soft Matter Biol. Phys., 2002, 9, 79–87 CrossRef CAS PubMed.
  44. P. Stansell, K. Stratford, J.-C. Desplat, R. Adhikari and M. Cates, Phys. Rev. Lett., 2006, 96, 085701 CrossRef CAS PubMed.
  45. T. Imaeda, A. Furukawa and A. Onuki, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 051503 CrossRef PubMed.
  46. T. Taniguchi and A. Onuki, Phys. Rev. Lett., 1996, 77, 4910 CrossRef CAS PubMed.
  47. J. Elgeti, M. Cates and D. Marenduzzo, Soft Matter, 2011, 7, 3177–3185 RSC.
  48. M. L. Blow, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2014, 113, 248303 CrossRef PubMed.
  49. M. Ravnik, G. P. Alexander, J. M. Yeomans and S. Žumer, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 5188–5192 CrossRef CAS PubMed.
  50. S. P. Thampi, R. Golestanian and J. M. Yeomans, Phys. Rev. Lett., 2013, 111, 118101 CrossRef PubMed.
  51. M. Linkmann, G. Boffetta, M. C. Marchetti and B. Eckhardt, Phys. Rev. Lett., 2019, 122, 214503 CrossRef CAS PubMed.
  52. V. Škultéty, C. Nardini, J. Stenhammar, D. Marenduzzo and A. Morozov, Phys. Rev. X, 2020, 10, 031059 Search PubMed.

Footnotes

Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2sm01188c
Note that in our geometry integration of the topological charge density, defined as in ref. 48, is conserved and equal to 0.

This journal is © The Royal Society of Chemistry 2023
Click here to see how this site uses Cookies. View our privacy policy here.