Active nematic-isotropic interfaces in channels

TL;DR

This study uses numerical simulations to explore the behavior of active nematic-isotropic interfaces in channels, revealing how activity influences interface stability, velocity, and defect dynamics in confined active liquid crystals.

Contribution

It introduces an advanced hybrid simulation model to analyze the hydrodynamics and stability of active nematic-isotropic interfaces under various activity levels and temperatures.

Findings

Interface velocity increases with activity, linearly for contractile and quadratically for extensile systems.

Active forces cause nematic expansion in contractile and contraction in extensile systems.

Defect nucleation and dynamics occur at high activities, with regular behavior observed.

Abstract

We use numerical simulations to investigate the hydrodynamic behavior of the interface between nematic (N) and isotropic (I) phases of a confined active liquid crystal. At low activities, a stable interface with constant shape and velocity is observed separating the two phases. For nematics in homeotropic channels, the velocity of the interface at the NI transition increases from zero (i) linearly with the activity for contractile systems and (ii) quadratically for extensile ones. Interestingly, the nematic phase expands for contractile systems while it contracts for extensile ones, as a result of the active forces at the interface. Since both activity and temperature affect the stability of the nematic, for active nematics in the stable regime the temperature can be tuned to observe static interfaces, providing an operational definition for the coexistence of active nematic and isotropic phases. At higher activities, beyond the stable regime, an interfacial instability is observed for extensile nematics. In this regime defects are nucleated at the interface and move away from it. The dynamics of these defects is regular and persists asymptotically for a finite range of activities. We used an improved hybrid model of finite differences and lattice Boltzmann method with multi-relaxation-time collision operator, the accuracy of which allowed us to characterize the dynamics of the distinct interfacial regimes.

Equations

$$f_{bulk}=$$

$$+\frac{A_0\gamma }{4}(Q_{\alpha \beta }{Q}_{\alpha \beta }{)}^2,$$

$$f_{el}=$$

$${\partial }_t{Q}_{\alpha \beta }+u_{\gamma }{\partial }_{\gamma }{Q}_{\alpha \beta }-S_{\alpha \beta }(W_{\alpha \beta },Q_{\alpha \beta })=\Gamma H_{\alpha \beta },$$

$$S_{\alpha \beta }=$$

$$-2\xi \left(Q_{\alpha \beta }+\frac{{\delta }_{\alpha \beta }}{3}\right)(Q_{\gamma ϵ}{\partial }_{\gamma }{u}_{ϵ}),$$

$$H_{\alpha \beta }=-\frac{\delta \mathcal{F}}{\delta Q_{\alpha \beta }}+\frac{{\delta }_{\alpha \beta }}{3}\mathrm{Tr}\left(\frac{\delta \mathcal{F}}{\delta Q_{\gamma ϵ}}\right).$$

$${\partial }_t\rho +{\partial }_{\alpha }(\rho u_{\alpha })=0,$$

$$\rho {\partial }_t{u}_{\alpha }+\rho u_{\beta }{\partial }_{\beta }{u}_{\alpha }={\partial }_{\beta }{\Pi }_{\alpha \beta }+\eta {\partial }_{\beta }\left({\partial }_{\alpha }{u}_{\beta }+{\partial }_{\beta }{u}_{\alpha }\right).$$

$${\Pi }_{\alpha \beta }^{\text{active}}=-\zeta Q_{\alpha \beta },$$

$${\Pi }_{\alpha \beta }^{\text{passive}}=$$

$$-\xi H_{\alpha \gamma }\left(Q_{\gamma \beta }+\frac{{\delta }_{\gamma \beta }}{3}\right)-\xi \left(Q_{\alpha \gamma }+\frac{{\delta }_{\alpha \gamma }}{3}\right)H_{\gamma \beta }$$

$$-{\partial }_{\alpha }{Q}_{\gamma \nu }\frac{\delta \mathcal{F}}{\delta ({\partial }_{\beta }{Q}_{\gamma \nu })}+Q_{\alpha \gamma }{H}_{\gamma \beta }-H_{\alpha \gamma }{Q}_{\gamma \beta }.$$

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)-f_i(\mathbf{x},t)$$

$$={\mathbf{M}}^{-1}\mathbf{RM}[f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t)]\Delta t+{\mathcal{S}}_i,$$

$$f_i^{eq}=\rho w_i\left[1+\frac{{\mathbf{c}}_i\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{c}}_i\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{{\mathbf{u}}^2}{2c_s^2}\right],$$

$$\rho =\sum _i{f}_i,\rho \mathbf{u}=\sum _i{\mathbf{c}}_i{f}_i+\frac{{\mathbf{F}}_i\Delta t}{2}.$$

$$S=\frac{S_n}{2}\left[\tanh \left(\frac{x-x_1}{2\lambda }\right)-\tanh \left(\frac{x-x_2}{2\lambda }\right)\right],$$

$$S_n=\frac{\gamma +\sqrt{3(3{\gamma }^2-8\gamma )}}{4\gamma }.$$

$$S=\frac{S_n}{2}\left[1-\tanh \left(\frac{x-x_i}{2\lambda }\right)\right].$$

$${\mathbf{F}}_{\perp }^{\text{active}}$$

$$\displaystyle+S\,\mathbf{m}\cdot(\mathbf{n}\cdot\boldsymbol{\nabla})\mathbf{n}-\frac{1}{3}\textbf{m}\cdot\boldsymbol{\nabla}S\Big{]}\textbf{m}.$$

