Pressure-driven changes to spontaneous ﬂow in active nematic liquid crystals

. We consider the eﬀects of a pressure gradient on the spontaneous ﬂow of an active nematic liquid crystal in a channel, subject to planar anchoring and no-slip conditions on the boundaries of the channel. We employ a model based on the Ericksen-Leslie theory of nematics, with an additional active stress accounting for the activity of the ﬂuid. By directly solving the ﬂow equation, we consider an asymptotic solution for the director angle equation for large activity parameter values and predict the possible values of the director angle in the bulk of the channel. Through a numerical solution of the full nonlinear equations, we examine the eﬀects of pressure on the branches of stable and unstable equilibria, some of which are disconnected from the no-ﬂow state. In the absence of a pressure gradient, solutions are either symmetric or antisymmetric about the channel midpoint; these symmetries are changed by the pressure gradient. Considering the activity-pressure state space allows us to predict qualitatively the extent of each solution type and to show, for large enough pressure gradients, that a branch of non-trivial director angle solutions exists for all activity values.


Introduction
Active fluids exhibit a continuous generation of internal energy, as in swimming bacteria and microtubuleforming suspensions [1], allowing for spontaneous flow generation [2][3][4][5]. This continuous energy production leads to the active agents (the bacteria or microtuble-forming motors) exerting a stress on the background fluid and creates a system that is always away from thermodynamic equilibrium. In active nematic liquid crystal fluids, the flow-generating agent is anisotropic (defined by, for instance, the long axis of a bacterium or microtubule), with a macroscopic symmetry of a liquid-crystalline phase. The internal active stress can then lead to a vast array of interesting effects such as hydrodynamic and distortional instabilities and non-equilibrium defect formation [6][7][8][9]. Recent experiments and theoretical investigations have also been undertaken into the control and designing of active nematics through a magnetic field [10][11][12], a technique also used in the manufacturing and design of liquid crystal displays. Reviews of the literature in this area can be found in [12][13][14][15]. In recent years, researchers have explored the possibility of activity-driven microfluidic systems [16][17][18][19][20][21], including the use of swimming bacteria to assist in targeted drug delivery [22], bacterial ratchet motors [23], active turbulence-powered rotators [20] and even active a e-mail: g.mckay@strath.ac.uk fluid-based logic circuits [16]. The effects of confinement in such systems is of crucial importance, with the frustration due to the fluid-boundary interaction leading to organisation [17][18][19]21,24]. In classical microfluidics or larger scale fluid processing applications, rather than internal activity driving flow, the fluid is often forced by the application of an external pressure gradient. While there are many studies of how low numbers of swimmers interact within a pressure gradient, there are relatively few investigations of how the activity-driven and pressure-driven flows will interact in a dense continuum of swimmers. A number of authors have considered models of dilute suspension of active swimmers in a pressure gradient [25][26][27] and dense continuum models of channel flow, albeit in a tumbling nematic regime [28]. Others have reported experimental studies of similar systems [29], while there are also general reviews of the rheology of active fluids in [13,30,31]. In this paper we add to these previous investigations, considering a dense continuum model of an active fluid subjected to an external pressure gradient along a channel.
The similarities between the macroscopic symmetries, flow and defect patterns generated in active fluids and those of elongated rod-like molecules in liquid crystals mean that continuum hydrodynamic models of nematic liquid crystals have commonly been adopted in the theoretical modelling of active fluids [2,3]. Many continuum models of active nematics are, therefore, based on modifying the Ericksen-Leslie equations for liquid crystal flow [32][33][34], in which orientational ordering and fluid flow are coupled and governed through a set of partial differential equations derived from balances of mass and momentum, both linear and angular. The effects of activity are built into the system by modifying the Ericksen-Leslie equations via additional stress terms which can drive the system out of thermodynamic equilibrium [35]. In both active and inactive nematic liquid crystals, the average direction of the constituent parts (i.e., molecules in a classic inactive liquid crystal, or, for instance, bacteria in an active fluid) is represented by a unit vector known as the director, n(x, t). Other authors have also considered an alternative mathematical description of active nematics which uses both average orientational ordering and a measure of the order about that average direction, the general form of which is a second-order tensor, commonly known as the Q-tensor [6,7,36,37]. This approach is particularly useful when investigating active systems with defects, which cannot be treated with director-based theories of liquid crystals.
In the present work we will consider non-polar active agents and so the active stress term we employ will also obey the symmetry n → −n. Other work has considered polar constituent parts, where it is possible for the macroscopic system to be polar. This reduction of symmetry allows additional activity terms to be present, such as selfpropelling velocity and active viscosity terms [4,5,38,39]. Various different sources of activity have previously been considered in the mathematical modelling of active fluids, with active terms added to the governing equations for non-active liquid crystals. For example, Yang and Wang [39] consider the influence of active viscous stress tensor and self-propelling speed terms for channel flows of active polar liquid crystals. The former is added to the stress tensor, whereas the latter enters the governing equation for the polar director field. More recently, Duclos et al. [40] and Hoffmann et al. [41] consider spontaneously induced flows due to an active stress term unique to active cholesteric liquid crystals.
The specific active stress contribution we consider, denoted by σ a , is the most common form of activity term. First introduced by Simha and Ramaswamy [42], it has also been considered in many subsequent papers, for instance [39][40][41]. It can be obtained by considering each individual active agent as a permanent force dipole, leading to local active stresses proportional to the second moment of the molecular orientation. (A derivation can be found in Thampi and Yeomans [1].) In terms of the director, the active stress can be written as where the outer product is defined as n ⊗ n = n i n j and n i (i = 1, 2, 3) is the i-th component of the director. The coefficient ζ is termed the activity parameter, which can be positive or negative, and has the dimensions of pressure. The magnitude of ζ quantifies the degree of activity or, equivalently and more specifically, the stress the active agents exert on the background fluid. The sign of ζ distinguishes how the active agents behave relative to the surrounding fluid, with the agents either pushing the fluid out or pulling the fluid in along their long axis. This simple description of "pushers" and "pullers" to describe active agents is commonly replaced by the terms "extensile" and "contractile", respectively. A schematic illustration of these two contrasting behaviours for active agents is shown in fig. 1, where for extensile agents ζ < 0 and for contractile agents ζ > 0 [4,38,39,43]. Note that the opposite sign of the activity parameter has been assumed in a number of publications, with a corresponding negative active stress in eq. (1) [1][2][3]36].
In this paper we explain the large-magnitude activity asymptotics of states discovered previously, for example in [2,3,37]. We also investigate the pressure-driven states in active nematics where we show how the stability of equilibria, in particular those of a certain symmetry, can be promoted through an applied pressure gradient, which, given the one-dimensional nature of the system, is constant. This work will aid future experiments, and theoretical investigations, in controlling and designing active nematic systems through an external pressure gradient.

