Active nematics in corrugated channels

Abstract

Active nematic fluids exhibit complex dynamics in both bulk and in simple confining geometries. However, complex confining geometries could have substantial impact on active spontaneous flows. Using multiparticle collision dynamics simulations adapted for active nematic particles, we study the dynamic behaviour of an active nematic fluid confined in a corrugated channel. The transition from a quiescent state to a spontaneous flow state occurs from a weak swirling flow to a strong coherent flow due to the presence of curved-wall induced active flows. We show that the active nematic fluid flows in corrugated channels can be understood in two different ways: (i) as the result of an early or delayed flow transition when compared with that in a flat-walled channel of appropriate width and (ii) boundary-induced active flows in the corrugations providing an effective slip velocity to the coherent flows in the bulk. Thus, our work illustrates the crucial role of corrugations of the confining boundary in dictating the flow transition and flow states of active fluids.

---

1 Introduction

Active and living systems are out of equilibrium due to continuous influx of energy at microscopic scales. Common examples of these systems are suspensions of cytoskeletal filaments powered by molecular motors,^1–4 bacterial suspensions,^5,6 cellular layers^7–9 and Janus catalysts.^10–12 Continuous energy injection not only leads to the motion of individual active particles but also results in a rather exciting phenomenon of collective motion. Collective dynamics in dense suspensions of active particles usually consist of jets and circulations – a state commonly called active turbulence.^13–15 Several experimental realizations of active turbulence have been demonstrated in the literature.^2,5,16–19

Active nematics are model fluids that capture the dynamics of some of the above mentioned systems.^13,20–22 By incorporating the active stress generated by microscopic active entities, active nematics build upon the theory of passive nematic liquid crystals.^23–25 Hence, mechanisms such as generation and dynamics of topological defects that sustain active turbulence, are inherent features of active nematic fluids.^13,26–29 Active stress generates hydrodynamic instabilities overcoming the elastic, viscous and frictional forces in the active fluid. In a confined active nematic fluid, the competition between activity and the opposing forces gives rise to a spontaneous flow transition at a critical activity.^3,21,30–32

Confinements play a major role in dictating the dynamic state of active systems, much different from those of driven systems. For example, when active nematics are confined in microchannels of square cross section, they exhibit a multitude of states, such as a state of no fluid flow, streaming states like unidirectional, oscillatory or double helix flows and swirling states including vortex lattices or turbulent flows. All these depend on varying a single parameter, namely the ratio of the channel width to the active length scale.^21,33,34 Despite the crucial importance of boundaries in dictating the dynamics of active fluids, most investigations so far have considered only channels with flat walls,^{3,31,32,35–38} cylindrical channels,³⁹ or annular rings.^33,40–42 The effect of non-uniformity of the channel walls on the dynamics of active fluids has not been addressed in the literature yet.

However, some of the previous investigations have hinted at the effect of non-uniform boundary walls on the dynamics of active nematic fluids. In the experiments of Wu et al.,³³ asymmetric corrugations of the channel walls were used to direct streaming flows generated by mixtures of microtubule bundles and kinesin molecular motors. When the active suspension was confined in a toroid of square cross section it produced a streaming flow in the channel, with equal probability of streaming in clockwise and counter clockwise directions. On the other hand, when the channel walls of the toroid were modified by providing saw-tooth like notches, it was possible to control the directionality of the streaming flows. Such a control on the directionality of active fluid flows is essential for any future engineering applications. While the exact mechanism that leads to directionality is not clear, it was shown that the directionality of the flow in the toroid was accompanied by swirls in the undulated sections of the channel wall. Being at low Reynolds number, the presence of swirling flows in the undulations of the channel is also intriguing. Asymmetric structures driving directionality have also been demonstrated in other active systems.⁴³ Another feature the curved boundaries provide is their capability to drive active fluid flows. The boundary curvature driven flows arise since active entities preferentially orient at the solid wall, as demonstrated using theoretical,^44,45 computational⁴⁶ and experimental⁴⁷ tools. However, the role of these boundary-induced active flows in dictating the flow transition, the directionality of flow or the flow state itself has not yet been analyzed.

In this work, we computationally investigate the flow generated by active nematic fluids when confined in corrugated channels. Since the spontaneous flow transition in an active nematic fluid is well studied theoretically,³⁰ experimentally^3,32 and numerically,^13,21 we primarily focus on the effect of corrugations on the flow transition from a no-flow state to a flow state of an active nematic fluid. The resulting flow state immediately after the transition is, typically, a streaming flow state and analysis of this flow further allows us to isolate the role of corrugations and boundary-induced active flows. The investigation is performed using a numerical framework based on multiparticle collision dynamics.⁴⁸

In Section 2, we discuss the system and the multiparticle collision dynamics algorithm for active nematic fluids. We first demonstrate the validity of the algorithm illustrating the spontaneous flow transition in flat-walled channels (Section 3). This is followed by an analysis of flow states generated in corrugated channels and contrasted with those in a flat-walled channel. Further, the role of geometry in establishing the flow transition and streaming-flow state is characterised and generalised. Finally, we provide a slip-based mechanism to understand the spontaneous flow transition in a corrugated channel and conclude the work (Section 4).

2 Methodology

2.1 Active nematic multiparticle collision dynamics (AN-MPCD)

The active nematic fluid confined in a corrugated channel is simulated using the active-nematic multiparticle collision dynamics (AN-MPCD) algorithm. The AN-MPCD algorithm is a particle based, mesoscopic algorithm proven to be useful to simulate various soft matter systems. By assigning a nematic stress to the particles, the AN-MPCD algorithm simulates the dynamics of active nematic fluids.^48–51 The AN-MPCD algorithm discretizes the continuum into N point-like particles, each with mass m, position r_i, velocity v_i and orientation u_i. The algorithm comprises two steps: streaming and collision. During streaming, each particle moves ballistically for a time δt to a new position

r_i(t + δt) = r_i(t) + v_i(t)δt. (1)

The collisions are accomplished by partitioning the domain into cubic cells. These collisions are governed by mesoscopic collision operator Ξ_i, which is stochastic and ensures the conservation of momentum within each cell. This momentum transfer among particles thereby updates the velocities according to

v_i(t + δt) = v(t) + Ξ_i (2)

where v(t) = 〈v_i〉 is the center of mass velocity of the cell, calculated at position r in the cell. Since all particles in the cell are assumed to have identical mass m, the center of mass velocity of the particles in the cell remains unchanged during the collision process. Further, the conservation of energy and angular momentum is also achieved by the constrained stochastic exchange of particle velocities. The hydrodynamic velocity field can be obtained as the average velocity v(t) at position r in the cell.

The collision operator injects energy into the system, corresponding to an extensile or contractile force dipole density. This is achieved by incorporating the local active stress proportional to the nematic tensor order parameter Q, such that the force dipole coaligns with the local director field n of the cell. Thus, the collision operator is a linear combination of passive and active contributions:

where Ξ^N_i is a nematic multi-particle collision operator (see Appendix A) which corresponds to the passive contribution^48,52 arising from the orientational order of nematic particles. The active contribution comprises two terms: (i) individual impulses (per unit mass) (aδt/m)κ_in for each particle i, representing the active force driving a change in momentum at each time step, and (ii) a term to ensure local conservation of momentum −(aδt/m)〈κ_j〉n, which is averaged over a cell. The first term (aδt/m)κ_in is composed of three factors: (i) a, which represents the local active dipole strength in the cell given by a = Ñ × α, where α is the activity of each particle and Ñ is the number of particles in a cell, (ii) κ_in = ±n, sets the direction of the active force acting on the particle i, and (iii) δt/m, ensures that a has units of force. As stated above, a is directly proportional to the number of local active agents in the cell, and positive (negative) particle activity α leads to extensile (contractile) active nematics. The factor κ_i is a parallel/antiparallel coefficient. For particles that are above the plane defined by the center of mass r and the director field n, κ_i = +1, indicating that the particle is driven “forward”, and particles below the plane are kicked “backward” (κ_i = −1).