$${\mathbf{F}}_{\perp }^{\text{active}}=-\zeta \left(\frac{|\nabla S|}{3}+\frac{S}{R}\right)\mathbf{m}.$$

$$\kappa =\frac{|x_i^{\prime \prime }(y)|}{{\left[1+(x_i^{\prime }(y){)}^2\right]}^{\frac{3}{2}}}.$$

$${\mathbf{F}}_{\perp }^{\text{active}}=\zeta \left(\frac{2|\nabla S|}{3}+\frac{S}{R}\right)\mathbf{m}.$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa, P-1749-016 Lisboa, Portugal; E-mail: rcvcoelho@fc.ul.pt*††footnotetext: b* Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, P-1749-016 Lisboa, Portugal. *

1 Introduction

Systems formed by units or particles capable of transforming the energy of the environment into directed motion are known as active matter. These intrinsically non-equilibrium systems are characterized by complex dynamical behavior. Examples that occur naturally are dense suspensions of bacteria, mixtures of microtubule-kinesin, and shoals of fish (see 1, 2 and references therein). While a detailed description of a particular system is daunting, particle based simulations of simple models exploded in the last decade and reported a range of new phenomena. Motility induced phase separation (MIPS) is one of the most striking 3.

A different approach focused on the phenomenology observed in dense systems, which results from the interplay of self-propulsion and alignment interactions due to the shape of the particles, and their hydrodynamic coupling to the suspending medium 4, 1, 2, 5. Here, we follow this route and consider momentum conserving (or wet 6) active systems with nematic order, which may be described in the hydrodynamic limit by the continuum equations of liquid crystals (LC) with an extra stress term that accounts for the activity. While passive LCs 7, 8 have been studied for decades, active nematics are currently a hot topic of fundamental research. Indeed, only recently, a second active stress term was proposed, which may have a profound effect on the behaviour and stability of active nematics 9.

One question that has been investigated concerns the conditions to observe directed coherent flow, which is key to applications in micromotors and micromachines 10, 11. While in the bulk active turbulence may be unavoidable for any degree of activity 4, 1, 2, 5, confinement in narrow channels results in a spontaneous symmetry breaking and opens a window to observe directed flow at low but non-zero activities 12, 13, 14. As the activity increases, the directed flow changes to oscillatory, dancing disclinations, and ultimately turbulent regimes, both in the deep nematic phase 15, 16, 17 and also at temperatures above the passive nematic-isotropic transition 18.

The fact that the activity, by itself, induces local nematic ordering and even active turbulence, raises the question of the coexistence between active nematic and isotropic phases and of the behaviour of the interfaces between them. Although this is a topic of intense research for MIPS in active systems with scalar order parameters 3, 19, 20, 21, it has not been addressed for active nematics, despite the fact that an intrinsic free energy has been proposed for these systems 22.

In addition to fundamental issues related to the existence and characterizaton of phase coexistence in active nematics, interfaces of LCs have played a prominent role in applications. For example, they can be controlled easily in order to modify their optical properties 23 and provide an efficient means to transport colloids 24. Interfaces of active LCs are expected to play similar roles. Previous studies of active and passive LCs have focused on the dynamical behaviour of droplets of one in the other, with the observation that passive droplets may be driven through the active fluid while active droplets may rotate in passive nematics  25, 26, 27, 28. The coexistence of active nematic and isotropic phases and the dynamics of their interfaces have not been investigated.

In this paper, we consider the dynamics of the interface between the nematic and the isotropic phases of an active LC confined in channels with strong homeotropic anchoring at the walls. The system exhibits a nematic ordering transition, driven by the temperature at fixed (low) activity. The velocity of the interface at the passive NI temperature increases from zero (i) linearly with the activity for contractile systems and (ii) quadratically for extensile ones. In the former, the nematic order increases while it decreases in the latter. This behavior is unexpected as extensile systems favor nematic order, but it is readily understood as a result of the active forces at the interface. The effect is reversed in planar channels as the active forces at the interface are also reversed, as we show here. At fixed activity, the velocity of the interface depends linearly on the temperature. In the stable regime the temperature can be tuned to observe a static interface in the channel providing an operational definition for the coexistence of confined active nematic and isotropic phases. Beyond the stable regime, extensile nematics exhibit an interfacial instability in hometropic and planar channels.

We performed computer simulations of the hydrodynamic equations of active liquid crystals using a hybrid method of lattice Boltzmann 29, 30 (LB) and a finite differences (FD) predictor corrector 31. There are challenges associated with the simulation of interfaces in multiphase LB models, such as the control of spurious velocities and the difficulty to set the physical values of both viscosity and surface tension 32. To overcome these problems, we introduced various improvements, as discussed in the following.

Our paper is organized as follows. In Sec. 2, we describe the equations of motion and the numerical method. In Sec. 3, we describe the results of numerical simulations of passive and active nematic interfaces in channels and show how the active interfaces may become static through a shift in the temperature. In Sec. 4, we discuss the interfacial instability observed in extensile nematics and mention briely the flow states that occur at higher activities. In Sec. 5, we investigate the effects of anchoring at the channel walls and check the effect of the interaction between interfaces. In Sec. 6 we summarize and conclude*†*††footnotetext: † Electronic Supplementary Information (ESI) available: details of the MRT model and videos of the interfacial dancing state and interfacial instability are given in the ESI. See DOI: 00.0000/00000000..

2 Methods

2.1 Equations of motion