Mathematical model
We consider an active nematic liquid crystal, confined between two parallel plates at z = 0 and z = d, and subject to a pressure gradient parallel to the x-direction (see fig. 2). The nematic director field n is constrained to lie in the (x, z)-plane and is described by n = (cos θ, 0, sin θ), where θ(z, t) is the director angle measured with respect to the x-direction. We assume that surface treatment of the bounding plates anchors the director such that it is forced to lie in the x-direction at z = 0 and z = d, a situation usually termed infinite anchoring. The velocity of the fluid, v(z, t), is assumed to be in the x-direction and satisfies the no-slip condition at the plates. The assumptions of a director confined to the (x, z)-plane and rectilinear flow are related to the nematic being "flowaligning" and strongly anchored in the x-direction at the boundaries. For such liquid crystals, rectilinear flow promotes alignment within the shear plane at the Leslie angle to the flow direction [44], while rotation within the plane, together with the presence of the solid boundaries, restrict the flow to a single direction. In such a system, director instabilities out of the shear plane could only occur if the nematic is non-flow-aligning [45], a situation that we do not directly consider in the present work. We will also assume that fluid inertia is negligible, an approximation that is valid when the Reynolds number is small or, equivalently in this situation, when the timescale of changes in the velocity are much smaller than the timescale of director rotation [46]. The dynamics of the fluid velocity and director angle can then be modelled by the balance of linear and angular momentum through the Ericksen-Leslie equations for nematic liquid crystals [33,34,44], with the inclusion of the active stress term introduced in eq. (1): where g(θ) = η 1 cos 2 θ + η 2 sin 2 θ + η 12 sin 2 θ cos 2 θ, In eqs. (2)-(5), ρ is the fluid density, K 1 and K 3 are the Frank elastic constants associated with, respectively, splay and bend of the director, p x represents the constant applied pressure gradient, γ 1 is the rotational viscosity, γ 2 is the torsional viscosity, while η 1 , η 2 , η 12 are Miesowicz viscosities [47]. The subscripts t, x or z denote the partial derivative with respect to that variable. Finally, as mentioned earlier, ζ is a measure of the activity exhibited by the active nematic liquid crystal. The infinite anchoring and no-slip assumptions on the boundary plates mean that eqs. (2) and (3) will be solved subject to the boundary conditions At this stage it is worth considering the possible symmetries of solutions of eqs. (2) and (3), with boundary conditions (6). (2) and (3) unchanged for any value of the pressure gradient p x , so we expect to obtain solutions for the director angle which are antisymmetric about the channel midpoint z = d/2, together with velocity solutions that are symmetric. For the opposite symmetries, (2) and (3) are unchanged only if p x = 0. We therefore expect to find symmetric director and antisymmetric velocity solutions only for zero pressure gradient. In addition to symmetry/antisymmetry about the channel midpoint, we notice that the governing equations are also unchanged under the transforma- Therefore, changing the sign of the pressure gradient will simply result in a change of sign of both the director angle and velocity. It is possible to decouple eqs. (2) and (3) using the same approach considered in Mottram et al. [46], leading to a single, non-local, dynamic equation for the director angle, namely where In this approach, the velocity is a function of the director orientation, namely The activity parameter ζ enters eq. (7) only through the fourth term whose form, sin θ cos θ, is similar to a magnetic or electric field term in the classic problem of director reorientation during a Freedericksz transition [44], albeit rescaled by the director dependent factor m(θ)/g(θ) and "shifted" by the non-local terms C+D/B. The fifth term in eq. (7) is the contribution from the pressure-driven flowdirector coupling. Non-trivial analytic solutions to the nonlinear, nonlocal partial differential equation in eq. (7) are not possible. However, we can still establish important information about the system, and further progress can be made using certain simplifying assumptions. We first consider two important asymptotic limits for the behaviour of the director in the bulk of the channel: the situation when the magnitude of the applied pressure gradient is large; and when the magnitude of the activity parameter is large.

