Active transport in a channel: stabilisation by flow or thermodynamics

TL;DR

This study uses continuum simulations to explore how active fluids behave in channels, revealing stable flow states and key parameters influencing flow, which aids the development of active microfluidic devices.

Contribution

It provides a systematic analysis of active fluid behavior in channels, including phase diagrams and the influence of flow alignment or tumbling.

Findings

Identification of multiple stable flow states

Development of a phase diagram for active flows

Key parameters include channel width and flow alignment behavior

Abstract

Recent experiments on active materials, such as dense bacterial suspensions and microtubule-kinesin motor mixtures, show a promising potential for achieving self-sustained flows. However, to develop active microfluidics it is necessary to understand the behaviour of active systems confined to channels. Therefore here we use continuum simulations to investigate the behaviour of active fluids in a two-dimensional channel. Motivated by the fact that most experimental systems show no ordering in the absence of activity, we concentrate on temperatures where there is no nematic order in the passive system, so that any nematic order is induced by the active flow. We systematically analyze the results, identify several different stable flow states, provide a phase diagram and show that the key parameters controlling the flow are the ratio of channel width to the length scale of active flow vortices, and whether the system is flow aligning or flow tumbling.

Figures

Tables

Table 1. Table 1: Parameters used in the simulations

Parameter	Range
$\lambda$	0.3 - 1.0
A	-0.05 - 0.015
B	0 - 0.3
C	0 - 0.3
$\zeta$	0.00015 - 0.07
$\Gamma$	0.34
$\mu$	0.667
$\rho$	1
K	0.01
$L$	100 - 400
$W$	20 - 70

Equations

$$Q_{ij}=\frac{q}{2}(3n_i{n}_j-{\delta }_{ij})$$

$$F_{\text{LDG}}=\frac{A}{2}{Q}_{ij}{Q}_{ji}+\frac{B}{3}{Q}_{ij}{Q}_{jk}{Q}_{ki}+\frac{C}{4}{\left(Q_{ij}{Q}_{ji}\right)}^2,$$

$$F_{\text{FO}}=\frac{K}{2}{\left({\partial }_k{Q}_{ij}\right)}^2$$

$${\partial }_i{u}_i=0,$$

$$\rho \left({\partial }_t+u_k{\partial }_k\right)u_i={\partial }_j{\pi }_{ij}$$

$${\pi }_{ij}={\pi }_{ij}^{\text{viscous}}+{\pi }_{ij}^{\text{active}}+{\pi }_{ij}^{\text{passive}}.$$

$${\pi }_{ij}^{\text{viscous}}=2\mu E_{ij}$$

$${\pi }_{ij}^{\text{active}}=-\zeta Q_{ij}$$

$$\begin{array}{r}{\pi }_{ij}^{\text{passive}}=-P{\delta }_{ij}+2\lambda \left(Q_{ij}+\frac{{\delta }_{ij}}{3}\right)\left(Q_{kl}{H}_{lk}\right)-\lambda H_{ik}\left(Q_{kj}+\frac{{\delta }_{kj}}{3}\right)\\ -\lambda \left(Q_{ik}+\frac{{\delta }_{ik}}{3}\right)H_{kj}-{\partial }_i{Q}_{kl}\left(\frac{\delta F}{\delta {\partial }_j{Q}_{lk}}\right)+Q_{ik}{H}_{kj}-H_{ik}{Q}_{kj}\end{array}$$

$$H_{ij}=-\frac{\partial F}{\partial Q_{ij}}+\frac{{\delta }_{ij}}{3}\text{Tr}\left(\frac{\partial F}{\partial Q_{kl}}\right).$$

$$\left({\partial }_t+u_k{\partial }_k\right)Q_{ij}-S_{ij}=\Gamma H_{ij}$$

$$\begin{array}{r}{S}_{ij}=\left(\lambda E_{ik}+{\Omega }_{ik}\right)\left(Q_{kj}+\frac{{\delta }_{kj}}{3}\right)+\left(Q_{ik}+\frac{{\delta }_{ik}}{3}\right)\left(\lambda E_{kj}-{\Omega }_{kj}\right)\\ -2\lambda \left(Q_{ij}+\frac{{\delta }_{ij}}{3}\right)\left(Q_{kl}{\partial }_k{u}_l\right)\end{array}$$