The collision operator also updates particle orientation u_i, thus computing the dynamics of the nematic tensor order parameter, defined as Q = 〈2u_iu_i − 1̂〉. Here 1̂ denotes the identity matrix. The tensor order parameter Q quantifies the extent and direction of orientational order through its largest eigenvalue S and the corresponding eigenvector n in the cell. The orientations are updated based on the local equilibrium distribution of the orientation field

u_i(t + δt) = n(t) + η_i (4)

where the noise η_i is drawn from the Maier–Saupe distribution ∼exp((US/k_BT)[u_i·n]²). The mean-field interaction constant U and inverse thermal energy 1/k_BT govern the width of the distribution around n(t). Consequently, when US/k_BT is small, all orientations are equally probable, leading to an isotropic phase. For large US/k_BT, the distribution sharply centers about n resulting in a nematic phase. The gradients in velocity are coupled to shear alignment through discretised Jeffery's equation (see Appendix A) of a slender rod with tumbling parameter λ.⁵² The value of the tumbling parameter sets the nematic fluid in either the shear-aligning or flow-tumbling regimes. Previous investigations^48,52 have shown that when λ < 1 the nematogens tumble with the flow. A relaxation parameter, χ, allowing averaging of Jeffery's equation over a small number of time steps, and a viscous rotation coefficient, γ_R, are incorporated in Ξ^N_i (see Appendix A) for the coupling of the velocity field to director dynamics. The coupling in necessary to address the backflow effects produced due to the torques on the MPCD particles. Thus, the fluid velocity and orientation of the nematic fluid are two-way coupled to produce active nematohydrodynamics.

2.2 Simulation details

We consider flow aligning nematics with extensile activity (α > 0) confined in a corrugated channel.⁵³ A schematic of the corrugated channel used for the simulations is shown in Fig. 1. Cartesian coordinates are adopted − x along the length of the channel and y in the perpendicular direction, i.e., across the channel, and the channel centerline is at y = 0. The profile of the corrugated walls in the channel is represented by cosine functions and the equations describing the top and bottom wall are y_top(x) = (W/2) + A cos(2πx/Λ) and y_bottom(x) = (−W/2) − A cos(2πx/Λ), respectively. Here, W is the mean width of the channel and A and Λ are the amplitude and wavelength of the corrugations respectively. Since the width varies along the channel, the minimum width, W_min, is defined as the distance between the trough of the top wall, W/2 − A, and the crest of the bottom wall, −W/2 + A, i.e., W_min = W − 2A. Similarly the maximum width in the corrugated channel is W_max = W + 2A.

Fig. 1 Schematic diagram of the corrugated channel used in the simulations: A and Λ are the amplitude and wavelength of the corrugations, respectively, L is the length, W is the mean width, W_min is the minimum width and W_max is the maximum width of the channel. The boxed equations in the figure describe the geometry of the top and bottom wall of the channel.

As discussed in Section 2.1, the total activity of the fluid in the channel is proportional to the number of active particles in the system. Therefore, for a fair comparison, (i) it is necessary for the volume of the fluid in the corrugated and flat channel to be identical and (ii) the fluid is of the same density and viscosity in both the cases. In this work, two-dimensional simulations are considered; therefore, the area of the corrugated channel occupied by the active fluid (volume per unit length) is given by

which is equal to the area of the flat-walled channel. This equivalence ensures that the number of active particles used to discretize the continuum in both flat and corrugated channels are equal. Secondly, the density and viscosity of the fluid in the AN-MPCD algorithm are determined by the average particle density of the cell.⁵⁴ Therefore, the average particle density of the cell is set to 〈N〉 = 20 in both cases, to ensure the same fluid properties in all channel geometries.

The results are reported in MPCD units of cell size l = 1, particle mass m = 1 and thermal energy k_BT = 1, which leads to units of time . The AN-MPCD particles are initialised with random velocities and oriented along the y axis. In the simulations, the mean-field interaction constant U = 10 is set to achieve the nematic phase.⁴⁸ The tumbling parameter λ = 2 is chosen for the active nematic fluid to be in the shear-aligning regime.⁴⁸ The heuristic coupling coefficient χ is fixed at 0.5 and the viscous rotation coefficient is fixed at γ_R = 0.01. The time step is set to δt = 0.1 and simulations are run for a total of 6 × 10⁵ time steps which includes a warm-up phase of 3 × 10⁵ time steps. The channels investigated are of the following dimensions: (i) a flat-walled channel (A = 0) with length L = 100 and width W = 20 and (ii) a corrugated channel with wavelength Λ = 20, amplitude A = 3 and the same length and width as in (i) unless specified otherwise. Periodic boundary conditions are imposed along the x direction for the flow velocity and orientation tensor (see Appendix B.3.1). The walls of the channel are impermeable with no-slip boundary conditions using phantom particles (see Appendix B.3.2). A strong homeotropic boundary condition is enforced for the director field ensuring that the director field is oriented perpendicular to the channel walls (see Appendix B.2).

2.3 Corrugated channel boundaries

The boundary conditions that simulate the active nematic fluid in a corrugated channel as shown in Fig. 1 comprise two parts: (i) the wall surface and (ii) the boundary rules for the transformation of the velocity during the collision event. The channel walls are sinusoid functions as described in Appendix B. The rules for velocity transformation during the collision event with a wall are bounce-back. The transformation of the orientation of the particle involves multiple operations to obtain the desired anchoring on the curved wall. The details of implementation of these boundary conditions are provided in Appendix B.3, and is closely adapted from Wamsler et al.⁵³

3 Results and discussion

We now present the results of two dimensional simulations of an active nematic fluid confined in flat and corrugated channels. We first investigate the active nematic fluid flow in a flat-walled channel by systematically varying the activity α and channel width W. The results confirm that AN-MPCD simulations are in agreement with previous studies^30,34,55 (Sections 3.1 and 3.2). Subsequently, by replacing the flat walls with corrugated walls, we identify the differences in the fluid flow behaviour and the underlying physical mechanisms involved (Sections 3.3 and 3.4). The amplitude, A, and the wavelength, Λ, of the corrugations are altered to understand the effect of geometrical parameters. Lastly, the influence of circulations trapped in the corrugations is discussed along with the results of a linear stability analysis where the effect of the swirls trapped in the corrugations of the channel is interpreted as an effective slip velocity (Section 3.5).

3.1 Spontaneous flow transition in a flat-walled channel

Theoretical^30,34,55 and experimental³² investigations have shown that a confined active nematic fluid undergoes a spontaneous flow transition at a threshold critical activity,³⁰ herein denoted as α_c. This flow transition can be characterised by a non-flowing state when α < α_c and a spontaneous flow state when α > α_c. To test this, using the AN-MPCD algorithm, we simulate an active nematic fluid confined in a flat-walled channel with systematic increments in the value of activity. Both unidirectional and bidirectional flows have been discussed in the literature following a spontaneous flow transition;^{13,30,32,34,55–57} however, in this manuscript, we focus on unidirectional flows and the scaling^30,55 relationships that connect the strength of the flow, activity and confinement length scale. The flow in the channel is characterized by the root-mean-squared velocity, , where v is the fluid velocity in the cell. The calculation of V_rms is a multi-step process; first the squared fluid velocity v² is computed within each cell and spatially averaged to obtain 〈v²〉(t). Subsequently, a steady-state is identified and the time average of is computed. It is evident that a no-flow state (negligible V_rms) is observed for α < α_c and a flow state is observed for α > α_c (Fig. 2(a)). Above α_c, V_rms continues to increase as a function of activity. The sudden change in V_rms is indicative of phase transition-like behavior. We determine the critical activity by curve fitting the V_rms data of both no-flow and flow states to the mathematical form (dashed curves in Fig. 2(a)). The critical activity α_c of the system is identified from the intersection point (diamond in Fig. 2(a)) of the fitted curves.

Fig. 2 Spontaneous flow transition in a flat-walled channel characterized by change in root-mean-squared (RMS) velocity V_rms as a function of (a) activity α for W = 20 and (b) channel width W for various activities. (c) The critical width W_c plotted as a function of activity for a flat-walled channel. (d) Plot of V_rms/α^1/2 as a function of the reduced distance from the critical point (W − W_c)/W_c for various activities, collapsing the data reported in (b) into a single curve.