Analysis: asymptotic solutions for large pressure gradient or large activity
In the case when |p x | ζ/d and , the dynamics of the director in the centre of the channel, i.e. away from any boundary or internal reorientation layers, can be found by neglecting elastic effects, so that The only equilibrium solutions of eq. (15), so that θ t = 0, will be those θ values that satisfy m(θ) = 0. Such solutions only exist in flow-aligning nematics and are θ = nπ ± θ L , where n ∈ Z and θ L = tan −1 ( (γ 2 + γ 1 )/(γ 2 − γ 1 )) is the Leslie angle [44]. For these constant director angle solutions the velocity is then, as is to be expected in this limit, the classic Poiseuille parabolic profile, v(z) The situation for highly active systems, where we take |ζ| p x d and |ζ| K i /d 2 (i = 1, 3), is more complicated and has not previously been considered analytically. In this case we might expect the dynamics in the bulk of the channel, where spatial gradients are negligible and θ = θ(t), to be governed by the equation We can now consider two possible symmetries of the director profile within the channel. If θ is antisymmetric with respect to the centre of the channel, then the integrals C and D in eqs. (10) and (11) are both zero. Therefore, in the bulk of the channel, away from boundary and internal reorientation layers where director distortions are relatively large, antisymmentric equilibria solutions for the director angle must satisfy m(θ) sin θ cos θ = 0, leading to the possibilities θ = ±θ L , 0 or π/2 rad (plus all π rotations of these angles).
For symmetric director profiles, C = 0 and the situation is more interesting as we find a novel solution in which the equilibrium orientation is determined through non-local effects. From eqs. (10) and (11), a symmetric director profile means that C = sin θ cos θ and D = 0, leading to zero on the right-hand side of eq. (16) for all values of θ. Rather than all director angles in the bulk of the channel being equilibria, as this result would suggest, it is necessary to consider the non-local behaviour in the system. To make analytic progress in this case, we use the relatively standard one-constant approximation, K 1 = K 3 = K, introduce a non-dimensionalised coordinate ξ via z = dξ and consider the static situation. The governing equation for the director angle is then where β = ζd 2 /K is the non-dimensionalised activity parameter, and Through an integration of eq. (17) from ξ = 0 to ξ = 1, it is straightforward to see that A + βD = 0, so we now consider solutions of which will be solved subject to the boundary conditions θ(0) = 0 = θ(1). We seek the solution to eq. (22) in the limit of large activity, i.e. as β → ∞. The integrodifferential equation (22) is, however, non-local through the integral (10). While we cannot derive a direct analytic solution, an approximate solution is possible. Consider first the Hamiltonian for the problem given by eq. (22) Note that the integral in eq. (23) can be determined analytically, but is too lengthy to include here. It is clear from eq. (23) that we would expect θ ξ ∼ √ β as β → ∞ and will therefore set θ ξ (0) = √ β p, where p is to be determined. We are considering symmetric director angle profiles, so that θ ξ (1/2) = 0, and we have θ(0) = 0 from the boundary condition. Therefore, setting θ(1/2) = θ m , we consider trajectories in the phase plane that connect the points (θ, θ ξ ) = (0, √ β p) to (θ m , 0). Since the Hamiltonian is constant along phase plane trajectories we have , which can be rearranged to calculate p in terms of θ m : Returning to eq. (22), and integrating between ξ = 0 and ξ = 1/2, we also obtain Guided by the form of the derivative of the director angle at ξ = 0, namely θ ξ = √ β p, and the expectation that the director angle will be constant in the bulk of the region, we now make an approximation for the form of θ(ξ) that allows for a linearly increasing director angle close to the boundary The integrations in eq. (25), including those within C, may be performed using the change of variables from ξ to .
(31) Equating eqs. (30) and (31) therefore provides an implicit equation for θ * , which is dependent on the material viscosities through m(θ) and g(θ). This is an analytical result, and valid for any nematic-forming liquid. The required viscosity measurements have not been undertaken for an active nematic, so, in order to estimate the value of θ * , we will use viscosity parameter values for the liquid crystal 5CB [44]. For these parameters we obtain the value θ * ≈ 1.212 rad. This approximate analytic value compares well to the numerically computed value of θ * ≈ 1.199 rad in sect. 4.2, and even better to the value θ * ≈ 1.207 rad which we find if we assume a one-constant approximation of K = (K 1 +K 3 )/2 while retaining the other 5CB parameters. Although in this paper we consider flow-aligning active nematics, we note that for non-flow-aligning nematics (as studied by, for instance, [28]) it is possible to make analytical progress from eq. (22). This situation is included in the appendix for reference.

