Active turbulence and spontaneous phase separation in inhomogeneous extensile active gels

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 gels^1–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.⁸

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,¹⁴ 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,²¹ such as marked activity-induced thinning or thickening,^{6,11,13,22,23} Darcy-like flow,²⁴ 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 interface^5,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,²⁷ 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,²⁸ through a mechanism reminiscent of path coalescence^29,30 or fluctuation-dominated phase ordering.³¹ 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,³⁵ and is explicitly given bywhere 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, A₀ 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.³⁵ 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

γ(ϕ) = γ₀ + ϕ(r,t)Δ, (2)

where γ₀ and Δ are appropriate constants. In our simulations described below, we fix γ₀ 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

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, , 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 = f_LC + 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 f_LC, eqn (1), through the γ(ϕ) term.

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

∂_αu_α = 0, (4)

ρ(∂_t + u_β∂_β)u_α = −∂_αp₀ + η∂_β²u_α + ∂_βΠ_αβ − ζ∂_β(ϕQ_αβ), (5)

where ρ is the fluid density, η is an isotropic viscosity, andis the stress tensor, where ζ is the activity,⁸ and measures the strength of active force dipoles. With the sign convention chosen here ζ > 0 means extensile rods and ζ < 0 means contractile ones.⁸ The molecular field H which provides the driving motion is given bywhere Tr denotes the tensorial trace.

The equation of motion for Q is taken to be³⁶

∂_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 termwhereare 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·∇ϕ = M∇²μ, (12)

where 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.³⁷ 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) scheme³⁸ 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 000 time units.

In what follows, we fix ρ = 2, η = 5/3, ξ = 0.7, A₀ = 0.1, K = 0.01, γ₀ = 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 nematics^11,12 and self-motile active gel droplets.³⁹

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 γ(ϕ₀) = γ₀ + ϕ₀Δ. 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 γ(ϕ₀)) 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 of⁴⁰ 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

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.

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 γ(ϕ₀) = 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.

Fig. 2 From top to bottom rows represent: snapshots of concentration field, {ϕ(x,y)|x,y∈ [0,200]}; director field, n̂ (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); |Q_xy|, 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 |Q_xy| in Fig. 2e,³⁵ 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.⁴² 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 |Q_xy| (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.

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.

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).

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 ϕ ≃ ϕ₀ = 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.²⁶

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).

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, f_hom(q_min(ϕ), ϕ), which, for a given value of ϕ, is minimized by q_min(ϕ) (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,

This homogeneous state has everywhere the same value of the conserved composition. Thus, for each point (x,y) in the system, . The magnitude of order that minimizes the total free energy for this given homogeneous value of ϕ isIf 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 f_hom(q_min(ϕ),ϕ). 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, 〈f_demixed〉, the value of which is given by the tangent line itself. Moreover, for sufficiently large Δ and certain values of ϕ, we find that f_hom(q_min(ϕ),ϕ) 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 f_hom(q_min(ϕ),ϕ) is concave down. Thus it is bounded by the inflexion points, where the second derivative of f_hom with respect to ϕ, f_ϕϕ(q_min(ϕ),ϕ), crosses the ϕ axis in Fig. 6a. The lower spinodal point is simply ϕ_c = (γ_c − γ₀)/Δ, 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_ϕ(q_min(ϕ),ϕ)) and a common pressure (Π = f_hom(q_min(ϕ),ϕ) − μϕ), 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 f_hom(q_min(ϕ),ϕ) 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 〈f_demixed〉 < f_hom(q_min(ϕ),ϕ) 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 mixtures^44–46 or in active model H.³⁹

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

---

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.

---