The critical activity for spontaneous flow transition of confined active nematics can be calculated from a linear stability analysis of the governing continuum equations⁵⁵ to be of the form

where, γ is the rotational viscosity, K is Frank's elasticity constant, μ is the fluid viscosity and [small klambda, Greek, circumflex] is the flow aligning parameter. We use [small alpha, Greek, circumflex]_c and [small klambda, Greek, circumflex] for critical activity and flow aligning parameter in continuum theory to differentiate the parameters from their AN-MPCD counterparts. From eqn (6), it is apparent that apart from the fluid properties, critical activity is also a function of the channel width, W. This indicates that the presence of a critical width W_c, is also a function of activity.

To quantify spontaneous flow transition using W_c, we systematically increase the width of the channel for various activities (Fig. 2(b)), where V_rms is plotted as a function of channel width W. The active nematic fluid undergoes the transition from no-flow to spontaneous flow as the channel width is increased beyond W_c for a specified activity. The critical width is computed using a similar methodology to determining critical activity by curve fitting the data using and extracting W_c. The critical width W_c decreases as activity increases (Fig. 2(c)). This analysis confirms that for a flat-walled channel, the critical width scales as , consistent with eqn (6) and previous investigations.³⁰

Further, it can be seen that the dependence of V_rms on W is similar for different activities (Fig. 2(b)). Therefore, we collapse the data by plotting V_rms/α^1/2 as a function of (W − W_c)/W_c (Fig. 2(d)). This collapse is explained through a scaling analysis from the linear stability results. At steady state, the active stress (∼αSθ) and the viscous stress (∼μv/W) balance

where θ is the characteristic angle of the director field. Upon spontaneous flow transition, we can consider θ to have the form

Thus, the angle scales with the reduced width, which acts as a non-dimensional function of the channel width akin to a second order phase transition. Since W ≈ W_c near the transition, combining eqn (6)–(8) gives the velocity of the active nematic fluid

explaining the collapse of the simulation results shown in Fig. 2(d). Therefore, the AN-MPCD simulations of a fluid confined within a flat-walled channel are in agreement with the investigations of the existing literature.

3.2 Coherent flows in a flat-walled channel

When confined, an active nematic fluid undergoes a spontaneous flow transition, and a coherent flow is usually observed immediately above the transition point.^31,36,58,59 Coherent flow is characterized by velocity fields that are dominant along the length of the channel v_x ≈ |v|. We now discuss the underlying physical mechanism associated with flow transition to a coherent flow in a flat-walled channel (Fig. 3). The flow field is plotted by superimposing streamlines on fluid vorticity, ω = (∇ × v)_z, where the colors represent the clockwise (blue) and anticlockwise (red) rotation of the fluid elements. The local director field n is shown by a black dashed line color shaded by the scalar order parameter, S. Active force, f_a, produced by the nematic fluid is calculated as f_a ∝ ∇·Q and is shown with black arrows depicting the direction and color shading depicting the magnitude of the force produced.

Fig. 3 Snapshots of the flow field (left column), director field (middle column) and active forcing (right column) in flat and corrugated channels at lower (α < α_c) and higher (α > α_c) activity. In the left column, the flow field is shown, which is represented by streamlines superimposed on the vorticity of the fluid. The colorbar is shown at the bottom with blue and red representing clockwise and anticlockwise rotation respectively. In the middle column, the director field is shown by a black dashed line color shaded by the scalar order parameter. In the right column, the active force is illustrated by black arrows which are color coded by the magnitude of the force. Figures (a)–(f) correspond to low activity and figures (g)–(l) correspond to higher activity at which coherent flows prevail. All channels have a width of W = 20 and corrugated channels have an amplitude of A = 3 and wavelength Λ = 20. For the flat-walled channel, α_c = 0.015; the cases of low and high activity correspond to α = 0.009 and α = 0.017 respectively. For the corrugated channel, α_c = 0.019; the cases of low and high activity correspond to α = 0.017 and α = 0.023 respectively.

For a lower activity (α < α_c), in a flat-walled channel a no-flow state is observed in the flow domain (Fig. 3(a)), with the director field aligned in the nematic phase (Fig. 3(b)). The active forcing produced by the fluid is insignificant (Fig. 3(c)). However, at a higher activity (α > α_c), a bend-like deformation of the director field is observed as shown in Fig. 3(h), in which the thick dashed line is a visual representation of the bend deformation and the arrow represents the direction of fluid flow. The active force produced by the fluid corresponding to the bend deformation is maximum at the center of the channel and is oriented along the x axis (Fig. 3(i)). Thus, the active force drives the flow along the length of the channel producing coherent flows (Fig. 3(g)).

3.3 Coherent flows in a corrugated channel

In the case of a corrugated channel, when the activity α is lower than a critical value α_c, swirling flows constituted by closed streamlines are observed. This is illustrated in Fig. 3(d) for an amplitude of A = 3 and wavelength Λ = 20, demonstrating a significant difference from the no-flow state observed in a flat-walled channel when α < α_c.

The origin of swirling flows can be explained by observing the director field (Fig. 3(e)). Due to the strong homeotropic anchoring boundary condition, the director field aligns perpendicular to the wall. However, due to the reflection symmetry about the channel centerline, the director field within the bulk of the channel remains oriented perpendicular to the centerline. Therefore, the nematic elasticity induces a bend-like deformation close to the channel wall which is visually depicted by thick dashed lines in Fig. 4(a). As a consequence of this bend, the fluid experiences an active force near the walls (Fig. 3(f)), where the magnitude of force is maximum at the bends. This active force drives the flow towards the centre of the corrugation, as indicated by the arrows in Fig. 4(a). Since the corrugation is symmetric with respect to the y-axis, two bends are observed near the wall, facing each other. As a consequence, the flow converges at the center of the corrugations and is subsequently driven towards the bulk of the channel. This dynamics leads to the formation of counter-rotating swirls. Due to the symmetry of the channel, a pair of bends are observed at both the top and bottom walls, resulting in four swirls within the fluid (Fig. 3(d)). The swirling flows prior to spontaneous flow transitions have been analyzed by Zumdieck et al.⁴⁴ Calculations in the limit A/Λ ≪ 1 and W ≪ W_c ≪ Λ showed that the velocity of the swirls depends on the geometry of the channel as .

Fig. 4 An enlarged view of the director field in the corrugated channel: (a) below critical value α < α_c from Fig. 3(e) and (b) above the critical value α > α_c from Fig. 3(k). The red dashed lines represent the bend like deformation and the blue arrows point in the direction of the active forcing and the subsequent fluid flow.

At higher activities (α > α_c) in the corrugated channels, coherent fluid flow is obtained alongside the swirls near the walls (Fig. 3(j)). From Fig. 3(k) we observe that a bend like deformation can be noticed in the bulk of the channel (Fig. 4(b)), which results in coherent flow in the channel. This is in contrast to lower activities (Fig. 4(a)). The combination of strong anchoring boundary condition and nematic elasticity induces a bend-like deformation in the director field close to the channel wall with opposite polarity to the bend in the bulk (Fig. 4(b)). The bend in the director field close to the wall generates an active force, driving the fluid towards the centre of undulation and then towards the bulk of the fluid resulting in swirling flows. From Fig. 3(l), it is evident that the active force corresponding to the bend in the director field close to the wall is pronounced compared to the bulk of the channel. It is these regions of high bend and strong active forcing that drive and sustain the swirling flow alongside the bulk coherent flow. The swirls generated remain trapped inside the corrugations and coexist along with coherent fluid flow along the center of the channel (Fig. 3(j)). It is interesting to note that in a flat-walled channel, for α > α_c, the fluid at the centerline of the channel experiences a large active force (Fig. 3(i)). In this case, although the active force near the walls is of comparable strength but opposite polarity, it is the active force along the centerline of the channel that dictates the direction of the coherent flow. Similarly, in the corrugated channel the active force is large around the centerline and drives the coherent flow along the length of the channel (Fig. 3(l)). However, the active force experienced by the fluid close to the walls has a higher intensity but with an opposite polarity.