$$\left({\partial }_t+u_k{\partial }_k\right)Q_{ij}-{\Omega }_{ik}{Q}_{kj}+Q_{ik}{\Omega }_{kj}=\Gamma H_{ij}^{\text{active}}$$

$$\begin{array}{r}{H}_{ij}^{\text{active}}=K{\nabla }^2{Q}_{ij}+\frac{\lambda \zeta }{3\mu }{Q}_{ij}+\frac{\lambda \zeta }{\mu }\left[Q_{ik}{Q}_{kj}-Q_{pq}{Q}_{qp}\frac{{\delta }_{ij}}{3}\right]\\ -\frac{\lambda \zeta }{\mu }{Q}_{ij}{Q}_{pq}{Q}_{qp}\end{array}$$

$${\varphi }_1={⟨\left|{⟨\frac{u_x(x,y)}{|\mathbf{u}(x,y)|}⟩}_x\right|⟩}_y,$$

$${\varphi }_2={⟨\left|{⟨\frac{u_y(x,y)}{|\mathbf{u}(x,y)|}⟩}_y\right|⟩}_x.$$

$${\lambda }_1=\left(\frac{3q+4}{9q}\right)\lambda .$$

Full text

I Introduction

Collective behaviour is present in biological systems at all length scales. Examples include cytoskeletal elements Kruse et al. (2004); Schaller et al. (2010), cell colonies Saw et al. (2017); Duclos et al. (2018), bacterial suspensions Wensink et al. (2012); Dunkel et al. (2013), locust swarms Buhl et al. (2006) and fish schools Lopez et al. (2012). Similar behaviour is seen in synthetic systems, such as catalytic colloids Harder and Cacciuto (2018); Buttinoni et al. (2013) and vibrated granular rods Narayan et al. (2007). These, active, materials are driven out of equilibrium by continuous energy input at the level of single particles. Despite the differences in the origin of their activity, several features such as long-ranged ordering patterns, hydrodynamic instabilities and turbulent-like flow dynamics are prevalent among different active systemsMarchetti et al. (2013); Koch and Subramanian (2011); Buttinoni et al. (2013); Doostmohammadi et al. (2018); Ramaswamy (2010); Opathalage et al. (2018). The extent and origins of this universality is not yet fully understood, and controlling the pattern formation and flow dynamics is a way to gain insight into the underlying physics that governs the behaviour of active matter.

Though active systems have been a subject of increasingly intense study in the past decade, limited progress has been made towards utilising them for real-life applications. One of the major hurdles is in controlling their flow behaviour, but manipulation of the geometry confining the active fluid may be a way forward. For example, a net mass flux through rectangular channels containing microtubule bundles powered by kinesin motors has been achieved by tuning the aspect ratio of the channels Wu et al. (2017), with potentials in microfluidics applications.

In this article we use continuum simulations to investigate the behaviour of active fluids confined in a two-dimensional channel. We consider active nematicsDoostmohammadi et al. (2018), a class of active fluids such as microtubule bundlesSanchez et al. (2012), elongated bacteria Volfson et al. (2008); Doostmohammadi et al. (2016); Shi et al. (2014) and fibroblastKemkemer et al. (2000); Duclos et al. (2014), epithelialSaw et al. (2017), and neural progenitor stem cellsKawaguchi et al. (2017), that manifest nematic liquid crystal features including nematic order and topological defects. These have been shown to undergo a spontaneous transition to a flowing state at a finite activityVoituriez et al. (2005); Duclos et al. (2018), and then to a ‘dancing’ state at higher activities Shendruk et al. (2017), where velocity vortices drive topological defects along oscillatory paths. Further increase in activity drives the system to an active turbulent state where a chaotic flow regime is accompanied by the continuous creation, movement and annihilation of $+1/2$ and $-1/2$ topological defects. The rate of defect creation and annihilation balance such as to maintain a statistically steady-state defect density Sanchez et al. (2012); Thampi et al. (2013). The defects are not simply advected, they have their own dynamics: $+1/2$ defects continuously move due to their innate self-propulsion Giomi et al. (2013), and both $+1/2$ and $-1/2$ defects are strong sources of active flows Thampi et al. (2014).