Numerical results
After investigating the limiting solutions for large pressure gradients and activity in the previous section, we now consider the behaviour of director and velocity profiles for finite parameter values. We will briefly recap the situation for zero pressure gradient, before examining the effect of applying a pressure-driven flow. We compute numerical steady state solutions (θ(z), v(z)) of the full equations (2) and (3). In order to obtain numerical solutions and continue along solution branches, we have employed both the finite-element package COMSOL [48] and the MATLABbased bifurcation analysis package MATCONT [49].
The material parameters, i.e. elastic constants and viscosities, of active nematics have not yet been fully characterised. Therefore, as mentioned previously, we use the material parameters measured for the liquid crystal 5CB [44] and a channel width of d = 10 μm. There is a significant difference in length scales between the constituent agents for 5CB (nanometre size) and an active nematic (micron size), and the actual values of the material parameters may be different. However, typical experimental measurements for much larger molecules in a nematic phase, i.e. polymer nematics, show that their material parameter values can be similar to those of low molecular weight nematics [50], at least in the ratios we consider in the non-dimensional quantities. We therefore use the well characterised 5CB values and suggest that the behaviour we find below will be qualitatively the same for an active nematic.

Extensile active nematic
As is well known [2,3,40], for zero pressure gradient the trivial state (θ, v) ≡ (0, 0) is a solution for all values of the activity parameter, but is unstable for activities ζ < ζ c , where ζ c = (8π 2 K 1 η 1 )/((γ 1 +γ 2 )d 2 ). The parameters used here correspond to ζ c = −13.87 Pa, while the Leslie angle is θ L ≈ 0.2 rad. Close to this value of activity, there are two modes of instability where the director profile is either antisymmetric or symmetric with respect to the centre of the cell. The respective velocity profiles exhibit the opposite behaviour. However, there are also branches of stable solutions that do not originate from bifucations from the trivial state. As a reference, in fig. 3 we show the three types of solutions found in the absence of a pressure gradient. The director angle profiles in fig. 3(a), (c) are the modes bifurcating from the trivial state as well as the trivial state itself. The equilibria in fig. 3(e) correspond to symmetric director angle profiles which do not bifurcate from the trivial branch and are associated with large elastic energies due to their spatial gradients. Note that, although not shown in fig. 3 but as suggested at the end of sect. 2, negative versions of all the solutions in fig. 3 also exist due to the symmetry θ(z) In fig. 4 we have plotted solution branches for the trivial solution and all three non-trivial modes shown in fig. 3 as the activity varies. In plotting these branches, we use two measures of the director angle solution in order to characterise the symmetry of θ(z), namely  where · represents integration across the channel from z = 0 to d, and (33) Here, the subscripts "o" and "e" refer to, respectively, odd and even (or antisymmetric and symmetric) contributions with respect to the centre of the cell. To help focus on the different types of mode, the equilibria branches in fig. 4 have also been projected onto the shaded horizontal and vertical planes.
Using the continuation package MATCONT for variable activity, we have calculated the eigenvalues of the Jacobian of the discretised numerical system of equations for each equilibrium solution. Solution branches for which all eigenvalues are negative, and thus the system is stable, are indicated by a solid (black) curve. If exactly one eigenvalue is positive, the branch is presented as a dashed (blue) curve, while a dotted (red) curve corresponds to two positive eigenvalues. For any branch with a positive eigenvalue the solution is unstable.
The modes in fig. 3(a), (b) and (c), (d) are all to be expected, with director distortion being linked to activityinduced flow, and have been found previously in, for example, Marenduzzo et al. [37]. In fig. 3(a), (c) we see that for large activities (in magnitude), the director aligns in the bulk of the channel, i.e. away from boundary and in-ternal reorientation regions, at ±θ L ≈ ±0.2 rad, which was predicted as one possibility in sect. 3. The mode in fig. 3(e), (f) is less well studied, perhaps because the solution branch is not connected to the trivial state and contains high gradients in θ (note the different scale of the vertical axis in fig. 3(e) compared to fig. 3(a), (c)). For all solutions in fig. 3, the director angle exhibits flow alignment to θ L in regions of positive shear, v z > 0, and flow alignment to −θ L or (π − θ L ) in regions of negative shear, v z < 0. As predicted in sect. 3, there will also be higher-order mode solutions similar to fig. 3(a) and fig. 3(c) at larger values of the activity parameter magnitude, but where the director angle alternates between θ L and −θ L an increasing number of times. We also predict that there will be equivalent higher-order mode solutions that are not connected to the trivial solution branch, similar to fig. 3(e) but with the director angle alternating between (nπ + θ L ) and (mπ − θ L ) for n, m ∈ Z. However, all these higher-order modes will involve large elastic distortions and may be unstable or metastable.
Having produced the reference bifurcation diagram for the zero pressure gradient case, fig. 4, we now consider how the introduction of a non-zero pressure gradient alters the (θ(z), v(z)) profiles and bifurcation structures. We will concentrate here on the trivial state and its bifurcations, so fig. 5 focuses on how a pressure gradient affects the equilibrium solutions seen previously in fig. 3(a)-(d) for a particular value of the activity parameter, ζ = −20 Pa. In a non-active Newtonian fluid, the addition of a negative pressure gradient, p x < 0, would lead to a parabolic velocity profile with a maximum velocity in the centre of the channel, similar to the flow shown in fig. 3(b). For an active nematic, therefore, such a pressure gradient reinforces the positive flow velocity associated with an antisymmetric director angle solution, thus enhancing the alignment with the Leslie angle and increasing the magnitude of the shear gradients near the boundaries, as seen in the curves for p x < 0 in fig. 5(a) and (b). The addition of a positive pressure gradient would, in a Newtonian fluid, lead to a parabolic velocity profile with a minimum velocity in the centre of the channel, in opposition to the activity-induced flow shown in fig. 3(b). Such a pressure gradient acts to negate the positive flow velocity associated with an antisymmetric director angle solution, thereby reducing the alignment with the Leslie angle and decreasing the magnitude of the shear gradients near the boundaries, as seen in the curves for p x > 0 in fig. 5(a) and (b). Further increases in the positive pressure gradient lead to reverse flow near the channel boundaries and, at sufficiently high p x values, force the active nematic to flow in the negative x-direction throughout the channel. For both p x < 0 and p x > 0, the pressure gradient-induced and activity-induced velocity profiles share the same symmetry about the centre of the channel, so the spatial symmetry of the final state is unaffected. The opposite antisymmetric θ(z) solutions (i.e. the negative of the solutions in fig. 5(a) and (b)) are not presented in fig. 5. However, as mentioned in sect. 2, we can establish their behaviour by recognising that the governing equations (2) and (3) are unchanged under the We can see in fig. 6 that a fixed pressure gradient (e.g. p x = −2.5 × 10 4 Pa m −1 ) will enhance one version of the antisymmetric θ(z) solution (as shown in fig. 5(a)) while the opposite (or negative) antisymmetric θ(z) will be diminished. This effect can be seen in the breaking of the φ o → −φ o symmetry of the bifurcation diagram resulting in a perturbed pitchfork bifurcation structure for the stable branches in the φ e = 0 plane, as shown in fig. 6.
For symmetric director angle profiles, i.e. θ(z) shown in fig. 3(c) and the equivalent opposite state −θ(z), the introduction of a pressure gradient will increase the velocity in one half of the channel while decreasing it in the other. Consequently, both positive and negative pressure gradients will induce a flow that will lead to an asymmetry with respect to the centre of the channel, as shown in fig. 5(c) and (d). As a result, the equilibrium branches move out of the plane φ o = 0, illustrated by the dashed curve in fig. 6. However, even when a pressure gradient is applied, the θ(z) profile in fig. 3(c) and the equivalent opposite state −θ(z) retain a symmetry with each other, namely θ(z) = −θ(d − z). Therefore, although the symmetry of the individual states is broken, the bifurcation diagram in fig. 6 retains the φ e → −φ e symmetry that occurred in the case of zero pressure gradient.