The existence of coherent flows with swirls aligns with the experimental observations reported by Wu et al.³³ When the active nematic fluid (suspension of microtubule and motor protein mixture) is confined in a ratchet like geometry, coherent flows were produced but were accompanied by swirls trapped in the corrugations. This observation also emphasizes that the channel geometry can influence the flow dynamics and behavior of an active nematic fluid confined within it.

3.4 Effect of changing the amplitude and wavelength of the corrugations

In a flat-walled channel, the active nematic fluid undergoes a flow transition from a state of no-flow to a state of coherent flow. In contrast, corrugated channels exhibit counter rotating swirls prior to the flow transition. Therefore, we now investigate the effects of change in the corrugation geometry on the spontaneous flow transition of an active nematic fluid to a coherent flow state. This is achieved by varying the amplitude, A, and wavelength, Λ, of the corrugations of the channel.

The results of change in amplitude are shown in Fig. 5(a), where root-mean-squared velocity, V_rms, of the fluid is plotted as a function of activity α for various amplitudes. The active nematic fluid undergoes a flow transition at all amplitudes. However, above the critical activity, the velocity of the fluid decreases with the increase in amplitude. Further, it is also observed that larger amplitudes shift the critical activity at higher values. Note that all simulations have been carried out at a fixed mean width W, and therefore, an increase in the amplitude A reduces the minimum width, W_min, of the channel. Therefore, the active fluid encounters a stronger geometric constriction with the increase in the amplitude of the corrugation, and thus causing the shift of the flow transition to higher activities.

Fig. 5 Effect of channel geometry on spontaneous flow transition to a coherent flow state in corrugated channels. Variation in V_rms with respect to activity for corrugations of (a) different amplitudes A but fixed wavelength Λ = 20 and (d) different wavelengths Λ but fixed amplitude A = 2. In both cases, the mean width is fixed at W = 20. The channel length is chosen as L = 100 in (a) and L = 200 in (d). Figures (b) and (c) illustrate the variation in V_rms as a function of α for a fluid within the bulk (region I) and corrugations (region II) of the channel, respectively. Figures (e) and (f) depict the change in critical activity in the corrugated channel as a function of amplitude and wavelength normalised by the mean width W of the channel, respectively.

To determine the critical activity of the active nematic fluid in the corrugated channel, the corrugated channel is split into two domains: region I represents the center of the channel, |y| < (W_min/2), and region II represents the fluid trapped in the corrugations, (W_min/2) < |y| < (W/2 + A) (see the schematic in Fig. 7(a)). The root-mean-squared velocity of the fluid calculated in the central region I, V^I_rms, and in region II within the corrugations, V^II_rms, is shown as a function of activity in Fig. 5(b) and (c) respectively. As in a flat channel, V^I_rms remains small below the critical activity but increases rapidly as beyond the flow transition. The critical activity α_c corresponding to the corrugated wall is determined from Fig. 5(b). On the other hand, the fluid velocity in the corrugations, V^II_rms, continues to gradually increase with increasing activity without a discernible indication of the spontaneous flow transition (Fig. 5(c)). The velocity of the fluid trapped in the corrugations V^II_rms increases with the increase in the amplitude of the corrugations. A larger amplitude results in a stronger bend in the director field that generates a stronger flow in the corrugations. The critical activity required for flow transition determined from Fig. 5(a) is shown in Fig. 5(e) as a function of amplitude A of the corrugations. The increase in critical activity with amplitude appears to show a quadratic dependence.

To study the effect of the wavelength of the corrugations Λ on the active nematic fluid flow, we vary Λ for a fixed amplitude and mean width of the channel. A comparison of the root-mean-squared velocity V_rms of the fluid as a function of activity α for various wavelengths (Fig. 5(d)) with V_rms as a function of amplitude (Fig. 5(a)) demonstrates that the influence of wavelength in the investigated range is weaker. Changing the wavelength of the corrugations affects the velocity of the coherent flow of the active nematic fluid, with larger wavelengths resulting in larger velocity. Further, increasing the wavelength of the channel reduces α_c, the critical activity required for the spontaneous flow transition (Fig. 5(f)). The critical activity decreases with wavelength, showing a strong power law decay.

A corrugated channel can be perceived as a geometry in between two flat-walled channels of different width, as schematically illustrated in Fig. 6. The two limiting cases Λ → 0 and Λ → ∞ both approach the geometry of flat-walled channels with W ≈ W_min and W ≈ W_max respectively (see Fig. 6). In the simulations analyzed above Λ/W ≈ 1 as indicated by the illustration in panel (c) in Fig. 6. Reducing Λ/W leads to a denser packing of corrugations (panel (b)). The limiting case Λ/W → 0 (panel (a)) corresponds to a flat-walled channel of width W_min. On the other hand, an increase in Λ/W leads to a channel geometry shown in panel (d) with limiting case Λ/W → ∞ (panel (e)), that again corresponds to a flat-walled channel but now of width W_max. Therefore, to first order, any corrugated channel can be thought of as a construction intermediate between the limits of two flat-walled channels of width W_min and W_max. As a consequence, the critical activity required for spontaneous flow transition of an active nematic confined in a corrugated channel falls between the critical activity of two flat-walled channels of width W_min and W_max. This was also evident in Fig. 5(b) where the root-mean-squared velocity data were plotted against activity for various wavelengths.

Fig. 6 Schematic illustrating the corrugated channel of mean width W as an intermediate construction between a flat-walled channel of width W_min = W − 2A and one of larger width W_max = W + 2A. The corrugated channel approaches the flat-walled channel of width W_min and W_max when Λ/W → 0 and Λ/W → ∞ respectively. The green arrow at the top shows the increase in critical activity for a corrugated channel of wavelength Λ when compared to that of a channel of width W_max. The yellow arrow at the bottom shows the decrease in critical activity for a corrugated channel when compared to that of a flat-walled channel of width W_min.

Since critical activity is inversely proportional to the square of channel width, critical activity required for spontaneous flow transition in a flat-walled channel of width W_min is larger than that in a channel of width W_max. Therefore, the critical activity required for the spontaneous flow transition in intermediate configurations will be smaller than that of a flat-walled channel of width W_min but greater than that of a flat-walled channel of width W_max.

Yet another way to perceive the coherent flows in corrugated channels when wavelength is increased from that corresponding to W_min to W_max is the reduction in the number of geometrical constrictions or correspondingly the increase in the volume occupied by the fluid between the crests, where swirls are present. Hence the reduction in critical activity with the increase in wavelength can also be thought of as providing an effective slip to the fluid in the bulk.

3.5 Effect of swirls on the coherent flow of active nematics in corrugated channels

We have seen that coherent flows in the center of a corrugated channel are accompanied by swirls trapped in the corrugations (Fig. 3(j)). The presence of swirls in region II (Fig. 7(a)) uniquely differentiates the coherent flows in corrugated channels compared to flat-walled channels. To elucidate the role of swirls in supporting coherent flows, we interpret the swirls as a mechanism to provide a slip velocity to coherent flows in the center of the channel (region I; Fig. 7(a)).

Fig. 7 (a) Schematic illustrating the demarcation of the active nematic fluid domain into two regions and the definition of slip length. Region I when |y| < (W_min/2) where the bulk of the fluid is located and region II when (W_min/2) < |y| < (W_max/2) where the swirls are present. The fluid in region I experiences a slip velocity due to the presence of fluid in region II. (b) Slip length calculated at the boundary of regions I and II. The boundary at y = ±W_min/2 is indicated by blue dashed lines, while the boundary accounting for slip length, i.e., y(x) = ±W_min/2 ± β(x), is shown by the red continuous line connecting the symbols. The circular arrows indicate the sense of rotation of the vorticity associated with the flow in the corrugations. (c) Variation in slip length β(x) for various channel wall amplitudes with Λ = 20 and (ii) average slip length 〈β〉 for corresponding amplitudes. (d) Variation in slip length β(x) for various wavelengths with A = 2 (x-axis is normalised by Λ) and 〈β〉 as a function of wavelength.