We employ the Landau-de-Gennes free energy, $\mathcal{F}=\int_{V}f\,d^{3}r$, to describe the passive liquid crystal at equilibrium. The free energy density, $f$, is the sum of two terms: $f=f_{bulk}+f_{el}$, where $f_{bulk}$ is the bulk free energy and $f_{el}$ is the elastic energy, i.e., the energy cost associated to distortions with respect to the uniform alignment of the director field. These energy densities are given by:

[TABLE]

The tensor order parameter is assumed to be uniaxial, both in the initial conditions and in the interpretation of the results: $Q_{\alpha\beta}=S(n_{\alpha}n_{\beta}-\delta_{\alpha\beta}/3)$, where $S$ is the scalar order parameter and $n_{\alpha}$ is the $\alpha$-component of the director. The parameter $\gamma$ controls the magnitude of the order and depends on the temperature for thermotropic or on the density for lyotropic LCs, $K$ is the elastic constant and $A_{0}$ is a constant. In this one-constant approximation the anchoring at the NI interface is planar 7. Here the Greek indices stand for Cartesian coordinates with sums over the terms with repeated indices. The coexistence between the bulk nematic and isotropic phases occurs at $\gamma_{NI}=2.7$ and $S=1/3$ as obtained by minimizing $f_{bulk}$. For $\gamma>\gamma_{NI}$ the bulk equilibrium phase is nematic while for $\gamma<\gamma_{NI}$ the bulk equilibrium phase is isotropic.

The tensor order parameter evolves according to the Beris-Edwards equation:

[TABLE]

where $\Gamma$ is the collective rotational diffusive constant. The co-rotational term, $S_{\alpha\beta}$, describes the response of the LC to gradients in the velocity field $\mathbf{u}$ and it is given by:

[TABLE]

where $W_{\alpha\beta}=(\partial_{\beta}u_{\alpha}-\partial_{\alpha}u_{\beta})/2$, $D_{\alpha\beta}=(\partial_{\beta}u_{\alpha}+\partial_{\alpha}u_{\beta})/2$. The aligning parameter $\xi$ depends on the particles shape: it is positive for rod-like particles and negative for disk-like particles. The molecular field stands for:

[TABLE]

The Navier-Stokes and continuity equations give the time evolution of the velocity field in terms of the stress tensor $\Pi_{\alpha\beta}$, which depends on the tensor order parameter:

[TABLE]

[TABLE]

The stress tensor is the sum of an active and a passive terms: $\Pi_{\alpha\beta}=\Pi^{\text{passive}}_{\alpha\beta}+\Pi^{\text{active}}_{\alpha\beta}$. In the simplest description of the hydrodynamics of active LCs, which accounts for the lowest order contribution with the appropriate symmetry, the active term is given by 33:

[TABLE]

while the passive one is:

[TABLE]

In these equations $\eta$ is the shear viscosity, $P_{0}$ is the hydrostatic pressure, and $\zeta$ is the activity parameter. Figs.1 A and B are a schematic representation of the velocity field for positive and negative $\zeta$. For $\zeta>0$, the velocity field resembles the one for pushers and the system is extensile, while for $\zeta<0$, it resembles the flow for pullers and the system is contractile.

2.2 Hybrid method

To perform the numerical simulations, we solved the hydrodynamic equations of liquid crystals using a hybrid method of LB and FD. This approach is common in the literature and has successfully reproduced analytical and experimental results 4, 34, 13. The LB has many advantages if compared to other computational fluid dynamics methods, namely the way complex boundary conditions are treated, the performance in parallel architectures and the possibility to include many physical models 30, 29, 35 which justifies its use to simulate liquid crystals. The first approaches were fully LB, meaning that both the Navier-Stokes and Beris-Edwards equations were solved using LB 36***In fact, the LB solves the Boltzmann equation and recovers the macroscopic equations of motion, e.g. the Navier Stokes equation, in the macroscopic limit as can be shown through the Chapman-Enskog expansion 30.. Later, the Beris-Edwards equations started to be solved with FD since this is more efficient in memory usage and furthermore eliminates spurious terms arising from the coupling between these two equations of motion in the Chapman-Enskog expansion 13. In the LB part (to solve the Navier-Stokes equation) of the models used in the literature, the stress tensor was implemented in the equilibrium distribution and the D3Q15 lattice was used to discretize the velocity space 36, 37, 13. There are additional challenges associated with simulation of interfaces in multiphase LB models as the control of spurious velocities and the difficulty to set the physical values of both viscosity and surface tension 32. To overcome these problems, we introduced three main improvements. First, we adopted the D3Q19 lattice since it is more isotropic than the D3Q15 and can resolve with finer resolution the angular distribution of the velocity field. Second, we use the multi-relaxation-time (MRT) collision operator, which has superior stability and accuracy and allows us to independently choose the relaxation rates of the hydrodynamic moments38. The MRT is known to reduce, by orders of magnitude, the spurious velocities in the pseudopotential LB models, which are a class of multiphase and multicomponent models that suffer badly from this problem 39, 40, 41. Third, in our model, the stress tensor is implemented in the force term, which reduces the spurious velocities in the free energy models 42.

The discrete form of the Boltzmann equation with the MRT collision operator reads as follows:

[TABLE]