Contractile active nematic
We now turn our attention to contractile agents, for which ζ > 0. A linear analysis shows that, in the absence of a pressure gradient, the trivial state is stable for all ζ > 0 [3], although other solutions of eqs. (2) and (3) do exist. Solution profiles for contractile active nematics are characterised by director angle configurations that exhibit large gradients close to the boundaries or the centre of the channel ( fig. 7(a), (c), respectively). These are associated with localised "jets" in the velocity ( fig. 7(b), (d)) which increase in magnitude and become increasingly sharp as the activity increases. From fig. 7(a), (c) we also see that our solutions match the predicted behaviour from sect. 3 for high values of the activity parameter, namely that the director angle may take the values θ = π/2 or θ * ≈ 1.212 rad in the bulk of the channel. While sect. 3 used an approximate form of the director angle and here we have obtained solutions numerically, we see that our numerically calculated value θ * ≈ 1.199 rad is very close to that determined analytically. The solution in fig. 7(a), (b), we believe, has been observed previously in the paper of Marenduzzo et al. [37], albeit for lower values of the activity so that the asymptotic value θ * was less readily observable. Furthermore, there have been no previous investigations into the bifurcation and stability of either the symmetric or anti-symmetric contractile director angle solutions, even for this zero pressure case. As with extensile agents, the symmetry θ(z) → −θ(z), v(z) → −v(z) means that there are also states of opposite sign in addition to the solutions in fig. 7.
In fig. 8 we show the reference zero pressure bifurcation diagram, observing that the non-trivial branches annhilate at fold (saddle-node) bifurcations. When we include the trivial state, there are up to nine possible equilibria for each activity parameter. In fig. 7 we have shown solutions on only two of the eight non-trivial solution branches: stable symmetric director profiles in fig. 7(a), (b), and antisymmetric director solutions in fig. 7(c), (d) that are un-  stable up to a single perturbation mode. The stable symmetric branches, and the equivalent solutions of opposite sign, are the non-trivial solid curves in the plane φ o = 0 in fig. 8, while the unstable antisymmetric branches are the dashed curves in the plane φ e = 0. These are stable to antisymmetric director angle perturbations but unsta-ble to symmetric director angle perturbations. The presence of non-trivial stable solutions, and the fold bifurcation points, is particularly interesting as it demonstrates multistability of solutions in the contractile situation. When a pressure gradient is introduced, the solutions in fig. 7 adapt in a similar way to the extensile case, as seen in fig. 9. For equilibria with antisymmetric velocity profiles (e.g. fig. 7(a), (b)), the introduction of a pressure gradient leads to asymmetry in both the director angle and velocity ( fig. 9(a), (b)). The positive and negative solution branches in the plane φ o = 0 in fig. 8 both adapt in the same way and are promoted to occur at lower values of the activity parameter. For a negative pressure gradient, equilibria for which the velocity is symmetric and positive ( fig. 7(c), (d)) will be enhanced and occur at smaller values of activity. Conversely, the pressure gradient-induced flow will retard negative velocity solutions meaning they can occur only for larger activities. This break in the φ o → −φ o symmetry is observed in the plane φ e = 0 in fig. 10.
As was shown in the extensile case, we have seen that for contractile active nematics, certain solutions are promoted. In other words, they exist over a wider range of activities, while the span of other solutions is reduced.
One key result here is that the stable non-trivial state in the contractile case is promoted through the application of a pressure gradient. The ranges of activity for which each solution exists, as the pressure gradient is changed, can now be investigated to understand the existence of each solution branch in the (ζ, p x )-parameter plane.