We model the effect of swirling fluid in the corrugations (region II) on the flow in the central region (in region I) as providing a slip at the interface of two regions, y = ±W_min/2. This model is schematically illustrated in Fig. 7(a). The extent of slip is characterised by a slip velocity v_s and the corresponding slip length β. The slip length corresponds to a fictitious distance at which the fluid velocity in region II decays to zero if the velocity gradient very close to the boundary in region I is linearly extrapolated as

The above equation assumes that viscous stress dominates the interface between region I and region II and hence utilises the Newtonian constitutive relation in evaluating the slip length β. Due to the variation in the geometry along the channel length, slip length is not a constant. The function β(x) accounts for the effect of the fluid that is present in the corrugations on the coherent flow in the center of the channel.

An example of the slip length β(x) extracted from eqn (10) for a coherent flow state in a corrugated channel is shown in Fig. 7(b). The horizontal dashed blue line demarcates region I from region II. The dashed red lines connecting the symbols indicate the slip length β(x) from the boundary y = ±W_min. As expected the slip length goes to zero at the constriction, where the interface separating region I and II coincides with the rigid boundary. The slip length is a periodic function but not mirror symmetric with respect to its zero value due to the asymmetry of the swirls in region II (Fig. 3(j)). The slip length is positive on the left hand side and negative on the right hand side of the corrugation when the coherent flow is from left to right, as shown in Fig. 7(b). The positive slip length indicates a support in maintaining the coherent flow. The negative slip length indicates the resistance offered by the swirls towards the coherent flow on the right hand side of the corrugation, which occurs due to the presence of the counter clockwise rotating swirls in this case. However, as evident from Fig. 7(b), the extent of positive slip on the left side is more than that on the right side, indicating that the fluid in region I experiences a net positive slip due to the presence of fluid in region II. Thus, the overall effect of active nematic fluid entrapped in the corrugations is to endow a positive slip to the coherent flow in the center of the channel.

Since the amount of fluid contained in the corrugations varies with the amplitude and wavelength which changes the behaviour of swirls in corrugations, the slip length β(x)/A is plotted for various amplitudes and wavelengths of the corrugations in Fig. 7(c) and (d). In Fig. 7c and d(i), the abscissa is normalised with Λ for comparison between various cases. It is clear that the slip length increases with amplitude.

The variation with wavelength is more intricate (Fig. 7d(i)). To better understand how the effective slip length calculated from eqn (10) varies with wavelength, we average the slip length over an integer number of wave, . The average slip length is plotted separately as a function of amplitude and wavelength in Fig. 7c and d(ii). Increasing the amplitude increases the average slip length monotonically. However, increasing the wavelength seems to show two distinct behaviours: when the wavelength is small (Λ/W < 1), slip length increases with the increase in wavelength, but it becomes independent of wavelength for large wavelength (Λ/W > 1).

To understand the implications of these variations we analyze the flow transition in corrugated channels through an analytical approach. Since the fluid in the center of the channel (region I) exhibits a flow transition similar to that in a flat-walled channel, the role of the slip in the spontaneous flow transition may be analyzed using a linear stability analysis (see Appendix C for detailed calculations). We consider an active nematic fluid contained in a flat-walled channel of width W. The channel walls are characterised by a slip length β, resulting in a slip velocity at the channel boundary. A quasi 1D channel geometry is assumed, i.e., flow occurs only along the channel length but neither the flow nor director field exhibits any variations along the channel length. As described above we assume that viscous stresses dominate the boundary and hence eqn (10) based on the Newtonian constitutive relation is used as a boundary condition on the channel walls. Following ref. 30 and 55, the linearised governing equations for the director field θ(y) and the stream-wise velocity field v_x(y) are

ρ∂_tv_x = ∂_y(μ∂_yv_x + [small alpha, Greek, circumflex]S sin 2θ). (12)

Subjected to the slip boundary conditions and strong homeotropic anchoring , we seek stationary solutions of eqn (11) and (12). The condition for the existence of flow for a nonzero activity is found to be

1 − cos TW + βT sin TW = 0, (13)

where (see Appendix C) is an inverse length scale. Together, TW represents a dimensionless activity number that represents the competition between confinement and active nematicity. The solution of eqn (13) prescribes the critical activity [small alpha, Greek, circumflex]_c,s required for spontaneous flow transition in a flat-walled channel endowed with slip boundary conditions. An explicit expression for critical activity cannot be written down due to the non-linearity of eqn (13), and must be found numerically. However, for small slip lengths (β/W → 0), eqn (13) can be further simplified to obtain an explicit expression for critical activitywhere [small alpha, Greek, circumflex]_c,0 is the critical activity required for spontaneous flow transition in a flat-walled channel with no-slip boundary conditions (eqn (6)). Clearly, eqn (14) recovers the correct limit of [small alpha, Greek, circumflex]_c,s → [small alpha, Greek, circumflex]_c,0 when the slip length β = 0. Further, the critical activity is reduced in a channel with slip boundary conditions compared to that with no-slip boundary conditions, and the reduction is proportional to square of the slip length at leading order.

We now compare the results of the linear stability analysis with that from AN-MPCD simulations in Fig. 8. Firstly, the root-mean-square velocity V_rms is plotted against activity for a corrugated channel of W_min = 16 and A = 2 in comparison with a flat-walled channel of width W = 16 (Fig. 8(a)) and in Fig. 8(b)V_rms is plotted against activity for a corrugated channel of W_min = 12 and A = 4 in comparison with that in a flat-walled channel of width W = 12. In each figure the fluid in the corrugated channel undergoes spontaneous flow transition at a lower activity compared to that of a flat-walled channel. This reduction is indicated in the figure using yellow arrows. The reduction in the value of critical activity for the corrugated channel can be attributed to the slip velocity provided by the fluid trapped in the corrugations, consistent with the predictions of the linear stability analysis. The fluid in the center of the channel (region I) is separated from the no-slip walls because of the presence of region II. On the other hand, in a flat-walled channel, the fluid is always in contact with the channel wall and the no-slip boundary condition on the channel walls mandates a higher activity for a spontaneous flow transition.

Fig. 8 Reduction in the critical activity for spontaneous flow transition in a corrugated channel of minimum width W_min compared to a flat walled channel of width W = W_min. V_rms is plotted as a function of α for (a) a corrugated channel of W_min = 16 and A = 2 in comparison with a flat-walled channel of W = 16 and (b) a corrugated channel of W_min = 12 and A = 4 in comparison with a flat-walled channel of W = 12. (c) Results from the linear stability analysis – the critical activity in a flat-walled channel with slip boundary conditions, α_c,s, normalised with that in a channel with no-slip boundary conditions, α_c,0, is plotted against the slip length as a continuous line (full solution given by eqn (13)). The approximate solution given by eqn (14) is shown as the dashed line. Critical activity for the corrugated channel, α_c, is normalised with that of a flat-walled channel α_c,0, obtained from AN-MPCD simulations, and is plotted against the normalised average slip length, 〈β〉/W_min which are shown by the markers in the figure.

We now compare the prediction of the linear stability analysis, eqn (13), quantitatively with the results from AN-MPCD simulations. The critical activity [small alpha, Greek, circumflex]_c,s in a flat-walled channel endowed with a slip velocity is plotted as a function of slip length β (Fig. 8(c)). In this plot, the critical activity is normalised by [small alpha, Greek, circumflex]_c,0, the critical activity in a rigid channel that imposes no-slip boundary conditions (eqn (6)). The slip length is normalised with the channel width W. As expected, both the full solution (eqn (13)) and the approximate solution (eqn (14)) of the critical activity decrease with increasing slip length. The AN-MPCD data confirm that the critical activity decreases with slip length. The significant match between the results from the linear stability analysis and the AN-MPCD simulations suggests that the effect of active swirling flows that are generated in the corrugations due to the distortions in the director field is equivalent to an effective slip for the active coherent flow that develops in the center of the corrugated channel.

4 Conclusions

In this work, we have performed numerical simulations using the active nematic multi-particle collision dynamics (AN-MPCD) algorithm to study the spontaneous flow transition and the resulting flow state of an active nematic confined in a corrugated channel. Active nematics confined between two parallel plates is a well studied system.^34,44,60 At low activities, the active nematic fluid is stationary with nematic elasticity dominating the system. At sufficiently high activity, the nematic fluid undergoes a spontaneous flow transition. The resulting velocity of the fluid is proportional to the square root of activity, while the critical activity is inversely proportional to square of the channel width. In this work, we first demonstrated these scalings for an AN-MPCD fluid.