However, results so far have been for a system where a thermodynamic free energy imposes strong nematic ordering. For example, spontaneous flow transitions in confined active nematics Marchetti et al. (2013); Ramaswamy (2010) and emergence of various flow states such as dancing state Shendruk et al. (2017) and its transition to turbulence Doostmohammadi et al. (2017) all assume the presence of nematic order in the system in the absence of activity. Many experimental systems have no nematic ordering in the absence of activity. We therefore consider temperatures above the nematic transition. We argue that local nematic order can be generated by an active flow field, even in the absence of thermodynamic ordering in such systems, and describe the flow configurations of such a flow-stabilised system in a channel. By comparing results for thermodynamically-stabilised and flow-stabilised active nematics we show that the key parameters controlling the flow states are the activity, the channel width and whether the active system is flow aligning or flow tumbling.

We describe our governing equations and simulation approach in Sec. II. We then summarise how active nematics can become ordered both below and above the nematic transition temperature. The possible flow configurations are characterised in Sec. III. In Sec.  IV we present phase plots showing the different flow states that are stabilised as activity, the flow alignment parameter and the temperature are varied.

II Equations of Motion

We first summarise the continuum description of an active nematic, and the simulation approach we use to solve the equations of motion. In particular we discuss how extensile flow fields can lead to local order in an active nematic even in the absence of thermodynamic ordering.

II.1 Governing equations

We consider the dynamics of an incompressible, active fluid of density $\rho$ and velocity $\mathbf{u}$ within the continuum approximation. The fluid is driven by active constituents which are assumed to be apolar in nature Marchetti et al. (2013). We define a tensor order parameter

[TABLE]

where $\mathbf{n}$ is the director field (of unit length) and $q$ gives the strength of nematic ordering of the active constituents de Gennes and Prost (1995). The free energy of the system is defined using a Landau-de Gennes approach in the same way as for a passive nematic liquid crystal. It consists of two contributions, a bulk term

[TABLE]

where the coefficients $A,B,C$ determine the order in the system, and a gradient term modelled using the Frank-Oseen free energy

[TABLE]

which describes the free energy cost of any deviations from uniform nematic ordering. In eqn (3) we have assumed a single elastic constant, $K$. The total free energy of the active nematic is $F=F_{\textnormal{LDG}}+F_{\textnormal{FO}}$.

The mass and momentum conservation equations for the fluid areYeomans (2016); Beris and Edwards (1994)

[TABLE]

[TABLE]

where the stress tensor

[TABLE]

The viscous stress is modelled using a Newtonian constitutive relation,Batchelor (2000)

[TABLE]

where $\mu$ is the shear viscosity of the fluid and $\mathbf{E}$ is the symmetric part of the velocity gradient tensor.

The active stress is taken to be proportional to the orientational order Simha and Ramaswamy (2002)

[TABLE]

where $\zeta$ describes the strength of the activity. When local order is developed in the fluid, either due to activity itself or due to lowering the temperature, this term can drive active flows. The sign of the activity coefficient $\zeta$ determines the symmetry of the stresslet flow field generated by the active constituents Simha and Ramaswamy (2002). $\zeta>0$ represents extensile systems which draw fluids in from their sides and push it out along the director field, e.g. bundled microtubules driven by motor proteins Edwards and Yeomans (2009); Ramaswamy (2010). Reversing the direction of the active stress gives rise to contractile systems, $\zeta<0$.

The elasticity due to orientational order in the fluid generates stresses, that also arise in a passive liquid crystal. These are given by de Gennes and Prost (1995); Beris and Edwards (1994)

[TABLE]

where $\lambda$ is the aligning parameter and $H_{ij}$ is the molecular potential field calculated from the free energy, $F$,

[TABLE]

Both the active and the passive stresses are dependent on the order parameter $\mathbf{Q}$. Any gradient in $\mathbf{Q}$, either in the orientation ($\mathbf{n}$) or strength ($q$), can result in a fluid flow. Moreover the fluid flow can, in turn, change the orientational order in the system. The dynamics of $\mathbf{Q}$ follows the evolution equation

[TABLE]

where $\Gamma$ is the orientational diffusivity and $S_{ij}$ is the generalized advection term

[TABLE]