where the index $i=\{0,\ldots,18\}$ runs over the velocity vectors. The transformation matrix $\mathbf{M}$, and the relaxation matrix $\mathbf{R}$ are given in the Electronic Supplementary Information. The source term, $\mathcal{S}_{i}$, is given by Guo’s forcing scheme, the equilibrium moments of which are also given in the Electronic Supplementary Information (see also Ref. 30 for a detailed description of the method). The equilibrium distribution is the Maxwell-Boltzmann distribution expanded up to second order in Hermite polynomials 43:

[TABLE]

where $\mathbf{c}_{i}$ and $w_{i}$ are the velocity vectors and the discrete weights of the lattice, and $c_{s}=1/\sqrt{3}$ is the speed of sound in the D3Q19 lattice. The density $\rho$ and the macroscopic velocity $\mathbf{u}$ are calculated from the distribution functions:

[TABLE]

Notice that the velocity passed to the FD method is the one calculated with Eq. (11), which is corrected by the term $\mathbf{F}_{i}\Delta t/2$ to ensure second order accuracy. The force is calculated as $F_{\alpha}=\partial_{\beta}(\Pi_{\alpha\beta}+\rho{c_{s}}^{2}\delta_{\alpha\beta})$, where $\Pi_{\alpha\beta}$ is given by Eqs. (7) and  (8). In the FD, Eq. (2) is discretized on the same grid as the LB and solved through the predictor-corrector algorithm using second order differences. Therefore, both methods (FD and LB) are second order accurate and should provide consistent solutions of the equations of motion.

3 Stable interface

In this section, we simulate active stable interfaces (low activities) in a channel and compare the results with those of passive interfaces. We show that the active interfaces move with a velocity that depends on the activity and that they may be stabilized by a shift in the temperature.

3.1 Channel setup and initial conditions

We start by simulating an open channel, with the nematic, between $x_{1}=3L_{X}/8$ and $x_{2}=5L_{X}/8$, and the isotropic fluid elsewhere. The fluid is confined between two flat walls, at a fixed distance $L_{Y}$, with infinite homeotropic anchoring (see Fig. 2) and periodic boundary conditions in the $x$ direction. The initial velocity is set to zero and the density to $\rho=1$ everywhere. The following parameters are used: $\tau=1.5$, $L_{x}\times L_{Y}\times L_{Z}=270\times 45\times 1$, $K=0.04$, $A_{0}=0.1$ and $\Gamma=0.34$. The other parameters change for each simulation and will be given in the caption of the corresponding figures. Our results are given in lattice units: the distance between nodes is $\Delta x=1$ and the time step is $\Delta t=1$ (see Refs. 44, 4 to transform these to physical units).

Since the interfacial profile follows approximately a hyperbolic tangent we initialize the scalar order parameter as:

[TABLE]

where $\lambda=\sqrt{27K/A_{0}\gamma}$ is the correlation length (equal to 2 at the bulk nematic-isotropic coexistence $\gamma=2.7$) and the scalar order parameter of the nematic phase, $S_{n}$, is computed through the minimization of the free energy:

[TABLE]

$S_{n}=1/3$ at bulk coexistence $\gamma=2.7$. This expression for $S_{n}$ is used in the initial conditions and its value is fixed at the walls. We find the position of the interface, $x_{i}$, by fitting

[TABLE]

to the interface on the right (the interface on the left is a mirror image in the stable regime). Other fields proportional to $S$, such as $Q_{yy}$, can also be used to obtain the position of the interface. In the passive case, the interface relaxes and becomes concave 24, indicating that the nematic wets the channel walls. Indeed, a thin film of nematic is observed between the channel walls and the isotropic phase, for all the parameters considered here. In the simulations of active nematics, we turn the activity after 20000 iterations, allowing for the relaxation of the passive interface.

3.2 Passive interface in a channel

In a liquid confined between two walls, the coexistence liquid-vapour temperature changes due to the interaction of the fluid with the walls, an effect known as capillary condensation 45. This also happens in nematics and has been called capillary nematisation. At the bulk NI transition temperature ($\gamma_{NI}$) the interface between the confined nematic and isotropic phases will move, leading to the expansion of the nematic if the walls favor the ordered phase 46. We calculate numerically the NI coexistence temperature under confinement $\gamma^{ch}_{NI}$, as the temperature where the passive interface is static. We start by evaluating the position of the interface at different times through a fit of $Q_{yy}$ using Eq. (14). Fig.3 shows the interfacial velocity at five different temperatures. Since the dependence is linear it is straightforward to calculate the NI coexistence temperature in the channel $\gamma^{ch}_{NI}\approx 2.68964$ as the temperature where the interfacial velocity vanishes. At temperatures $\gamma>\gamma^{ch}_{NI}$ the velocity is positive for the interface on the right (nematic expansion) and for $\gamma<\gamma^{ch}_{NI}$ the velocity is negative (nematic contraction). The interface is stable for a range of temperatures $|\Delta\gamma|<0.02$ around $\gamma^{ch}_{NI}$. At temperatures higher or lower than these, one of the two phases becomes mechanically unstable, the interface disappears and the system becomes nematic ($\gamma>\gamma^{ch}_{NI}$) or isotropic ($\gamma<\gamma^{ch}_{NI}$). For the narrow channels considered here the strong anchoring at the walls dominates the interfacial anchoring and the anchoring at the interface is dictated by the wall anchoring (see Fig.1 C).

3.3 Active interface and stabilization by temperature

In order to proceed we consider the interface of an active nematic at the coexistence temperature of the passive system, $\gamma=\gamma^{ch}_{NI}$.