We extended channel-confined active fluids to study the behaviour of active nematics confined in corrugated channels. In this study, the wavy walls on the top and bottom of the corrugated channel are out-of-phase, but qualitatively similar conclusions are reached in a snaking channel in which the bottom and top walls are in-phase (Appendix D). We find that the distinguishing feature of the spontaneous flow transition in corrugated channels is the transition from a weak flow state to a strong flow state. This is in contrast to a flat-walled channel where the transition is from a no-flow state to a flowing state. The weak flows before the transition originate from distortions in the director field arising from the preferred anchoring of the active nematic to the curved boundaries. The critical activity required for flow depends upon the channel geometry. Increasing the amplitude of the corrugations while keeping the mean width of the channel fixed increases the constrictions in the channel and thus increases the critical activity required for flow transition. On the other hand, changing the wavelength of the corrugations can be understood by regarding the corrugated channel as intermediate between two flat-walled channels of different widths, W_min and W_max. The critical activity for the corrugated channel is smaller than that for a flat-walled channel of width W_min and higher than that for a flat-walled channel of width W_max. The current study limited itself to the transition to unidirectional flow. However, it is well-known that bidirectional flow^32,56,57 and higher modes frequently occur^13,34,55 in flat channels. Future work should consider how undulations, slip length or roughness affect the stability of bidirectional flow states.

At activities above the flow transition, the boundary-induced swirling flows are confined to the corrugations. These swirls can be contrasted with Newtonian flows in corrugated micro-channels: (i) Newtonian fluids at low Reynolds number do not generate swirls in symmetric corrugated channels due to reversibility⁶¹ and (ii) swirls generated by a Newtonian fluid at high Reynolds number are driven by the flow in the center of the channel.^62,63 In contrast the swirls observed in the active nematics in corrugations are generated by the active stress arising from the distortions in the director field. Our analysis shows that the effect of these swirls can be modelled as an effective slip to the coherent fluid flow that occurs in the bulk. Using this principle, we performed a linear stability analysis to determine the effect of the slip mechanism on the spontaneous flow transition, and matched the results from the simulations of corrugated channels with predictions from the stability analysis. Our study reveals the physical mechanisms at play when active nematics are confined in corrugated channels and further highlights the significant role of the boundaries in determining the dynamics of confined active systems.

Finally, our studies also suggest the possible mechanism by which asymmetric notches/teeth on the channel wall can direct coherent flows of channel confined active nematics, as seen in experiments.³³ Active swirls generated inside the notches due to deformations in the director field provide a slip to the coherent flow, and thus, the chirality of the stronger and bigger swirl of the asymmetric notches dictates the direction of the flow in the center of the channel. Future studies should focus on validating this hypothesis and investigating the role of asymmetry of the channel corrugations in dictating the direction of coherent flows.

Data availability

The MPCD code used in the simulations reported in the article is available at https://github.com/Shendruk-Lab/MPCD.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Nematic multi-particle collision operator

The nematic collision operator Ξ^N_i is a version of the Anderson-thermostatted collision operator^52,64where ξ_i is a random velocity drawn from the Maxwell–Boltzmann distribution for a thermal energy k_BT and 〈ξ_j〉 is the average of the random velocities generated for all the particles within each cell. Subtracting the average velocities ensures conservation of momentum when all MPCD particles have the same mass. The instantaneous moment of inertia for a cell is given as for point particles relative to the cell's center of mass r^cm, where . The third term in eqn (15) removes spurious angular momentum introduced by collision operation or by the rotation of the nematogens δ[script L]_ori. Including δ[script L]_ori in Ξ^N_i ensures conservation of angular momentum and accounts for the nematic back flow. The rotation of the nematogens arises from the stochastic orientational collision operation (eqn (4)) and from the response of nematogens to velocity gradients, which is accounted for through Jeffery's equation

δu^J_i = δtχ[u_i·Ω + λ(u_i·E − u_iu_iu_i : E)] (16)

where λ is a bare tumbling parameter and χ is the heuristic shear coupling coefficient in a flow with shear rate tensor E = [∇v + (∇v)]/2 and vorticity tensor Ω = [−∇v − (∇v)]/2. The heuristic shear coupling coefficient χ tunes the hydrodynamic susceptibility of orientation to velocity gradients. MPCD particles rotated by torques are balanced by applying equal and opposite torque to fluid by a change in angular momentum , where γ_R is the viscous rotation coefficient.

B Corrugated boundary implementation

In AN-MPCD, boundary conditions comprise two pieces: boundary surfaces, which are sinusoids, and the boundary rules which are applied to fluid particles when they violate the boundary surfaces. To implement wavy walls we follow Wamsler et al.⁵³

B.1 Boundary surfaces

The boundaries are represented by surfaces [scr S, script letter S]_b(r) = 0 expressed in the formwhere b is the index of the surface. The first term is the equation of a plane, where q_b sets the position of the bth boundary and A_b denotes the plane's normal vector. The second term incorporates the corrugations with amplitude B_b and wavelength Λ_b. Even though the simulations are in 2D, the cross product and projection on ẑ inside the cosine are written for conciseness. More complicated functions are required for non-planar wavy surfaces in 3D.⁵³

When the ith MPCD particle resides on the bth surface, [scr S, script letter S]_b(r_i) = 0. At all other points, [scr S, script letter S]_b(r_i) ≠ 0. Whenever the MPCD particle is inside the control volume, [scr S, script letter S]_b(r_i) > 0, but whenever the MPCD particle has impinged on the wall, [scr S, script letter S]_b(r_i) < 0. If an MPCD particle impinges on a wall, its ballistic trajectory is ray-traced back to the surface and a collision event is initiated. The collision is represented by a set of boundary rules.

B.2 Boundary rules

The surface boundary defined by eqn (17) exists at [scr S, script letter S]_b(r) = 0. Therefore for a particle with position r_i, if [scr S, script letter S]_b(r) ≥ 0, then it is defined as to be inside the channel. Whenever an MPCD particle is outside a boundary, a set of rules defines how the position r_i(t), velocity v_i(t) and orientation u_i(t) are transformed.

The boundary rules are applied in directions that depend on the local surface normal

at the collision point.

Upon crossing the boundary b, particle i's position is updated to r_i → r_i + [scr D, script letter D]_bn_b, where [scr D, script letter D]_b represents a shift in particle's position in the normal direction. Periodic boundary conditions (Section B.3.1) set [scr D, script letter D]_b equal to the system size, whereas impermeable walls set [scr D, script letter D]_b = 0.

The particle velocity upon colliding with the boundary b is updated using parameters [scr M, script letter M]_b,n (normal) and [scr M, script letter M]_b,t (tangential) and is transformed as v_i → [scr M, script letter M]_b,n(n_b ⊗ n_b)·v_i + [scr M, script letter M]_b,t(1 − n_b ⊗ n_b)·v_i. Here, (n_b ⊗ n_b) and (1 − n_b ⊗ n_b) are normal and tangential projection operators. Periodic boundary conditions do not change a particle's velocity and so [scr M, script letter M]_b,n = [scr M, script letter M]_b,t = +1. On the other hand, impermeable no-slip boundaries require bounce-back rules with [scr M, script letter M]_b,n = [scr M, script letter M]_b,t = −1.

The orientation of the particle also has multiplicative operators and is transformed to u_i → μ_b,n(n_b ⊗ n_b)_b·u_i + μ_b,t(1 − n_b ⊗ n_b)·u_i, which is rescaled back to the unit vector. Here, μ_b,n and μ_b,t are normal and tangential orientation parameters used to update particle orientation upon crossing the boundary. Free anchoring is implemented through μ_b,n = μ_b,t = 1 such that the orientation is unchanged. Setting μ_b,t = 0 imposes homeotropic anchoring, while μ_b,n = 0 results in planar anchoring. In this study, we consider only homeotropic channel walls (μ_b,t = 0).

B.3 Boundary conditions