with $\bm{\Omega}$ the antisymmetric part of the velocity gradient tensor. $S_{ij}$ describes the response of the fluid to spatial gradients in the velocity field. Equation 12 indicates that the aligning parameter $\lambda$ determines the relative importance of the response of $\mathbf{Q}$ to the extensional and rotational parts of the flow field.

II.2 Simulation details

The equations of motion summarised in Sec. II.1 are solved using a hybrid lattice Boltzmann algorithm. In this method, the fluid dynamics is obtained by solving a discretised Boltzmann equation which yields the Navier-Stokes equations in the continuum limit Denniston et al. (2002). The evolution of the $\mathbf{Q}$ tensor is followed using the method of lines wherein the spatial derivatives are discretised using second order central differences and the temporal integration of the resulting algebraic equations is performed using a fourth order Runge-Kutta method Doostmohammadi et al. (2016).

A schematic of the two-dimensional domain used for the simulations is shown in Fig. 1. The active fluid is sandwiched between two parallel plates of length $L$ separated by a distance $W$. Typically $L>>W$ so that the results are independent of the channel length. To facilitate the discussions, we adopt Cartesian coordinates, $x$ along the channel length, and $y$ in the perpendicular direction, across the channel. Periodic boundary conditions for the flow and orientation tensor are imposed along the $x$ direction. No-slip boundary conditions are used for the velocity field at the walls. Homogeneous boundary conditions are used to align the director field parallel to the walls with a magnitude equal to that of the passive liquid crystal at the chosen simulation temperature. Hence $\mathbf{Q}=0$ on the walls for temperatures above the passive nematic ordering temperature. To mimic the experimental conditions, all the simulations are started with a quiescent fluid and with a random initial configuration for the orientation tensor. The range of parameters used in the simulations are listed in Table 1.

II.3 Mechanisms for nematic ordering

In a passive, equilibrium liquid crystal, the degree of orientational order is directly related to the inter-particle interactions and the temperature of the system. As the temperature is lowered the liquid crystal becomes more orientationally ordered to increase its positional entropyFrenkel (2014). This dependence of the order on temperature can be modelled phenomenologically through the material coefficients, $A$, $B$ and $C$ in the Landau-de Gennes bulk free energy de Gennes and Prost (1995) (eqn 2). In the absence of non-equilibrium effects, i.e. zero activity, the parameter $\overline{T}=27AC/B^{2}$ controls the state of the system in the isotropic-nematic phase-space as shown schematically in Fig. 2(a).

In an active nematic, the activity itself can provide a second route to nematic ordering, even in the absence of any thermodynamic tendency towards order i.e. at infinite temperature Thampi et al. (2015)(see Fig. 2(b)). This occurs because small fluctuations in the director field can lead to local nematic alignment which generates a net stresslet flow. For an extensile system, and $\lambda>0$, the flows can rotate and orient the nematogens, further strengthening the order in the system.

This mechanism is apparent from the continuum model of nematic liquid crystals when the momentum balance is primarily between active and viscous stresses, $2\mu E_{ij}=-\zeta Q_{ij}$. Then eqn. 11 may be rewritten as

[TABLE]

where

[TABLE]

for $A=B=C=0$. Terms in $\mathbf{H}^{\textnormal{active}}$ on the right hand side of eqn. (13), which arise from the generalized advection term $\mathbf{S}$, resemble those that appear in the molecular field $\mathbf{H}$ due to the Landau-de Gennes free energy $F_{\textnormal{LDG}}$ (see eqns. (2) and (10)). Thus, activity acts as an effective free energy with the parameter $\frac{\lambda\zeta}{3\Gamma\mu}>0$ determining the strength of the ordering.

In general, both thermodynamic and flow ordering will contribute to determining the strength of the nematic order in the active system. In the following when we need to distinguish between the two cases we will use the terminology thermodynamically-stabilised active nematic to refer to $\overline{T}<\frac{27}{24}$ where the nematic state is (meta)stable in the passive system, and flow-stabilised active nematic to refer to $\overline{T}>\frac{27}{24}$ where activity-induced flow is needed to stabilise the nematic order.

III Active flows in a channel

III.1 Flow states