The orientation of the directors close to the interface is determined by a competition between the anchoring at the walls 47 and the interfacial passive and active anchorings 27. The latter tends to align the directors parallel to the interface in extensile systems and normal to it in contractile ones. In narrow channels with strong anchoring, as the ones considered here, the effect of the wall anchoring is dominant and the director field is uniform (Fig. 1 C). In wider channels and/or for weak anchoring at the walls the directors are expected to align parallel to the interface, both in passive LCs with a single elastic constant $K$ as well as in extensile nematics with the same elasticity (Fig. 1 D).

At low activities a stable interface is formed (Fig. 4 A and B). The interface moves with constant velocity and constant shape, while the bulk remains at rest. The domain of the nematic phase expands for contractile systems ($\zeta<0$) and contracts for extensile ones ($\zeta>0$). This behavior is counterintuitive since bulk extensile systems favor nematic ordering while contractile ones disfavor it.

To understand this result, we calculate the active forces at the interface. Projecting the active force $\mathbf{F}^{\text{active}}=-\zeta\nabla\cdot\mathbf{Q}$ on the outward normal to the NI interface, $\mathbf{m}$, we have:

[TABLE]

An estimate of the force is easily obtained by assuming that the interface is circular, with radius $R$, and that the director is parallel to it, $\mathbf{n}\cdot\mathbf{m}=0$ (Fig. 1 E), which is reasonable in the center of the channel. The normal force is then given by the sum of two terms:

[TABLE]

The first of these terms is the active force due to the gradient of the nematic order, while the second is due to the interfacial curvature. For extensile nematics, the normal active force acts inwards, while for contractile ones the force acts outwards, revealing the mechanism for the contraction and the expansion of the nematic phase observed in the simulations. In addition, the force is linear in $\zeta$ and it increases as $R$ decreases. We note that the second term of the active force, which increases the curvature of a curved interface with planar anchoring, acts in a similar way to the active forces that increase bend fluctuations in extensile nematics and are ultimately responsible for the bend instability of unconfined extensile nematics. By contrast, the first term of the force, proportional to $|\boldsymbol{\nabla}S|$, is absent deep in the nematic, where the scalar order parameter is nearly uniform. At the interface, however, this term becomes increasingly important as the curvature decreases, and it is the only active force at a flat interface. Similar active forces were reported in Ref. 25 for interfaces in active nematic emulsions, characterized by an additional conserved order parameter. One important difference between our results and those of Ref. 25 is that the interfaces considered here are stable due to the confinement in the channel.

In what follows we show that the active interfacial motion can be stopped by a shift in temperature, in a way reminiscent of the behaviour of passive nematic-isotropic interfaces near coexistence.

We start by noting that, for extensile nematics $\zeta>0$ the interfaces move towards each other and, when they are sufficiently close, the velocity field of one pushes the other rendering the interfaces static without the need for a temperature shift (see Fig. 5). This effect does not occur for contractile systems $\zeta<0$ as the two interfaces move in opposite directions. We are interested in the dynamics of a single interface, far from surfaces or other interfaces. Thus, in Fig. 6 we plot the interfacial velocity measured when the interfaces are far apart and the interaction between them is negligible.

For each system (or activity), we changed the temperature around $\gamma^{ch}_{NI}$ and measured the interfacial velocity. We found that the interfacial velocity is linear in $\Delta\gamma$, for both contractile and extensile nematics. It is then straightforward to calculate the shift $\Delta\gamma$ required to observe a static interface for a given activity (Fig. 7A). In Fig. 7B we plot the phase diagram for active nematics in a channel, as the temperature shift, $\Delta\gamma$, that renders the interface static. For contractile nematics the coexistence line is approximately linear. For extensile ones, the linear relation is observed in a narrow range of activities becoming quadratic afterwards. This quadratic behaviour correlates with the quadratic dependence of the interfacial velocity on the activity observed in extensile nematics (Fig. 6).

In order to check that the combined effects of temperature and activity lead to a static interface, we performed longer simulations (up to $t=10^{7}$ time steps) at the values of $\Delta\gamma$ given in Fig. 7B for different activities. We found that the interfacial velocities are indeed very close to zero at these temperatures (which could be refined iteratively).

The shape of the static interface changes with the activity and is plotted in the inset of Fig. 8. Although there is no net flow under these conditions, the vortices near the interface (see Fig. 4) become stronger as the activity increases. We quantify the dependence of the static interfacial shape on the activity through the curvature at the center of the channel ($y=L_{Y}/2$). The interfacial curvature is defined as:

[TABLE]

and is plotted in Fig. 8, where $x_{i}(y)$ is the position of the interface. For extensile nematics the curvature increases with the activity. Assuming that the interface is circular, we can estimate its radius and find: $R(\zeta=-0.0025)\approx 33$ and $R(\zeta=0.0009)\approx 12$.

We can now understand why the interfacial velocity has a linear regime and a non linear one, Fig. 6. When the curvature is small (large $R$), the second term in Eq. (16) is small and the force is approximately linear in $\zeta$. When the curvature increases (small $R$), the second term in Eq. (16), which depends on the activity (Fig. 8), is relevant and the active force becomes non-linear.

4 Flow states and interfacial instability