Two-parameter continuation
We can summarise the effect of varying the activity parameter and applied pressure gradient on the various equi-libria, for both extensile and contractile active nematics, by considering the critical bifurcation points in the (ζ, p x )space, i.e. the location of fold and pitchfork bifurcations. By constructing a (ζ, p x )-space diagram for this multistable system, we are able to determine which of the various solutions found in the previous section exist over a range of pressure gradients and activity strengths. There are seven possible critical bifurcation points: the single fold/pitchfork point in fig. 6 fig. 9 and their negative versions. However, for clarity we do not plot the locations of the two fold points of the branches of unstable solutions in the φ e = 0 plane in fig. 10 because, being unstable, these solutions are unlikely to be observed in reality. Of the remaining five critical points, the two pairs of fold points occur at the same value of the activity parameter and so the locations of the critical points will lead to only three loci in the (ζ, p x )-space. We plot these three loci in fig. 11.
Curve (A) in fig. 11(a) indicates the location of the fold bifurcations where the four solution branches meet in fig. 6 for p x = 0 (or the pitchfork bifurcation point in fig. 4 when p x = 0). The fold bifurcation points in the plane φ o = 0 in fig. 4 for p x = 0 correspond to symmetric director angle profiles. The position of these bifurcation points as p x varies, coinciding with the director angle losing symmetry, is represented by curve (B). Similarly, curve (C) corresponds to the fold bifurcation points lying in the plane φ o = 0 in fig. 8. In the region between curves (B) and (C), the only stable equilibria will therefore be on the stable branch of the perturbed pitchfork bifurcation that exists for all values of ζ, as seen in figs. 6 and 10. This region in the (ζ, p x )-space is characterised by solutions of relatively low director distortion and velocity, for which the activity parameter is too small in magnitude to dominate either the elastic or the pressure gradient effects. In fact, when the pressure gradient is absent, the solutions revert to the trivial case (θ, v) ≡ (0, 0). Therefore, in fig. 11 we denote the low distortion regime as the "trivial state", though it is generally more accurate to describe it as the perturbed trivial state. For (ζ, p x ) to the left of curve (A) there are four stable equilibria, characterised by antisymmetric flow-alignment close to ±θ L and asymmetric flowalignment close to ±(π − θ L ) rad. (The asymmetry is the result of a pressure gradient on a symmetric director angle structure similar to the behaviour in figs. 5(c) and 9(a).) Between curves (A) and (B) for extensile active nematics, and to the right of (C) for a contractile agent, there are two stable, asymmetric director angle solutions as well as the trivial state. These correspond, respectively, to an alignment close to ±(π−θ L ) rad in much of the channel for extensile activities or at the angle approximated by ±θ * as found in sect. 3 for contractile activities. The process of transition between these states in the connected region surrounding the shaded area between curves (B) and (C) is now investigated.
The extent of the region between (B) and (C), where the perturbed trivial state is the only equilibrium, decreases as the magnitude of the pressure gradient increases, with the two curves meeting at approximately ζ = −2.05 Pa. Figure 11(b) focuses on the region close to this activity parameter for p x < 0 highlighted by the grey box in fig. 11(a) and demonstrates a swallowtail catastrophe [51]. Curves (B) and (C) are symmetric with respect to the pressure gradient, so a similar swallowtail feature occurs for p x > 0. The process is expanded upon in fig. 12 where we show bifurcation plots in the (ζ, φ e )-plane as |p x | increases. The φ o measure for equilibria θ(z) will be non-zero for these solutions, but the transition can be illustrated by restricting the attention to φ e . The horizontal dotted lines in fig. 11(b) correspond to the pressure gradients used in fig. 12. The perturbed trivial branch is not shown in fig. 12, but it does also exist for this range of pressure gradients and activities.
For p x = −7.7 × 10 5 Pa m −1 in fig. 12(a), there is a small interval of activity over which the trivial state is the only equilibrium, between the left-hand fold point (curve (B) in fig. 11) and the right-hand fold (curve (C) in fig. 11). In fig. 12(b), the right-hand fold has transformed at a cusp catastrophe and this branch now has three fold points associated with it. In fig. 12 fig. 12(c), the range of activity for these extra states has widened. Also, there are no longer any activities for which the trivial state is the unique equilibrium. Between fig. 12(c) and fig. 12(d), at p x = −7.76 × 10 5 Pa m −1 , the two uppermost fold points have merged creating a stable branch of solutions that exists for all activity parameter values, and with two unstable branches now linked through two fold points and a stable branch. Finally, for large enough |p x |, in fig. 12(e), the final two fold points have annihilated leaving one stable and one unstable branch, which both exist for all activity parameter values. It is the formation of this stable branch that allows the continuous transition between solutions with alignment close to ±(π − θ L ) rad for extensile activity, to solutions with alignment close to ±θ * for contractile flow.
In figs. 13 and 14, we consider this transition between ±(π − θ L ) rad, and ±θ * alignments for a value of |p x | greater than the value considered in fig. 12(e). Clearly pressure gradient-induced flow dominates for small activity parameter values ζ. However, activity still dominates the pressure gradient effects for sufficiently large |ζ| values. Figure 14(a) shows the transformation from a solution in which there is alignment with (π − θ L ) rad (e.g., ζ = −200 Pa), to a solution in which there is alignment with θ * (ζ = 200 Pa). The corresponding flow profiles are plotted in fig. 14(b), (c). Figure 13 shows, therefore, that for sufficiently high pressure gradients it is possible to transition from a director state where alignment with θ * dominates to a state where alignment with the Leslie angle dominates through changes in activity.  for θ(z). The five pressure gradients correspond to those indicated in fig. 11(b). Solid curves are stable equilibria solution branches, while unstable equilibria solution branches with one positive eigenvalue are dashed curves. Note that, for these values of the activity parameter, the trivial solution is also stable.  . 13. Stable equilibrium branch for pressure gradient px = −10 6 Pa m −1 , using the measure φe for θ(z). As the activity parameter ζ increases, the asymmetric director angle solutions transition from activity-induced alignment at (π − θ L) rad in much of the channel to alignment at θ * .
It is worth noting that the special case ζ = 0 Pa in fig. 14 does not represent the classic Poiseuille flow for a nematic with a parabolic velocity profile. Poiseuille flow would be the result of applying a pressure gradient to the trivial state in the absence of activity, leading to a symmetric velocity and an antisymmetric director angle solution which vanishes in the centre of the channel. The states examined in figs. 13 and 14 are inherently asymmetric with the director angles non-zero except at the boundaries.