We next describe the different flow configurations that can be generated when an active nematic is confined within a long channel. For small activity there is no motion in a channel of finite width; the fluid undergoes a spontaneous flow transition at a threshold activity $\zeta_{c}$Voituriez et al. (2005). Depending upon parameters, four distinct types of flow fields can be obtained as activity is increased beyond this threshold. These states, (i) unidirectional (ii) oscillatory (iii) ‘dancing’ flow and (iv) turbulent flows, are shown in Fig. 3l. Change in activity is reported in terms of the activity number, $W\sqrt{{\zeta}/{K}}$ which is the ratio of the two characteristic length scales, the channel width and the length scale associated with active vortices, $\sqrt{{K}/{\zeta}}$.

Unidirectional flow, illustrated in Fig. 3l(a) for an activity slightly larger than $\zeta_{c}$, is characterised by streamlines aligned parallel to the channel walls. The velocity field is translationally invariant along $x$ and ${u}_{x}={|\mathbf{u}|}$ varies from zero at the channel walls to a maximum in the centre of the channel. Topological defects are not generated in unidirectional flows.

Oscillatory flow: As the activity is increased the component of velocity normal to the channel walls becomes non-zero. However, $u_{y}$ is typically smaller than $u_{x}$. The corresponding streamlines, vorticity field and director field are shown in Fig. 3l(b). An alternating pattern of vorticity patches with opposite sign starts to develop, and the corresponding director field shows locally nematic regions. Topological defects are occasionally present.

Dancing flow:Shendruk et al. (2017) A one-dimensional flow vortex lattice, with topological defects following regular oscillatory paths along the channel, is observed at higher activities when the length scale of active vortices, which decreases with increasing activity, becomes of order the channel width (Fig. 3l(c)). The channel is filled by a row of velocity vortices which alternate in direction. In contrast to the unidirectional and oscillatory states, here the director field develops a regular pattern of motile topological defects. $-\frac{1}{2}$ defects are confined close to the walls due to elastic interactions with the boundaries. However, because of their innate self-propulsion, $+\frac{1}{2}$ defects continuously move along the channel, in either direction, along well-defined, oscillatory trajectories. Fig. 4(a) compares the distribution of $+\frac{1}{2}$ and $-\frac{1}{2}$ defects across the channel. Once the oscillation is established there is no nucleation/annihilation of topological defects. As the defects move there is an accompanying oscillation in the relative strengths of the positive and negative vortices. Hence ${u}_{x}$, and to a lesser extent ${|\mathbf{u}|}$ vary with time.

Turbulent flows, with properties summarised in Fig. 3l(d), are seen when the typical size of the active vortices is smaller than the channel width. Flow vortices, which are not uniform in size Guillamat et al. (2017); Giomi (2015), are short-lived and move around in the channel. The resulting streamlines appear very similar to those of high Reynolds number turbulent flows, hence the name active or mesoscale turbulence Wensink et al. (2012). $\pm\frac{1}{2}$ defects are continuously created in the bulk of the fluid, move around and then undergo annihilation. The defect distribution, illustrated in Fig. 4(b), shows less asymmetry in the locations of $\pm\frac{1}{2}$ defects, although the -1/2 defects still tend to reside near the walls.

III.2 Flow order parameters

In the previous section we have seen that field variables are different for the different flow states. Therefore, to distinguish between the different flow configurations in the simulation, it is helpful to define order parameters that estimate the relative net flow along $x$ and along $y$

[TABLE]

[TABLE]

where $\langle\cdots\rangle_{x}$, $\langle\cdots\rangle_{y}$ denote averages along the length and width of the channel. For unidirectional flows, $u_{y}=0$ everywhere and thus $\phi_{1}=1,\phi_{2}=0$. In case of oscillatory flows, $u_{x}$ varies throughout the channel, but it dominates over $u_{y}$. We identify $0.7<\phi_{1}<1$ and $\phi_{2}<0.2$ as oscillatory flows in the simulations. For the dancing state of alternating vortices neither $u_{x}$ nor $u_{y}$ is negligible. However, $u_{y}$ is persistent across the channel and $u_{x}$ changes direction along the channel. Visual inspection shows that small $\phi_{1}$ and $0.45<\phi_{2}<0.65$ identifies the dancing states. For turbulent flow, both $u_{x}$ and $u_{y}$ fluctuate throughout the channel yielding $\phi_{1}$ small and $\phi_{2}<0.45$. These identifications are summarised in Fig. 5.