In Fig. 4 we give an overview of the different states observed as the activity increases, from contractile to extensile nematics. The top pannels (A and B) correspond to contractile (A) and extensile (B) nematics where a stable interface, with fixed shape, was found to move at constant velocity in the channel. These states are observed at low activities. At higher activities, the interface of extensile nematics exhibits an instability. A complex interface with a persistent saptio-temporal structure (panel C) is observed for a range of activities beyond the instability. At even higher activities the interface cannot be sustained, as the isotropic phase in the channel becomes unstable, and complex flow states, dancing D and active turbulence E, are observed in line with earlier reports  4, 48.

Stable interfaces have a shape that is invariant in time and mirror symmetric with respect to the centre of the channel. An interfacial instability occurs when the system breaks one or both of these symmetries (Fig. 9) at a specific value of the activity, $\zeta_{inst}$. For the initial conditions described above, the interfaces move towards each other, interact through the velocity field and become static as illustrated in Fig. 5. Fig. 9B shows that at $\zeta_{inst}\approx 0.0012$ the interfaces break the mirror symmetry and one of the two vortices close to the interface becomes stronger. At $\zeta=0.0013$, the system becomes isotropic and the interface disappears. This suggests that the system is multistable in this range of parameters and that the final state may be affected by finite size effects. At $\zeta=0.0014$ (see Fig. 9 C) a strongly distorted interface reappears. The instability is driven by the active force that increases with the interfacial curvature which in turn increases the curvature of the interfaces with planar anchoring. This instability is closely related to the bend instability for extensile systems reported previously in unconfined active nematic emulsions 25, 27. With one important difference: the active force that drives the interfacial instability is also responsible for the re-appearence of the nematic phase and the interfacial state persists at long times.

A distinct instability in active nematics, where the fluid starts moving coherently across the channel, was reported deep in the ordered phase, at higher activities that scale with the inverse of the square of the channel height 12. This scaling was verified numerically in the non-linear regime in Ref. 13 and is also approximately verified for the interfacial instability. In Fig. 10 we plot the activity at threshold $\zeta_{inst}$ for channels with different heights. For each channel $L_{Y}$, we calculated $\gamma^{ch}_{NI}$ as the passive NI temperature depends on the channel height. We found that the value of the threshold activity decreases approximately with ${L_{Y}}^{2}$ in line with the scaling reported previously. This is to be expected as the active length associated with the active vortices scales with the ratio of the elastic to the active forces, $\sqrt{K/\zeta}$ and the activity number, which is the ratio of the channel height to the active length, is nearly constant, $A=L_{Y}\sqrt{\zeta/K}\approx 7.8$, at the interfacial instability threshold.

At higher activities, in the slab geometry, an interfacial dancing state (Fig. 4C) is observed. This state is characterized by two NI interfaces and a pair of $+1/2$ defects moving in complex oscillatory trajectories. The defects leave one interface into the nematic and disappear from it at the other interface. This state resembles the bulk dancing state, confined by the NI interfaces. The interfaces are no longer concave with constant curvature, as the stable interfaces discussed in Sec. 3, but they are not transient. The interfacial dancing state is stable over a finite range of activities: for a flow aligning nematic with $\xi=0.7$ and a channel with $L_{Y}=45$, the state is observed for activities in the range $0.0015<\zeta<0.0021$.

At even higher activities, the interface disappears (the isotropic phase becomes unstable) and the previously reported dancing (Fig. 4D) and active turbulent states (Fig. 4E) are observed 4, 48.

We used the spatial average of the vorticity norm in the steady state, calculated at $t=2\times 10^{6}$, as an order parameter to characterize the states of extensile nematics: $\langle\omega\rangle$, where $\boldsymbol{\omega}=\nabla\times\mathbf{u}$. In Fig. 11A, we plot the average vorticity as a function of the activity for nematics with $\xi=0.7$. The results reveal that the average vorticity $\langle|\mathbf{\omega}|\rangle$ increases as the activity increases. More interestingly the vorticity exhibits jumps between dynamical states. The first jump, at $\zeta\approx 0.0012$, occurs between the stable interface and the interfacial dancing states. As discussed previously, we have observed multi-stability and a sequence of intermediate states (see Fig. 9) including isotropic states with zero vorticity, in the transition region. States with zero vorticity are represented by the black squares in the right panel of the figure, Fig. 11B. The second jump occurs between the interfacial dancing and the dancing states where the nematic fills the channel and the interface is no longer stable. In Fig. 11B we plot the vorticity at $t=2\times 10^{6}$ for active nematics with different alignment parameters and activities. At the bulk NI temperature the flow is tumbling if the alignment parameter $\xi<0.43$ 13 but the stability of the flow states does not depend strongly on it. However, most of the black squares (isotropic) are in the flow tumbling regime.

5 Other channels: planar anchoring and closed channel

Finally, we consider other channels to investigate the effects of the anchoring at the walls and the open boundary condition.

We start by considering planar channels with strong anchoring and scalar order parameter given by Eq. 13. For narrow channels the confined nematic is uniform and the alignment by the walls dominates the interfacial and the active anchorings. We initialize the director field in this configuration. We consider a slab geometry with the nematic in the center of the channel surrounded by isotropic. All other parameters are as in Sec. 3.1. Fig. 12, illustrates a sequence of states that resemble those found for homeotropic channels.

At low activities, the interface is stable and propagates with constant velocity depending on the activity (see Fig. 12 A and B). However, the nematic phase contracts for contractile systems and expands for extensile ones. As before we use Eq. (15), to calculate the active force for a circular interface with a normal director field, $\mathbf{n}\cdot\mathbf{m}=1$ (Fig. 1 F). The active force is now:

[TABLE]