The system is enclosed by four boundaries. Two of these are the left periodic boundary condition (PBC) b = 0 and the right PBC b = 1. The other two are the impermeable no-slip walls at the bottom b = 2 and top b = 3.

B.3.1 Periodic boundary condition. In this study, the PBCs are implemented along the x-direction. This PBC is achieved by updating the particle location by the length of the simulation domain in the direction normal to the boundary through which the particle exits. Velocity and orientation remain unchanged. Therefore, [scr M, script letter M]_b,n = [scr M, script letter M]_b,t = μ_b,n = μ_b,t = 1 and [scr D, script letter D]_b,n = L, where L ∈ {100, 200} is the channel length. For the PBC on the left of the system, A₀ = x̂, q₀ = 0 and B₀ = 0. For the PBC on the right, A₁ = −x̂, q₁ = Lx̂ and B₁ = 0.

B.3.2 Impermeable, no-slip wavy walls. The bottom wall is defined by the surface parameters A₂ = ŷ and q₂ = 0, with variable corrugation parameters B₂ = B and Λ₂ = Λ. Likewise, the bottom wall is defined by A₃ = −ŷ and q₂ = Wŷ (for channel width W = 20), with the same variable corrugation parameters B₃ = B and Λ₃ = Λ as the bottom wall.

On these solid walls, the no-slip boundary condition is applied by enforcing the bounce-back rule [scr M, script letter M]_b,n = [scr M, script letter M]_b,t = −1 for b ∈ {2, 3}. The particle is rewound to the location where it crossed the boundary so that it streams with update velocity for rest of the time step. Solid walls do not impose a translational shift, [scr D, script letter D]_b,n = 0. The intersection of boundaries with MPCD cells results in lower particle density as part of the cell is inaccessible to MPCD particles, which in turn lowers the local viscosity⁶⁴ and causes an effective slip. Therefore, to enforce a perfect no-slip boundary condition on solid walls, ‘phantom’ particles are added to ensure correct average density in the cell.^65,66

Here, the walls are chosen to be homeotropic, so μ_b,t = 0 for b ∈ {2, 3} However, for the no-slip boundary condition, the boundary rules are insufficient for imposing strong anchoring because the particles that cross the boundary are immediately mixed with particles that have not.⁶⁷ To ensure strong anchoring, the particles in the cells intersecting the boundary [scr S, script letter S]_b are all aligned along n_b. Since the intersecting cells are subjected to anchoring conditions, this results in the generation of strong anchoring.⁶⁷

C Linear stability analysis in a flat walled channel with slip boundary conditions

Neglecting the variation in the channel geometry and slip length along x, we calculate the critical activity for onset of spontaneous flow transition in a flat walled channel endowed with slip boundary conditions.

In addition to the equation for conservation of mass, the generalised Navier-Stokes equation for the active nematic fluid is given as

ρ(∂_t + v·∇v) = ∇·σ (19)

where ρ is the density of the fluid and σ is the total stress tensor and is given as σ = 2μE − Pδ + σ^(e) + [small alpha, Greek, circumflex]Q. In the definition of the total stress tensor (σ), μ corresponds to the fluid viscosity, E is the strain-rate tensor , P is the pressure and σ^(e) is the elastic stress tensor. The elastic stress tensor is given as σ^(e)_ij = −[small klambda, Greek, circumflex]SH_ij + Q_ikH_kj − H_ikQ_kj, where H_ij is the molecular tensor and [small klambda, Greek, circumflex] is the flow aligning parameter. The molecular tensor is given as H_ij = −δF/δQ_ij.

In a corrugated channel, coherent flow in region I is predominantly unidirectional and thus assuming a translationally invariant flow in the x direction, the flow field is described by the velocity field v_x = v_x(y) and v_y = 0. The strain-rate tensor E has only one no-zero component E_xy = ∂_yv_x, and thus eqn (19) is reduced to the following form:

ρ∂_tv_x = ∂_yσ_xy. (20)

On solving for the total stress tensor σ_xy and neglecting non-linear terms and elasticity, the aforementioned equation can be written as

ρ∂_tv_x = ∂_y(μ∂_yv_x + [small alpha, Greek, circumflex]S sin 2θ). (21)

The hydrodynamic equation for the nematic tensor order parameter Q can be written in the form

[∂_t + v·∇]Q = [small klambda, Greek, circumflex]SE + Q·Ω − Ω·Q + γ⁻¹H (22)

where Ω is the vorticity tensor which is given as . Following ref. 55 and decoupling the equation for the orientation angle θ of the nematogens from that of the scalar order parameter S, we obtain

Assuming that the viscous stress dominates at the boundary, and hence using the Newtonian constitutive relation, the boundary conditions can be written as

For the linear stability analysis, we consider φ(y,t) = φ₀ + εφ₁(y,t), where φ = {θ,v_x}, φ₀ = {π/2,0} and ε ≪ 1. Substituting θ and v_x in eqn (21) and (23) and considering the system at steady state yields the linearised system of equations:

μ∂_y²v₁ − 2[small alpha, Greek, circumflex]S∂_yθ₁ = 0 (26)

Solving the coupled homogeneous differential equations we obtain

where C₁ and C₂ are constants to be determined and

Applying the boundary conditions we obtain the condition for the critical activity for spontaneous flow transition:

This equation is further simplified by considering ε = β/W as a small parameter, and the critical activity is obtained as

D In-phase corrugated channel

Here we present some results of active nematic confined in a corrugated channel, with the top and bottom walls in-phase, unlike the out-of-phase channel walls described in the main manuscript. The boundary conditions are as mentioned in Section 2.2: homeotropic anchoring for the director field and no-slip for the velocity field on the walls. The channel walls have an amplitude A = 1, wavelength = 20 and a mean width of W = 20.

When the activity α is less than critical activity α_c, the channel with in-phase walls shows a lattice of swirling flow structures centred on the centreline (Fig. 9(a)). Since the curvature of the top and bottom channel walls is the same, the anchoring boundary conditions (namely the geometry) will result in variations in the director field primarily along the channel length. The variations are in an anti-symmetric fashion such that the corresponding flow generated will be in opposite directions. Hence, vortices span the width of the channel (Fig. 9(a) in contrast to Fig. 3(d)). It is interesting to note that the structure of the flow is similar to the well-known dancing flow state in confined active nematics, but now with a major difference: here, the lattice structure of the flow vortices arises purely from the geometry without the presence of topological defects^31,35,36,58 (Fig. 9(b)).

Fig. 9 Snapshots of the flow field (left column), director field (middle column) and active forcing (right column) in channels with in-phase corrugated walls. In the left column, the flow field is shown, which is represented by streamlines superimposed on the vorticity of the fluid. The colorbar is shown at the bottom, with blue and red representing clockwise and anticlockwise rotation respectively. In the middle column, the director field is shown by a black dashed line color shaded by the scalar order parameter. In the right column, the active force is illustrated by black arrows which are color coded by the magnitude of the force. Figures (a)–(c) correspond to low activity and figures (d)–(f) correspond to higher activity at which coherent flows prevail. All channels have a mean width of W = 20, amplitude A = 1 and wavelength Λ = 20.

Upon spontaneous flow transition, symmetry along the channel length is broken, and a coherent flow is observed. The sinusoidal variations in the active force field will follow the geometry of the channel with in-phase walls as shown in Fig. 9(f), which is unlike the reflection symmetry in the active force field observed in the case of channels with out-of-phase walls (Fig. 3(l)). A comparison of root mean square velocity of the active nematic fluid as a function of activity is shown in Fig. 10(a) for both in-phase and out-of-phase channels. In accordance with the discussion above, the strength of the flow before the flow transition is different in the two channels, but the strength of the flow is comparable in the two channels after the flow transition. Fig. 10(b) shows a comparison of the slip length as a function of channel length in both types of channels. It may be seen that the magnitude and variation of slip length β in both out-of-phase and in-phase channels are similar. It is clear from the above discussions that the qualitative features and the mechanisms at play in channels with walls of broken symmetry remain the same.