IV Results

We now present results showing how the flow configurations change as the activity number and the aligning parameter are varied. The aligning parameter $\lambda\propto(r^{2}-1)/(r^{2}+1)$ is related to the aspect ratio of the active particles, $r$. We consider $\lambda>0$, corresponding to rods. Results for thermodynamically-stabilised active nematics will be presented first, followed by a comparison to those for flow-stabilised fluids. We will then discuss the relationship between the two cases.

The flow configurations identified at each point in the phase plot are obtained from $10$ different realisations of each simulation. Different realisations were achieved by changing the seed that prescribes the random distribution of the director field in the initial conditions. We found that the final flow state is often strongly dependent on the initial conditions, and we indicate this by using two or more superimposed symbols. (This metastability was less pronounced in reference Shendruk et al. (2017) because of the nematic initial conditions of the director field.) Moreover unidirectional flow states only exist over a very narrow range of parameters. Therefore, for clarity, they are not distinguished from oscillatory flow states in the phase plots.

IV.1 Flow states of a thermodynamically-stabilised active nematic

Fig. 6 shows the different flow states of a thermodynamically-stabilised active nematic fluid as a function of activity number and aligning parameter. For sufficiently high aligning parameter ($\lambda>\sim 0.7$), the phase sequence with increasing activity number is no flow $\rightarrow$ oscillatory flow $\rightarrow$ dancing flow $\rightarrow$ active turbulence. However, for smaller $\lambda<\sim 0.6$, the quiescent state transitions directly to the dancing state and then, at a higher activity number, to turbulence.

At many points in the phase diagram, and particularly close to the borders between different flow states, the final steady state depends sensitively on the initial conditions. Interestingly, at the borders, the flows obtained are not always the ones pertaining to the neighbouring regions: thus the smearing of the boundaries is not just a consequence of hysteresis associated with the transitions. For example, turbulent flow occurs between the oscillatory (region I) and dancing (region IIB) flow states; and the oscillatory state, which corresponds to a net flow, occurs between the dancing (region IIB) and turbulent (region III) regimes, which do not generate any net flow.

IV.2 States of a flow-stabilised active nematic

We now compare results for an active nematic where the ordering is predominantly from the extensile flow field ($\overline{T}=\frac{27}{20}$). A phase plot showing the stable flow configurations as a function of activity number and aligning parameter is presented in Fig. 7. For $\lambda>\sim 0.7$ the phase sequence with increasing activity number is no flow $\rightarrow$ dancing flow $\rightarrow$ turbulence. For $\lambda<\sim 0.6$ the quiescent state transitions directly to turbulence.

There are several differences between the behaviour in the thermodynamically-aligned and flow-aligned nematics. In the latter case the stationary state (no flow state) becomes unstable at a considerably larger value of the activity number, as the activity has to drive first ordering, and then the instability to flow. Once flow is possible the dancing or turbulent states are more likely to be stabilised directly and no unidirectional or oscillatory flow states were observed within the range of parameters studied. Another striking difference is that, unlike the thermodynamically-stabilised regime, the final flow states in the flow-stabilised case were independent of initial conditions for all the parameters we tested.

By varying the effective temperature, $\overline{T}$, it is possible to crossover from a system where the nematic ordering is predominantly thermodynamically driven to one where it is predominantly due to the flow. Fig. 8 shows how the flow states vary with temperature for different activity numbers for alignment parameter $\lambda=0.7$. The threshold activity number for flow generation is $\approx 10$ for $\overline{T}<\sim 1$. The threshold sharply increases around $\overline{T}=1$ and then continues to increase mildly with further increase in temperature. The oscillatory states become unstable for higher temperature, and metastability decreases with increasing temperature. However, the pattern of states is complex. Therefore, in the next subsection, we show that a different choice of variables allows us to present a more unified picture of the results for varying temperatures.

IV.3 Unifying the observations at all temperatures

Motivated by the dependence of the flow configurations on both the $\lambda$ and the degree of nematic ordering that is induced by thermodynamics and flow, we recall that the response of a nematic rod to flow gradients depends on both $\lambda$ and on the magnitude of the order parameter $q$. If the order parameter is assumed to be of a constant magnitude then we may write Marenduzzo et al. (2007)