Clearly, the active forces at interfaces with planar and homeotropic anchoring point in opposite directions (see Eq. (16)). In planar channels, the active force at the interface is also larger, in line with the larger velocities observed in the simulations.

In Fig. 12 C, we illustrate the interfacial dancing state observed in the planar channel. We note the formation of islands of nematic close to the boundaries and the presence of $+1/2$ defects that move in directions opposite to those observed in the homeotropic channel (as a result of the opposite signs of the vorticities). The dancing state, shown in Fig. 12 D, is similar in both channels.

Finally, we consider a closed channel with only one interface (see Fig. 13). The left boundary is nematic with planar anchoring and the right one is isotropic, with no slip boundary conditions at the four walls. The system is initialized with the left half of the channel in the uniform nematic state with directors parallel to the interface and the right half isotropic. We note that the stable interface behaves in the same way as described in Sec. 3. Furthermore, the interfacial dancing state appears even when the channel is closed (Fig. 13 C) but manifests itself through a different route (Fig. 13 D to G): The interface does not break the mirror symmetry but becomes increasingly curved until a defect is formed and is ejected into the nematic phase. When the interface approaches the left rigid boundary it finally breaks the mirror symmetry and becomes asymmetric as in Fig. 9.

6 Summary and conclusion

To summarize, we studied the behavior of nematic-isotropic interfaces in active liquid crystals confined in a narrow channel. We found that, at high activities, the interface disappears and the system becomes nematic in the previously reported dancing or active turbulent flow states. At lower activities, there are two flow states where an interface is present: the stable interface and the interfacial dancing state. In the stable interface regime with homeotropic anchoring, the right interface propagates forwards for contractile nematics and backwards for extensile ones, as a result of the active force at the interface.

The moving active interfaces may be stopped by a shift in the temperature. We calculated this shift for a wide range of activities and checked the robusteness of this stabilization mechanism by performing longer simulations. The static active interfaces are flatter than the passive ones for contractile nematics and more strongly curved for extensile ones. The curvature increases with the activity and, above a critical curvature, the interface becomes unstable. In addition, we observed an interfacial dancing state similar to the unconfined dancing state. In this state, +1/2 defects are continuously formed at the IN interface and are ejected into the bulk nematic. We also explored other channels and initial conditions to verify the robustness of the results.

To perform the simulations, we used an improved hybrid model based on the lattice Boltzmann method and on a predictor-corrector FD scheme. Among the improvements are: the use of a more isotropic lattice (D3Q19); a multi-relaxation-time collision operator, which is more accurate and stable since it allows the choice of different relaxation rates for the hydrodynamic moments; implementation of the stress tensor in the force term, eliminating spurious velocities. Those improvements were required to obtain a reliable description of the interfacial dynamics since these systems are prone to spurious numerical effects.

A final word on experiments. Although the motivation of the work reported here was theoretical, the relevance of interfaces in active nematic experiments is clear. As mentioned in the introduction the active turbulent state may be observed at temperatures above the passive NI transition rather than deep in the bulk nematic phase. Under these conditions, the role of interfaces may affect transient or steady dynamical states, as described here. Indeed, recently an experimental study of the propagation of active-passive interfaces in bacterial swarms has been reported 49. The methods and models proposed in our work may be used to address a number of questions raised by these and related experiments.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge financial support from the Portuguese Foundation for Science and Technology (FCT) under the contracts: PTDC/FIS-MAC/28146/2017 (LISBOA-01-0145-FEDER-028146) and UID/FIS/00618/2019. Margarida Telo da Gama (MTG) would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program ”The mathematical design of new materials” where most of this work was carried out. This program was supported by EP-SRC Grant Number: EP/R014604/1. MTG participation in the program was supported in part by a Simons Foundation Fellowship. We thank Prof. Mykola Tasinkevych for the fruitful discussions.

Notes and references

• Ramaswamy 2017

S. Ramaswamy, Journal of Statistical Mechanics: Theory and Experiment, 2017, 2017, 054002.

• Bechinger et al. 2016

C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Review of Modern Physics, 2016, 88, 045006.

• Cates and Tailleur 2015

M. E. Cates and J. Tailleur, Annual Review of Condensed Matter Physics, 2015, 6, 219–244.

• Doostmohammadi et al. 2018

A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nature Communications, 2018, 9, 3246.

• Marchetti et al. 2013

M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Review of Modern Physics, 2013, 85, 1143–1189.

• Doostmohammadi et al. 2016

A. Doostmohammadi, M. F. Adamer, S. P. Thampi and J. M. Yeomans, Nature Communications, 2016, 7, 10557.

• de Gennes and Prost 1995

P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, Clarendon Press, 1995.

• Beris and Edwards 1994

A. Beris and B. Edwards, Thermodynamics of Flowing Systems: with Internal Microstructure, Oxford University Press, 1994.

• Maitra et al. 2018

A. Maitra, P. Srivastava, M. C. Marchetti, J. S. Lintuvuori, S. Ramaswamy and M. Lenz, Proceedings of the National Academy of Sciences, 2018, 115, 6934–6939.

• Wu et al. 2017

K. 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, year.

• Opathalage et al. 2019

A. Opathalage, M. M. Norton, M. P. N. Juniper, B. Langeslay, S. A. Aghvami, S. Fraden and Z. Dogic, Proceedings of the National Academy of Sciences, 2019, 116, 4788–4797.

• Voituriez, R. et al. 2005

Voituriez, R., Joanny, J. F. and Prost, J., Europhys. Lett., 2005, 70, 404–410.