Conclusions
We have analysed the effects of introducing an applied pressure gradient on the flow in a channel containing an active nematic. As expected, in the limit of large pressure gradients, the nonlinear equation for the director angle has equilibrium solutions that demonstrate flow-alignment associated with the Leslie angle, θ L , with an associated classic Poiseuille parabolic flow. However, less intuitively, in the limit of large activity, we found analytically that there were other possible director orientations in the bulk of the channel, specifically where θ = 0, π/2 or θ * ≈ 1.212 rad.
Considering steady-state solutions of the full nonlinear system, we then determined the effect of an applied pressure gradient on the various symmetric and antisymmetric director angle solutions. For an extensile active nematic, the application of a pressure gradient leads to a perturbed pitchfork bifurcation in which the pressure introduces elements of directional bias in the system. We observed that the antisymmetric director angle solutions coinciding with zero pressure gradient remained antisymmetric when a pressure gradient is applied, whereas the corresponding symmetric director structures became asymmetric. For a contractile active nematic, we found nine possible solution branches, including the trivial state, although only two of the non-trivial branches contain stable solutions. The pressure gradient has a similar effect on the symmetry of the solutions as in the extensile case. Crucially, it promotes some non-trivial solutions which then occur over larger ranges of activity. We have therefore discovered a rich bifurcation structure in this model of pressure-driven flow of an active nematic within a channel. We have found that director angle solutions that align at a non-intuitive angle θ * can, under sufficiently high pressure gradients, be induced to occur for a wide range of activity values, extensile or contractile. In fact, beyond a critical value of |p x |, these solutions form a continuous transition to flow-aligning solutions that extends over all activity values. The presence of two stable states, for instance the trivial state, θ = 0, and the θ = θ * state, is important when considering the rheology of active nematics. For instance, in microfluidics applications the effective viscosities of these two states, g(0) = η 1 and g(θ * ) ≈ η 2 , can be very different (in 5CB, they vary by a factor of five, η 1 = 0.0204 Pa s and η 1 = 0.1052 Pa s) and so the volumetric flow rate will also be very different. One way of manipulating the flow in a channel will be to influence the stable director orientation of the active nematic. Given the current interest in designing microfluidic systems in which active fluids play a role in driving flow and guiding behaviour at channel junctions, we would hope that the state diagram in fig. 12 would be useful in using pressure gradients to design specific microfluidic devices.