[TABLE]

The director field aligns at a constant angle to a simple shear flow when $\lambda_{1}>1$, otherwise it tumblesEdwards and Yeomans (2009). We will now investigate whether $\lambda_{1}$ provides a more suitable control parameter for the flow behaviour.

To obtain a representative value of the magnitude of order parameter for each simulation, we define $q_{avg}=\langle\langle q\rangle_{x}\rangle_{y}$, averaging over the simulation domain and time. No-flow configurations are excluded from this calculation as $q_{avg}=0$ for these states for $\overline{T}>1$.

A phase plot showing the stable flow configurations as a function of activity number and $\lambda_{1}$ is presented in Fig. 9. This figure includes data from Figs. 6 - 8. Note that the results from higher temperatures lie on the right-hand side of the plot because spontaneous flow only occurs for larger activity numbers.

To a good approximation the flow states divide into four regions. For flow aligning systems, $\lambda_{1}>\sim 1$, there is a transition from oscillatory to dancing flow with increasing activity number, whereas for flow tumbling nematics, $\lambda_{1}<\sim 1$, the transition is from dancing to turbulent flow. This is reasonable as a tumbling configuration will more easily form flow vortices.

Thus we are able to present the results for thermodynamically-stabilised and flow-stabilised active nematics within the same framework. Our results show that the key parameters controlling the transitions are the activity number, and whether the nematogens are flow-aligning or flow-tumbling. The origin of the nematic ordering (thermodynamic- or flow-stabilisation) is not important.

V Summary and discussion

Typically flow instabilities in active nematics lead to a chaotic flow field characterised by high vorticity and motile topological defects. In the previous work a flow lattice was observed either in a flow stabilized active nematic laid on a substrate Doostmohammadi et al. (2016) or in a thermodynamically stabilized, tumbling active nematic confined between two plates Shendruk et al. (2017). Here we show that for both thermodynamically and flow stabilised active nematics regular flows can be obtained by confinement. In two dimensions these include laminar flow, oscillatory flow and a ‘dancing’ flow state where +1/2 topological defects move along the channel. In this paper we have presented results showing the region of stability of the different flow configurations, concentrating in particular on comparing systems below and above the isotropic-nematic transition, $\overline{T}=1$, where the nematic ordering is predominantly due to thermodynamic forces or active flow stabilisation, respectively. This was motivated by noting that many experimental systems are not nematic until active forcing is introduced. Similar flow states are seen at all temperatures, but there is a much lower probability of obtaining net flow for $\overline{T}>1$. We found that the final flow configuration is strongly dependent on the choice of initial conditions for $\overline{T}<1$. However, this metastability is absent at higher temperatures. Moreover the critical activity required to trigger hydrodynamic instabilities is much higher above $\overline{T}=1$. This is because the activity has to first stabilise local nematic order before a spontaneous flow transition can occur.

By varying the temperature, activity and flow alignment we found that the key parameters controlling the transition are the activity number and a scaled alignment parameter. The activity number, $W\sqrt{{\zeta}/{K}}$, is a dimensionless variable characterising the ratio of the channel width to the vortex length scale associated with the active fluid. Larger activity numbers correspond to smaller vortices and hence increases in the activity number favour first the dancing state and then active turbulence. The scaled alignment parameter, $\lambda_{1}$, takes into account both the shape of the nematogens and the magnitude of the nematic order and controls whether the system is flow aligning or flow tumbling. The dancing and turbulent states appear at lower activity numbers in flow tumbling systems because it is easier to form vortices.

The transition to laminar flow has been observed in confined confluent cell layers Duclos et al. (2018), and there have been indications of vortex structures in experiments studying both cellsNier et al. (2016); Duclos et al. (2018)and bacteriaWioland et al. (2016). It would be of interest to look at other active systems in channels, although experiments may be difficult due to metastability effects. Recent experiments on microtubule gels in three-dimensional channels with rectangular cross sections have shown that the existence of net flow depends on the aspect ratio, but not the dimensions, of the channel cross sectionWu et al. (2017). This system is in the flow-aligning regime, above the isotropic-nematic transition temperature, and extending our simulations to three dimensions might help to explain this surprising result.