Fig. 10 (a) Comparison of spontaneous flow transition observed in channels with in-phase and out-of-phase walls by plotting root mean square velocity as a function of activity. (b) Variation in slip length β(x) for the in-phase and out-of-phase corrugated channels at α = 0.024. The parameters used in the simulations are same as those mentioned in the caption of Fig. 9.

Acknowledgements

SPT acknowledges the support by the Department of Science and Technology, India, via the research grant CRG/2023/000169. This research has received funding (T. N. S.) from the European Research Council under the European Unions Horizon 2020 research and innovation programme (Grant Agreement No. 851196).

Notes and references

1. B. Martnez-Prat, J. Ignés-Mullol, J. Casademunt and F. Sagués, Nat. Phys., 2019, 15, 362–366.
2. T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431–434.
3. P. Chandrakar, M. Varghese, S. A. Aghvami, A. Baskaran, Z. Dogic and G. Duclos, Phys. Rev. Lett., 2020, 125, 257801.
4. P. Chandrakar, M. Varghese, S. Aghvami, A. Baskaran, Z. Dogic and G. Duclos, Phys. Rev. Lett., 2020, 125, 257801.
5. H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 14308–14313.
6. S. Liu, S. Shankar, M. C. Marchetti and Y. Wu, Nature, 2021, 590, 80–84.
7. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216.
8. J.-M. Armengol-Collado, L. N. Carenza, J. Eckert, D. Krommydas and L. Giomi, Nat. Phys., 2023, 19, 1773–1779.
9. O. Hallatschek, S. S. Datta, K. Drescher, J. Dunkel, J. Elgeti, B. Waclaw and N. S. Wingreen, Nat. Rev. Phys., 2023, 5, 407–419.
10. Z. Lin, C. Gao, M. Chen, X. Lin and Q. He, Curr. Opin. Colloid Interface Sci., 2018, 35, 51–58.
11. M. N. Popescu, Langmuir, 2020, 36, 6861–6870.
12. J. Wang, M.-J. Huang, R. D. Baker-Sediako, R. Kapral and I. S. Aranson, Commun. Phys., 2022, 5, 176.
13. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246.
14. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
15. R. Alert, J. Casademunt and J.-F. Joanny, Annu. Rev. Condens. Matter Phys., 2022, 13, 143–170.
16. R. K. Gupta, R. Kant, H. Soni, A. Sood and S. Ramaswamy, Phys. Rev. E, 2022, 105, 064602.
17. P. Arora, A. Sood and R. Ganapathy, Phys. Rev. Lett., 2022, 128, 178002.
18. R. Alert, J. Casademunt and J.-F. Joanny, Annu. Rev. Condens. Matter Phys., 2022, 13, 143–170.
19. I. S. Aranson, Rep. Prog. Phys., 2022, 85, 076601.
20. L. Balasubramaniam, R.-M. Mège and B. Ladoux, Curr. Opin. Genet. Dev., 2022, 73, 101897.
21. S. P. Thampi, Curr. Opin. Colloid Interface Sci., 2022, 61, 101613.
22. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
23. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.
24. 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.
25. D. L. Koch and G. Subramanian, Annu. Rev. Fluid Mech., 2011, 43, 637–659.
26. S. P. Thampi, R. Golestanian and J. M. Yeomans, Europhys. Lett., 2014, 105, 18001.
27. S. P. Thampi and J. M. Yeomans, Eur. Phys. J.: Spec. Top., 2016, 225, 651–662.
28. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. Cristina Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365.
29. L. C. Head, C. Doré, R. R. Keogh, L. Bonn, G. Negro, D. Marenduzzo, A. Doostmohammadi, K. Thijssen, T. López-León and T. N. Shendruk, Nat. Phys., 2024, 20, 492–500.
30. R. Voituriez, J. F. Joanny and J. Prost, Europhys. Lett., 2005, 70, 404.
31. A. Samui, J. M. Yeomans and S. P. Thampi, Soft Matter, 2021, 17, 10640–10648.
32. G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost and P. Silberzan, Nat. Phys., 2018, 14, 728–732.
33. K.-T. Wu, J. B. Hishamunda, D. T. N. Chen, S. J. DeCamp, Y.-W. Chang, A. Fernández-Nieves, S. Fraden and Z. Dogic, Science, 2017, 355, eaal1979.
34. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921.
35. T. N. Shendruk, A. Doostmohammadi, K. Thijssen and J. M. Yeomans, Soft Matter, 2017, 13, 3853–3862.
36. S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, Phys. Rev. Lett., 2020, 125, 148002.
37. H. Wioland, E. Lushi and R. E. Goldstein, New J. Phys., 2016, 18, 075002.
38. M. Varghese, A. Baskaran, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2020, 125, 268003.
39. M. Ravnik and J. M. Yeomans, Phys. Rev. Lett., 2013, 110, 026001.
40. C. Joshi, Z. Zarei, M. M. Norton, S. Fraden, A. Baskaran and M. F. Hagan, Soft Matter, 2023, 19, 5630–5640.
41. S. Chen, P. Gao and T. Gao, J. Fluid Mech., 2018, 835, 393–405.
42. C. Joshi, Z. Zarei, M. M. Norton, S. Fraden, A. Baskaran and M. F. Hagan, Soft Matter, 2023, 19, 5630–5640.
43. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
44. A. Zumdieck, R. Voituriez, J. Prost and J. F. Joanny, Faraday Discuss., 2008, 139, 369–375.
45. A. J. Houston and G. P. Alexander, New J. Phys., 2023, 25, 123006.
46. S. P. Thampi, A. Doostmohammadi, T. N. Shendruk, R. Golestanian and J. M. Yeomans, Sci. Adv., 2016, 2, e1501854.
47. S. Ray, J. Zhang and Z. Dogic, Phys. Rev. Lett., 2023, 130, 238301.
48. T. Kozhukhov and T. N. Shendruk, Sci. Adv., 2022, 8, eabo5788.
49. H. Híjar and A. Majumdar, Soft Matter, 2024, 20, 3755–3770.
50. R. R. Keogh, T. Kozhukhov, K. Thijssen and T. N. Shendruk, Phys. Rev. Lett., 2024, 132, 188301.
51. T. Kozhukhov, B. Loewe and T. N. Shendruk, arXiv, 2024, preprint, arXiv:2401.17777 DOI:10.48550/arXiv.2401.17777.
52. T. N. Shendruk and J. M. Yeomans, Soft Matter, 2015, 11, 5101–5110.
53. K. Wamsler, L. C. Head and T. N. Shendruk, Soft Matter, 2024, 20, 3954–3970.
54. G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, in Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids, ed. C. Holm and K. Kremer, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, pp. 1–87.
55. L. Giomi, L. Mahadevan, B. Chakraborty and M. F. Hagan, Nonlinearity, 2012, 25, 2245.
56. Z.-Y. Li, D.-Q. Zhang and B. Li, Phys. Rev. Res., 2021, 3, 023253.
57. V. J. Pratley, E. Caf, M. Ravnik and G. P. Alexander, Commun. Phys., 2024, 7, 127.
58. S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, Soft Matter, 2019, 15, 1597–1604.
59. R. R. Keogh, S. Chandragiri, B. Loewe, T. Ala-Nissila, S. P. Thampi and T. N. Shendruk, Phys. Rev. E, 2022, 106, L012602.
60. 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.
61. L. G. Leal, Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes, Cambridge University Press, 2007.
62. T. Nishimura, Y. Ohori and Y. Kawamura, J. Chem. Eng. Jpn., 1984, 17, 466–471.
63. F. Oviedo-Tolentino, R. Romero-Méndez, A. Hernández-Guerrero and B. Girón-Palomares, Int. J. Heat Fluid Flow, 2008, 29, 1233–1239.
64. H. Noguchi, N. Kikuchi and G. Gompper, Europhys. Lett., 2007, 78, 10005.
65. A. Lamura, G. Gompper, T. Ihle and D. M. Kroll, Europhys. Lett., 2001, 56, 319.
66. D. S. Bolintineanu, J. B. Lechman, S. J. Plimpton and G. S. Grest, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 066703.
67. L. C. Head, Y. A. G. Fosado, D. Marenduzzo and T. N. Shendruk, Soft Matter, 2024 10.1039/D4SM00436A.

---