Author contribution statement
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A. Large activity asymptotics for non-flow-aligning nematics
There is one particular parameter regime for which we are able to make progress by assuming a local form of eq. (22). This parameter regime exists if the nematic viscosities are such that m(θ)/g(θ) = λ, a non-zero constant. This is the situation, or an approximation to the situation, for nonflow-aligning nematics, i.e. when m(θ) = 0 does not have a solution, and m(θ) is of one sign. We now consider the local form of eq. (22) by taking C to be a constant, which we set to be I, so that Through solving the local problem in eq. (A.1) directly, we find that the function β min (I) is monotonically increasing with β min → ∞ as I → I * (see Schaaf [52]). In this limit (θ → θ * , I → I * , β min → ∞) we therefore have θ ξ → 0 and so θ * must satisfy H(θ * , 0) = H(0, 0), or equivalently Given the definition of I * , namely I * = sin θ * cos θ * , eq. (A.3) is equivalent to tan θ * = 2θ * . This equation may be solved numerically to give θ * ≈ 1.166 rad. This is the predicted director angle in the bulk of the region, away from any boundary or reorientation regions, for symmetric θ solutions in the limit of large positive (contractile) activity parameter for a non-flow-aligning nematic. We have also determined numerically the director profile for a contractile active nematic in the situation when m(θ)/g(θ) is a constant. When we choose parameter values K 1 = K 3 = 10 −11 N, η 1 = η 2 = 0.0204 Pa s, η 12 = 0, γ 1 = 0.0777 Pa s and d = 10 μm, the behaviour of the director angle in the centre of the cell, for increasing activity parameter, is similar to that in the flow-aligning case, with the asymptotic value exactly equal to the value determined analytically, namely θ * ≈ 1.166 rad.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.