• Marenduzzo et al. 2007

D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Physical Review E, 2007, 76, 031921.

• Edwards, S. A. and Yeomans, J. M. 2009

Edwards, S. A. and Yeomans, J. M., EPL, 2009, 85, 18008.

• Shendruk et al. 2017

T. N. Shendruk, A. Doostmohammadi, K. Thijssen and J. M. Yeomans, Soft Matter, 2017, 13, 3853–3862.

• Giomi et al. 2012

L. Giomi, L. Mahadevan, B. Chakraborty and M. F. Hagan, Nonlinearity, 2012, 25, 2245–2269.

• Wensink et al. 2012

H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proceedings of the National Academy of Sciences, 2012, 109, 14308–14313.

• Chandragiri et al. 2019

S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, Soft Matter, 2019, 15, 1597–1604.

• Solon et al. 2018

A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri and J. Tailleur, Physical Review E, 2018, 97, 020602.

• Solon et al. 2018

A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri and J. Tailleur, New Journal of Physics, 2018, 20, 075001.

• Tjhung et al. 2018

E. Tjhung, C. Nardini and M. E. Cates, Physical Review X, 2018, 8, 031080.

• Thampi, Sumesh P. et al. 2015

Thampi, Sumesh P., Doostmohammadi, Amin, Golestanian, Ramin and Yeomans, Julia M., EPL, 2015, 112, 28004.

• Blow and da Gama 2013

M. L. Blow and M. M. T. da Gama, Journal of Physics: Condensed Matter, 2013, 25, 245103.

• West, J. L. et al. 2006

West, J. L., Zhang, K., Glushchenko, A., Andrienko, D., Tasinkevych, M. and Reznikov, Y., The European Physical Journal E, 2006, 20, 237–242.

• Blow et al. 2014

M. L. Blow, S. P. Thampi and J. M. Yeomans, Physical Review Letters, 2014, 113, 248303.

• De Magistris et al. 2014

G. De Magistris, A. Tiribocchi, C. A. Whitfield, R. J. Hawkins, M. E. Cates and D. Marenduzzo, Soft Matter, 2014, 10, 7826–7837.

• Blow et al. 2017

M. L. Blow, M. Aqil, B. Liebchen and D. Marenduzzo, Soft Matter, 2017, 13, 6137–6144.

• Fialho et al. 2017

A. R. Fialho, M. L. Blow and D. Marenduzzo, Soft Matter, 2017, 13, 5933–5941.

• Succi 2018

S. Succi, The Lattice Boltzmann Equation: For Complex States of Flowing Matter, Oxford University Press, 2018.

• Krüger et al. 2016

T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method - Principles and Practice, Springer International Publishing, 2016.

• Vesely 2001

F. Vesely, Computational Physics: An Introduction, Springer US, 2001.

• Blow and da Gama 2013

M. L. Blow and M. M. T. da Gama, Journal of Physics: Condensed Matter, 2013, 25, 245103.

• Hatwalne et al. 2004

Y. Hatwalne, S. Ramaswamy, M. Rao and R. A. Simha, Physical Review Letters, 2004, 92, 118101.

• Sengupta et al. 2013

A. Sengupta, U. c. v. Tkalec, M. Ravnik, J. M. Yeomans, C. Bahr and S. Herminghaus, Physical Review Letters, 2013, 110, 048303.

• Coelho et al. 2017

R. C. V. Coelho, M. Mendoza, M. M. Doria and H. J. Herrmann, Physical Review B, 2017, 96, 184307.

• Denniston et al. 2002

C. Denniston, G. Tóth and J. M. Yeomans, Journal of Statistical Physics, 2002, 107, 187–202.

• Denniston et al. 2004

C. Denniston, D. Marenduzzo, M. E. Orlandini and J. M. Yeomans, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2004, 362, 1745–1754.

• Coveney et al. 2002

P. V. Coveney, S. Succi, D. d’Humières, I. Ginzburg, M. Krafczyk, P. Lallemand and L.-S. Luo, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 2002, 360, 437–451.

• Yu and Fan 2010

Z. Yu and L.-S. Fan, Physical Review E, 2010, 82, 046708.

• Ammar et al. 2017

S. Ammar, G. Pernaudat and J.-Y. Trépanier, Journal of Computational Physics, 2017, 343, 73 – 91.

• Shan 2006

X. Shan, Physical Review E, 2006, 73, 047701.

• Wagner 2003

A. J. Wagner, International Journal of Modern Physics B, 2003, 17, 193–196.

• Coelho and Doria 2018

R. C. V. Coelho and M. M. Doria, Computers & Fluids, 2018, 165, 144 – 159.

• Thampi et al. 2015

S. P. Thampi, A. Doostmohammadi, R. Golestanian and J. M. Yeomans, EPL (Europhysics Letters), 2015, 112, 28004.

• Doi 2013

M. Doi, Soft Matter Physics, OUP Oxford, 2013.

• Croxton 1986

C. Croxton, Fluid interfacial phenomena, Wiley, 1986.

• Batista et al. 2015

V. M. O. Batista, M. L. Blow and M. M. Telo da Gama, Soft Matter, 2015, 11, 4674–4685.

• Chandragiri et al. 2019

S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, Soft Matter, 2019, 15, 1597–1604.

• Patteson et al. 2018

A. E. Patteson, A. Gopinath and P. E. Arratia, Nature Communications, 2018, 9, 5373.