Self-propulsion of an active polar drop

TL;DR

This paper analyzes how active polar drops spontaneously break symmetry and self-propel due to active stresses, revealing mechanisms and differences between contractile and extensile stresses, including chaotic motion under extensile stress.

Contribution

It provides an analytical theory explaining the active stress-induced instability and self-propulsion in polar drops, complementing previous numerical observations.

Findings

Both contractile and extensile stresses cause instability and propulsion.

Extensile stress leads to chaotic motion of the drop.

Analytical model explains symmetry breaking mechanisms.

Abstract

We investigate the self-propulsive motion of a drop containing an active polar field. The drop demonstrates spontaneous symmetry breaking from a uniform orientational order into a splay or bend instability depending on the types of active stress, namely, contractile or extensile, respectively. We develop the analytical theory of the mechanism of this instability, which has been observed only in numerical simulations. We show that both contractile and extensile active stress result in the instability and self-propulsive motion. We also discuss asymmetry between contractile and extensile stress, and show that extensile active stress generates chaotic motion even under a simple model of the polarity field coupled with motion and deformation of the drop.

Equations

$$\frac{\partial \mathbf{P}}{\partial t}+\left(\left(\mathbf{v}+v_p\mathbf{P}\right)\cdot \nabla \right)\mathbf{P}+\mathbf{\omega }\cdot \mathbf{P}$$

$${\kappa }_{ij}$$

$${\omega }_{ij}$$

$$\eta {\nabla }^2\mathbf{v}-\nabla p+\mathbf{f}$$

$$\mathbf{f}$$

$${\sigma }_{ij}^{(a)}$$

$${\sigma }_{ij}^{(e)}$$

$$\frac{\partial \mathbf{P}}{\partial t}+\omega \cdot \mathbf{P}$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{v}\cdot \nabla \varphi$$

$$=$$

$$\delta$$

$$F_{\mathrm{GL}}=$$

$$+\frac{\alpha g}{12}{\left(V_0-\int {\varphi }^2(3-2\varphi )d{V}^{\prime }\right)}^2.$$

$$F$$

$$F_{\mathrm{p}}[\mathbf{p}(\mathbf{x})]$$

$$F_{\mathrm{cp}}[\varphi (\mathbf{x}),\mathbf{p}(\mathbf{x})]$$

$$\tilde{\mathbf{v}}(\mathbf{k})$$

$${\tilde{T}}_{ij}(\mathbf{k})$$

$$\mathrm{Pe}$$

$$\mathbf{P}$$

$$\psi$$

$$\mathbf{v}$$

$${\partial }_t\psi ={\omega }_{xy}-\left(v_x{\partial }_x+v_y{\partial }_y\right)\psi$$

$$+\xi \left[\frac{{\kappa }_{yy}-{\kappa }_{xx}}{2}\sin 2\psi +{\kappa }_{xy}\cos 2\psi \right]+{\Gamma }^{-1}K\Delta \psi .$$

$$\mathrm{Er}$$

$$\Delta \psi =\mathrm{Pe}({\partial }_t\psi -{\omega }_{xy}+\mathbf{v}\cdot \nabla \psi$$

$$-\xi \left[\frac{{\kappa }_{yy}-{\kappa }_{xx}}{2}\sin 2\psi +{\kappa }_{xy}\cos 2\psi \right])$$

$$\Delta \mathbf{v}-\nabla p+\mathbf{f}=0$$

$$\mathrm{div}\mathbf{v}=0$$

$$\mathbf{f}=\mathrm{div}\left(\pm \mathbf{P}\mathbf{P}+\mathrm{Er}{\sigma }^{(e)}\right).$$

$$\mathbf{n}\cdot {\mathbf{v}}^{(i)}$$

$$\mathbf{t}\cdot {\mathbf{v}}^{(i)}$$

$$\mathbf{n}\cdot {\mathbf{\sigma }}^{(o)}$$

$$\mathbf{u}$$

$$v_i$$

$$F_{i_1{i}_2\cdots i_l}$$

Full text

I Introduction

Cell motility is a fascinating phenomenon in biological systems. Depending on the cell type, cells may move in a fluid, on a substrate, in a tissue consisting of other cells, or many other environments Bray (2000). It is evident that this phenomenon is a result of the complex cooperation of biomacromolecules and therefore it is difficult to understand its mechanism. Nevertheless, simple models may give an insight into the essential aspects of the mechanism responsible for cell motility, and assist in the classification of its mechanism into several universality classesCamley and Rappel (2017). This direction of research has been met with great success, for example, for micro-organisms in a fluid, which is described by a squirmer modelLighthill (1952); Blake (1971). The model imitates the motion of flagella and cilia on a cell body as an effective slip boundary on a solid (but possibly deformable) body, and all the complexity of propulsion falls into a functional form of the effective slip boundary condition of a fluid dynamics problem. Despite the reasonable generality of the model, the problem is analytically tractable for a simple geometry, and for this reason, extensive studies have been carried out to clarify that it can reproduce realistic motion near a wall, interaction between swimmers, and collective behaviorsLauga and Powers (2009).

The extension of the idea toward much more complex motility such as keratocyte or amoeboid movement on a substrate is in its infancy. Such motility is driven by the activity of a cytoskeleton, the mechanics of which is controlled by the collective behaviors of actin filaments and myosin molecular motors. The activity of a cytoskeleton falls into two classes: (i) actin polymerization and (ii) contractility of actin filaments mediated by molecular motors. The latter may be classified into two sub-categories: (ii-a) contractility on cortex on the surface of a cell and (ii-b) contractility in cytosol in the bulk. Recently, several models have been proposed to describe each mechanism such as polarity-driven cell motility Shao et al. (2010); Ziebert et al. (2012); Aranson (2015); L ber et al. (2015); Camley et al. (2017) for (i) 111 See also Callan-Jones et al. (2008); Blanch-Mercader and Casademunt (2013) discussing another mechanism of polarity-driven self-propulsion, namely, deformation-induced motion.

, and motility induced by active stress Tjhung et al. (2012); Recho et al. (2013); Tjhung et al. (2015) for (ii-b). The surface contractility results in blebbing, which is considered as another type of motility Charras et al. (2008); Kapustina et al. (2016).

The main ingredients of these models are the field that describes position and shape of the moving cell, and the polarity field that represents orientation and its magnitude of actin filaments inside the cytoskeleton. In this respect, they are regarded as active hydrodynamics confined in a drop whose interface is moving by the internal force generated by active fluidMarchetti et al. (2013). This combination of the polarity field and free-boundary problem is one of the main obstacles to analytical treatment. This is, in fact, another motivation for this study, namely, we would like to compare the models of cell motility and self-propulsion of a chemically driven dropYoshinaga (2017). In the latter model, a drop produces or consumes a concentration field such as surfactants by chemical reaction, resulting in inhomogeneous surface tension and motion of the drop through the Marangoni effectYabunaka et al. (2012); Yoshinaga et al. (2012); Yoshinaga (2014); Yoshinaga and Yabunaka (2018). This model shares the free-boundary problem of the drop with the model of contractility-driven cell motility, but the polarity field is replaced by a scalar field describing the concentration of chemical molecules. Owing to the simplicity of the scalar field, various analytical treatments have been proposedYabunaka and Yoshinaga (2016); Thomas Le Goff (2017); Seyboldt and J licher (2018). The polarity field is, on the other hand, very difficult to treat because of its vectorial nature. This drawback of its complexity may be compensated for by a richer structure in the model; as we will show, the active polar drop may demonstrate more various motion than the motion of chemically driven drops where only straight and helical motion have been observed. By analyzing the model with the polar field, we hope to clarify the similarity and differences between these models.

Active polar or nematic fluids exhibit distinct features compared from simple isotropic fluid and also from conventional liquid crystal, even under unconfined systems. First, the fluid may show topological defects (disclination), which plays an important role in its dynamics. Although the defects are also seen in conventional liquid crystals, the motion of the defect is strongly influenced by its topological charge (sign); for example, only the $1/2$ disclination of the active nematic fluid is able to move spontaneouslyPismen (2013); Giomi et al. (2013, 2014). Another important feature of the active fluid is that uniform orientation becomes unstable as activity increases. Eventually, it demonstrates turbulence-like flow containing a large number of topological defectsThampi et al. (2014). This instability occurs by coupling between fluid flow and the polarity field. It was first proposed theoretically in Simha and Ramaswamy (2002), followed by detail theoretical studiesVoituriez et al. (2005); Edwards and Yeomans (2009) as well as numerical analysis Giomi and Marchetti (2012). Depending on the sign of the active stress, contractile and extensile stress results in splay and bend instability, respectively, in a flow-tumbling regime. On the other hand, in a flow-aligning regime for a rod-like shape, only the extensile stress is linearly unstableEdwards and Yeomans (2009).

Intuitively, the active polar drop uses the mechanism of this instability to move spontaneouslyTjhung et al. (2012). Our claim is that this is true for the extensile stress, whereas for the contractile stress, the effect of the interface confining the active fluid inside the drop must be considered. We discuss how the interface gives rise to additional flow. We also show that after the instability, the distortion of the polarity field produces source dipole flow resulting in self-propulsion.

In this work, we focus on active polar drops where the motility is driven by bulk contractility generated by active stress. To the best of the author’s knowledge, most previous studies have focused only on numerical simulations in two Tjhung et al. (2012); Giomi and DeSimone (2014); Marth et al. (2014); Fialho et al. (2017) or three dimensions Tjhung et al. (2015). There have been scarce theoretical studies on active drops. In A. Whitfield et al. (2014), the speed of an active polar drop is computed by assuming an ansatz of the polar field A. Whitfield et al. (2014). The flow field inside a thin film of a drop was analytically calculated in Khoromskaia and Alexander (2015). These works focused on fluid flow and the speed of a drop under a given polarity field designed for their analysis. However, it is not clear how transition (bifurcation) occurs between the stationary state and the self-propelled state. The main focus in this work is to investigate the mechanism of the transition (drift bifurcation) using as a simple setting as possible.

II Model

We consider an active polar field composed in a drop with a radius $R$ in two dimensions. The model for polarity distribution, ${\bf P}({\bf x})$, is written asTjhung et al. (2012)

[TABLE]

where the symmetric and anti-symmetric parts of the velocity gradient tensors are denoted by $\kappa_{ij}({\bf x})$ and $\omega_{ij}({\bf x})$, respectively. They are explicitly given by the gradient of the velocity field ${\bf v}({\bf x})$ as

[TABLE]

The term $\xi\bm{\kappa}\cdot{\bf P}$ is called flow alignment, and imitates the shape of a filament such as rod-like ($\xi>0$) or disk-like $\xi<0$. It also demonstrates shear-alignment ($|\xi|>1$) or shear-tumbling ($|\xi|<1$) depending on the parameter $\xi$. Here, $v_{p}$ is the polymerization velocity and describes advection. We neglect this term here and set $v_{p}=0$ because this is not the main mechanism of self-propulsion in this model.

As we are interested in small systems with the characteristic length scale $l\simeq 0.1-10\mu$ m, we use the Stokes equation:

[TABLE]

where the force acting on the fluid is

[TABLE]

The force acting on the fluid has two parts: active stress and elastic stress. The active stress is given by Kruse et al. (2005); Tjhung et al. (2012)

[TABLE]

The sign of the parameter $\zeta$ demonstrates contractile stress ($\zeta>0$) and extensile stress ($\zeta<0$). The elastic stress $\sigma^{(}e)$ arises from the Frank elasticity of the polar field and under a constant approximation

[TABLE]

where $\bf h=-\delta F/\delta{\bf p}$.

As we are interested in the universal aspects of the model, we simplify it so that it can reproduce spontaneous motion. It is often the case that friction of the velocity field is included in (4) to express friction between a cell and substrate. This term is given by $-\kappa{\bf v}$, which introduces another length scale $\sqrt{\eta/\kappa}$ into the model. The length scale sets screening of the hydrodynamic flow, and when the friction is stronger, the velocity gradient is more localized near the interface between the drop and surrounding fluid. Although this term changes a critical activity to obtain self-propulsion, it does not change the structure of bifurcation and therefore it is not the main mechanism of motility in this model. We, thus, neglect this term to consider the minimal ingredients of self-propulsion.

The free energy is chosen as (13) so that the polar state $|{\bf P}=1|$ is stable and thus the equation of the polarity field is given by

[TABLE]

Together with (9), we have the closed equations.

II.1 Numerical Simulation

To consider translational motion and deformation of the drop, the density field $\phi({\bf x})$ is introduced using the phase-field approach

[TABLE]

where

[TABLE]

The equation has the two stable points $\phi=0$ and $\phi=1$, and they indicate inside ($\phi=1$) and outside ($\phi=0$) the drop. Although the conventional Ginzburg-Landau (GL) equation has an infinitely growing domain, the modified version suppresses the growth to describe a drop with a finite size. When the volume is equal to $V_{0}$, $\phi=1/2$ is an unstable point as in the standard time-dependent Ginzburg-Landau (TDGL) model. The difference is found in $\delta(\phi)$, which controls the unstable point so that the total volume is approximately conserved as $V_{0}$. The advantage of this approach is that the model has a free energy Nonomura (2012)

[TABLE]

The additional term makes a contribution to the hydrodynamic equation as a gradient term, and this does not modify the form of the force due to the phase field.

The free energy of this model is given by the following three terms

[TABLE]

The free energy of the polarity field ensures that the amplitude of the polarity is $|{\bf P}|=1$ almost everywhere except in the region in which the polarity field is largely deformed. The deformation of the polarity field is penalized by the Frank elastic constant, $K$. The indices in the elastic term assumes a matrix norm. The coupling between $\phi$ and ${\bf p}$ is chosen so that $|{\bf p}|=1$ inside the drop where $\phi=1$ and $|{\bf p}|=0$ outside. For this purpose, we set $\beta=1$.

The velocity field is calculated in the Fourier space. The velocity field is expressed as

[TABLE]

where the Oseen tensor in the Fourier space is

[TABLE]

III Spontaneous Motion and Deformation

First, we discuss our numerical results. The active polar drops are stationary when their activity is low and their size is small. As the activity increases, the stationary state becomes unstable and the drop starts to move. The motion is straight with a constant speed as shown in Fig. 1(A). This motion occurs both in extensile and contractile drops. This result has already been obtained in Tjhung et al. (2012). However, as the activity increases, the extensile drop produces different motion to the contractile drop. The extensile drop has a spinning motion where the center of mass does not move, but the polarity field changes by rotating the drop (Fig. 1(C)). Rotational motion also occurs in which the translational motion follows a closed path (Fig. 1(A)). Then, as the activity increases, the motion of the extensile drop changes from a zigzag motion to random motion (Fig. 1(B)). These complex motions do not occur in the contractile drop. Such chaotic motion has not been obtained by the hydrodynamic model, but has been observed in extensile drops using the kinetic modelGao and Li (2017).

The different behaviors of spontaneous motion between the contractile and extensile drops are evident in the phase diagram shown in Fig. 2. The overall tendency is that for higher activity and larger drops, the self-propulsive motion becomes complex. The transition between the stationary state and the straight translational motion is characterized by a non-dimensional number of the Péclet number, Pe,

[TABLE]

In fact, the transition between stationary and self-propelled states occurs at ${\rm Pe}\simeq 25$ as shown in Fig. 2. At higher activity, the boundary between two types of the motion becomes less clearer. In the phase diagram, we define each state as stationary for $\langle|u|\rangle\leq 0.001$ and $\langle|\omega|\rangle\leq 0.01$, translation for $\langle|u|\rangle>0.001$ and $\langle|\omega|\rangle\leq 0.01$, spinning for $\langle|u|\rangle\leq 0.001$, $\langle|\omega|\rangle>0.01$.When $\langle|u|\rangle>0.001$ and $\langle|\omega|\rangle>0.01$, the motion is called zigzag for $({\rm max}\omega-{\rm min}\omega)/{\rm max}|\omega|\geq 1$ and rotation for $({\rm max}\omega-{\rm min}\omega)/{\rm max}|\omega|<1$. This criteria ensures whether the drop periodically turns left and right (zigzag), or moves with constant angular velocity (rotation). When the speed and angular velocity changes irregularly, its state is called chaotic.

To summarize the numerical results, several observation can be made. (i) Spontaneous motion occurs irrespective of the sign of active stress. The contractile and extensile drops both exhibit self-propulsion, and seem to be unstable above the critical activity $\zeta_{c}$. For an extensile drop, this is not surprising as the bend instability should occur in unconfined systems for $\xi>1$. However, the contractile active polar fluid in large systems is linearly stableEdwards and Yeomans (2009), while both our numerical results and those in Tjhung et al. (2012); Marth et al. (2015) indicate the pitchfork bifurcation. (ii) Both the hydrodynamic modelTjhung et al. (2012); Marth et al. (2015) studied here using partial differential equations of density and polarity fields and the kinetic model Gao and Li (2017) using a Smoluchowski equation of a probability distribution of a position and orientation of filaments reproduce self-propulsion including chaotic motion. This implies that the mechanism of various types of self-propulsion is encoded in the hydrodynamic model, while to have quantitative features, one has to consider the kinetic model, which includes more information about higher-order moments.

IV Theoretical Analysis

To study the mechanism of self-propulsion, we consider the polarity field inside a disk with a radius $R$. When there is no disclination, the polarity is a unit vector and thus is expressed by

[TABLE]

The phase $\psi({\bf x})$ may be expanded in polar coordinates as

[TABLE]

The advantage of this approach is that the bend (splay) instability is expressed by $m=1$ and $\psi=Ar$ ($\tilde{\psi}=Ar$). It is also convenient to expand other vector fields such as velocity ${\bf v}$ and force ${\bf f}$ fieldsYoshinaga and Liverpool (2018), for example,

[TABLE]

where ${\bf n}$ and ${\bf t}$ are unit normal and tangential vectors, respectively (see Fig. 4), and we define two unit vectors as ${\bf r}_{m}\equiv\left(\cos m\theta,\sin m\theta\right)$ and ${\bf t}_{m}\equiv\left(-\sin m\theta,\cos m\theta\right)$. For the special case, ${\bf r}_{1}={\bf n}$ and ${\bf t}_{1}={\bf t}$.

By projecting (1) onto $(-\sin\psi,\cos\psi)$, we obtain the equation of the phase as

[TABLE]

We choose the unit length scale as $R$, the unit time scale as $R/U_{0}$, the characteristic velocity as $U_{0}=\frac{\zeta R}{\eta}$, and the characteristic pressure as $p_{0}=\zeta$. The non-dimensional numbers in this system are the Peclet number (17) and the Ericksen number

[TABLE]

The non-dimensionalized equations are obtained as

[TABLE]

Here, the sign in the active stress is chosen such that positive ($\zeta>0$) corresponds to contractile stress whereas negative ($\zeta<0$) corresponds to extensile stress. Without topological defects, the amplitude of the polarity field is almost constant $|{\bf P}|\simeq 1$ and therefore the elastic stress is $\sigma^{(e)}\sim\left(\partial{\bf P}\right)\left(\partial{\bf P}\right)$. This term is neglected in our analysis as it is not necessary to reproduce self-propulsion. We have confirmed that the numerical simulations without elastic stress also demonstrate self-propulsion of the active polar drop.

Boundary conditions for (24)-(25) at $r=1$ are

[TABLE]

Here, ${\bf u}$ is self-propulsive velocity. The first two equations demonstrate continuity of the velocity fields inside and outside the drop. The last condition implies force balance on the surface of the drop. There is another boundary condition for $\psi$ in (23). This is set by the condition that there is no mechanical force acting in the system (force-free condition) in (68).

IV.1 Self-propulsive Motion

In this section, we consider how self-propulsion occurs in the active polar drop. The velocity of a drop is computed asYabunaka et al. (2012); Yoshinaga (2014)

[TABLE]

where ${\bf R}$ is a vector pointing to the surface of the drop. The normal velocity on the surface is $v_{n}={\bf v}\cdot{\bf n}$. For a circular drop, the velocity given by (30) is non-zero only when $v_{n}\sim\cos\theta$ or $v_{n}\sim\sin\theta$. We expand the velocity field using force multipoles, $F_{i_{1}i_{2}\cdots i_{l}}$, as

[TABLE]

where $T_{ij}({\bf x},{\bf x}^{\prime})$ is the Oseen tensor, and the $l$th force multipoles are expressed as

[TABLE]

Under the force-free condition, the self-propulsive velocity (30) is nonzero only when there is a source dipole. The velocity field of the source dipole in two dimensions is obtained by taking trace ($\delta_{kl}$) of $\partial_{k}\partial_{l}T_{ij}$ as

[TABLE]

where the source dipole is obtained from the distribution of the force as

[TABLE]

For a given source dipole, the velocity is obtained as

[TABLE]

This aspect of the self-propulsion is shared by other phenomena such as squirmersLauga and Powers (2009), phoretic motion and self-phoresisJülicher and Prost (2009), and self-propulsion driven by the Marangoni effectYoshinaga (2014). All these motions are driven by the source dipole of their flow fields. In contrast to the active polar drop, the flow field in these models is driven by interfacial (or surface) forceAnderson (1989). In the active drop studied here, the flow is driven by active stress acting in bulk. Nevertheless, in terms of force-free motion, only the source dipole is associated with self-propulsion.

The active stress is expressed in terms of polynomials of $\psi$

[TABLE]

We subtract a trace of the active stress as it merely modifies the pressure. By expansion of (19), the self-propulsive velocity is also expanded in terms of its coefficients $\psi_{m}$ and $\tilde{\psi}_{m}$ as

[TABLE]

When $\psi(r)=a_{1}^{1}r$ and $\tilde{\psi}(r)=b_{1}^{1}r$ with a given constants $a_{1}^{1},b_{1}^{1}$, the velocity is ${\bf u}\sim\pm(b_{1}^{1},a_{1}^{1})$. This analysis suggests that, for small distortion of the polarity field, the velocity of the drop is proportional to the first mode of the expansion in (19). The contractile active stress gives rise to self-propulsion in the $+x$ direction for $\tilde{\psi}_{1}>0$, whereas the extensile active stress results in motion in the $-y$ direction for $\psi_{1}>0$ (see Fig. 3).

IV.2 Perturbation around Stationary States

When ${\rm Pe}=0$, the stationary polarity field is uniform inside the drop. As ${\rm Pe}$ increases, the flow field is generated by the active stress. When ${\rm Pe}$ is small, the flow may perturb the polarity field, but it is not strong enough to generate a source dipole. The flow is dipolar and does not lead to self-propulsion. The drop remains stationary. To study the mechanism of self-propulsion, we need to clarify when this stationary state becomes unstable. Therefore, we linearize the model around the stationary state.

In contrast to self-propulsion driven by chemical reactions through the Marangoni effect Yabunaka et al. (2012); Yoshinaga (2014) in which the stationary state is trivial, the stationary state of the active polar drop is not simple owing to the flow alignment term. The uniform polarity field does not satisfy (23) owing to the stress acting on the surface of the drop. Nevertheless, it is important to observe that (19) is divided into two terms: odd and even $m$. The stationary state results in dipolar flow, which generates perturbation only in even $m$ terms. This implies that when the initial phase does not contain odd modes, there is no odd mode in the final state. On the other hand, self-propulsion is associated with the first mode of (19). As we are interested in transition between the stationary state and self-propelled states, it is natural to use a following scaling

[TABLE]

We denote by $\psi^{*}$ a stationary solution of (23), which contains only even modes in (19). We linearize a set of the model (23) in terms of $\epsilon$. Each variable is expanded as

[TABLE]

and following (19), each order in (39) is expanded in modes denoted by $m$. Note that $\psi^{*}$ contains only even $m$ whereas $\psi^{(1)}$ contains $m=1$ and all even $m$. Higher-order terms such as $\psi^{(2)}$ contain, in general, all $m$. The active stress is expanded accordingly as

[TABLE]

At $\mathcal{O}(\epsilon)$, the model becomes,

[TABLE]

where the velocity gradient and rotation tensors are proportional to $\psi^{(1)}$ as $\kappa_{ij}^{(1)}=\frac{1}{2}\left(\partial_{i}v_{j}^{(1)}+\partial_{j}v_{i}^{(1)}\right)$ and $\omega_{ij}^{(1)}=\frac{1}{2}\left(\partial_{i}v_{j}^{(1)}+\partial_{j}v_{i}^{(1)}\right)$. The velocity field perturbed by $\psi^{(1)}$ is the solution of (24) under the active stress (42). We are interested in $\psi^{(1)}_{1}$: when it becomes non-zero as ${\rm Pe}$ is varied and how it grows. Thanks to the scaling (38), the even modes in $\psi^{(1)}$ do not generate $\psi^{(1)}_{1}$ because $\psi^{*}$ consists only of even modes. Therefore, (43) is formally rewritten as

[TABLE]

We analyze (44) by using the Zernike expansion of $\psi({\bf x})$ is given by

[TABLE]

where

[TABLE]

for even $n-m$. For odd $n-m$, then $R_{n}^{m}=0$. Note that $n\geq m$. Clearly, for given $n\in[0,\infty)$ and $n\geq m$, $R_{n}^{m}(r)$ is expressed by power series in $r$.

Taking inner product of the linearized equation (44) with the Zernike polynomials for $m=1$ in the expansion of (45), we obtain the following linear algebraic equation using the boundary condition (68):

[TABLE]

where each block is a square matrix and is given by

[TABLE]

and

[TABLE]

The other two matrices are expressed in a similar way. The first term in (48) corresponds to the left-hand side in (44) owing to the Frank elasticity (see (82)), whereas the second term corresponds to the right-hand side from the coupling with fluid flow. The first row in the matrices represents the boundary condition (68), namely,

[TABLE]

Other rows in the matrices are obtained by taking an inner product of (44) for each basis in (45), namely, multiplying $R_{n}^{1}(r)$ and integrating over $r$ using the orthogonal relation (81). From the second row, the projection onto $R_{1}^{1}(r)$, $R_{3}^{1}(r),\ldots$ is performed. When ${\rm Pe}=0$, the matrix (48) is clearly invertible and therefore the solution is $a_{n}^{1}=b_{n}^{1}=0$ and the stationary solution is stable. The diagonal term in (48) is given by $2(k+2)$ for $k=1,3,\ldots$ corresponding to the second, third, $\ldots$ row. This suggest that the possible scenario for zero eigenvalue arises from the smallest $n$.

The stability of the stationary state $\psi^{*}$ is obtained by the determinant of the linearized matrix in (47). A non-trivial solution appears only when the determinant vanishes, and the matrix is not invertible. This corresponds to find zero eigenvalues in which the stationary state looses it stability. This aspect is explicitly shown in Appendix B.

IV.3 Stability around Uniform Orientation

To analyze (44), or equivalently (47), we need to calculate stationary solution of (23)-(25). The concrete form of $\psi^{*}$ is available only numerically. Nevertheless, it is reasonable to use an approximation

[TABLE]

where the polarity field is aligned with the $x$-axis (see Fig. 4). Here, the lowest-order term is chosen to be zero, though any spatially uniform constant is possible by rotation. Then, the inhomogeneous term is approximated as

[TABLE]

Note that when $\psi^{*}=\pi/2$ is chosen

[TABLE]

and, as we show later, the sign of the dipolar flow becomes opposite. As a result, the stability of (52) does not change. The physical meaning of (52) is decomposed into two parts. The first part consists of the first and last terms, which describe the flow generated by distortion of the polarity field. The second part consists of the second and third terms. Flow is generated by the stationary polarity field (51), which is coupled with the distortion of the polarity field and stabilizes/destabilizes the polarity field.

First, we compute steady velocity ${\bf v}^{*}$ and resulting shear flow $\kappa^{*}$. The active stress for the polarity field (51) is expressed as

[TABLE]

where $\Theta(x)$ is a step function $\Theta(x)=1$ for $x>0$ and $\Theta(x)=0$ otherwise. The positive and negative signs denote contractile and extensile active stress, respectively. The active force generated by the active stress is obtained as

[TABLE]

The steady velocity field is given by

[TABLE]

where

[TABLE]

Owing to incompressibility and the boundary condition of tangential stress at the interface, we obtain

[TABLE]

Note that the balance of normal stress at the interface is satisfied by deformation. We assume the uniform surface tension is large so that deformation is negligibly small and therefore (58) approximates the velocity field under deformation. Schematic pictures of the velocity field are given in Fig. 4. The extensile active stress generates dipolar flow pulling from the sides of the drop and pushing from top and bottom. The contractile active stress generates the flow in the opposite direction. The flow of contractile active stress is qualitatively the same as was observed in Tjhung et al. (2012). We finally obtain

[TABLE]

and

[TABLE]

Next, we consider a velocity field generated by distortion of the polarity field $\psi^{(1)}_{1}$. The active stress generated by $\psi^{(1)}_{1}$ is expressed as

[TABLE]

and the resulting force is

[TABLE]

with $\bar{\nabla}=(\nabla_{y},\nabla_{x})$. The distortion is expressed by the first mode in (19) which gives rise to the force .

[TABLE]

where

[TABLE]

In order to ensure the force-free condition, it requires

[TABLE]

We may express the velocity field at the first mode in (20), and calculate it by the method outlines in Appendix C. Once we obtain ${\bf v}_{r}^{(1)}$ and ${\bf v}_{t}^{(1)}$, the shear and rotational flow is expressed as

[TABLE]

Here, we have used incompressibility to simplify the rotational flow. Using the force generated by active stress (67), the $\cos\theta$ terms in (69) and (70) result in $L^{(a,a)}$ in the linearized matrix (47) for the coefficients of the expansion (45), whereas $\sin\theta$ terms result in $L^{(b,b)}$.

Owing to the effect of rotational flow and shear alignment, the determinant of the linearized matrix (47) becomes zero when the Péclet number becomes the critical Péclet number, ${\rm Pe}_{c}$. This instability occurs both for $a_{n}^{1}$ and $b_{n}^{1}$, suggesting that both splay and bend instability is induced by this effect. The physical meaning of the effect of the dipolar flow induced by the stationary polar field is sketched in Fig. 5. When the active stress is contractile, the orientation along the $x$-axis is unstable, and a small perturbation of the first mode $m=1$ in Fig. 3 will grow. The extensile active stress exhibits an opposite flow, and stabilizes orientation along the $x$-axis. The flow generated by the perturbed polarity field also gives rise to perturbation of the first mode as (69) and (70). The stability of the uniform polarity field is shown in Fig. 6. This is obtained by zeroth of the determinant of (47) with (63), (69), and (70). We truncate the Zernike expansion of $r$ in (45) at the lowest order necessary to find the transition. The advection term is found to be small and is neglected. As activity increases, instability occurs both for contractile and extensile stress in the rod-like flow-alignment regime ($\xi>1$).

We may also consider the system without boundary. In this case, (23)-(25) are linearized around the stationary state ${\bf P}^{*}=\hat{\bf x}$, that is, $\psi^{*}=0$. Note that in contrast to the active drop, there is no flow ${\bf v}^{*}=0$ at the stationary state of the bulk system. In the Fourier space, the system becomes unstable when

[TABLE]

where $\theta_{k}$ is the angle of the wave vector ${\bf k}$, and $L$ is the system size. The bend instability corresponds to $\theta_{k}=0$ whereas the splay instability corresponds to $\theta_{k}=\pi/2$. The result is shown in Fig. 6. Similar to the nematic caseEdwards and Yeomans (2009), the splay instability occurs when $\xi<1$ for the contractile stress whereas the bend instability occurs when $\xi>-1$. The mechanism of this asymmetry between contractile and extensile stress is that the rotational flow destabilizes the uniform orientation both by contractile (splay instability) and extensile (bend instability) stress, whereas the shear flow in the rod-like flow-alignment regime suppresses the instability only for the contractile stress. Our result suggests that the effect of the boundary enhances the instability for the contractile stress for $\xi>1$.

V Chaotic Motion and Topological Defects

When the activity is high, the self-propulsive motion is no longer periodic. Its long-time behavior is diffusive whereas the short-time behavior is ballistic, as shown in root mean-square displacement $\sqrt{\rm MSD}=\sqrt{\langle|{\bf x}(t)-{\bf x}(0)|^{2}\rangle}$ in Fig. 7(C). During the diffusive motion, distortion of the polar field is accumulated at lines, which are dynamically deforming (Fig. 7(B)). The shape of the drop is also fluctuating correlated with the motion of the line defects. As activity is increased, the density of the line defects increases, and at high extensile activity, disclinations appear.

In the case of an active nematic drop, the polarity field around a topological defect is expressed by

[TABLE]

where $\theta=\tan^{-1}y/x$. The plus and minus signs denote $+1/2$ and $-1/2$ disclinations, respectively. The active stress is then given by

[TABLE]

where we subtract the trace of the active stress because it merely modifies pressure in incompressible systems.

The force dipole in (32) generated by the active stress can be obtained by integrating (73) inside the drop. This clearly vanishes as $F_{ij}=0$. The source dipole of $+1/2$ disclination is from (34)

[TABLE]

and for the $-1/2$ defect the source dipole vanishes ${\bf q}=0$. As the source dipole is associated with the self-propulsive velocity, the $+1/2$ disclinations move whereas $-1/2$ disclinations do not. The same conclusion was reached by directly solving the Stokes equation in Pismen (2013); Giomi et al. (2014). Our analysis decomposes the self-propulsion into multipoles of forces and associated flow.

In the case of defects of the polar field, the field is expressed as

[TABLE]

After a calculation similar to that mentioned previously, we find there is neither force dipole nor source dipole for the defects.

The position of the topological defect is extracted by the method in Mazenko (1997). The signed defect density is obtained as

[TABLE]

At the position of defects with positive topological charge $\rho\gg 0$, whereas for negative topological charge $\rho\ll 0$. The number density of the defects with positive and negative charges is statistically the same. In Fig. 7(D), the density of the positive defects is shown. As activity increases, the density linearly increases.

VI Discussions and Summary

In this work, we have analyzed the active polar drop, and showed that it exhibits different type of self-propulsion depending on whether stress is contractile or extensile. The contractile stress results in translational motion, and, at higher activity, rotational and zigzag motion. In addition to these motions, chaotic diffusive motion occurs only with extensile active stress. The origin of the complex motion seems to be turbulence-like behaviors in bulk owing to the active stress, which has been discussed in the system without boundaries or with solid boundaries. On the other hand, the self-propulsion under contractile stress has a different origin: the flow generated by a surface force similar to inhomogeneous surface tension in the Marangoni effect. This extra effect, only existing in the confined drop, gives rise to instability of the uniform polarity field. This is in contrast to bulk systems where the contractile active stress of rod-like molecules does not linearly destabilize the uniform polar field.

Our result of asymmetry between the contractile and extensile stresses has some similarities to the results in Ramaswamy and J licher (2016), in which a contractile active polar fluid with frictional boundaries exhibits spontaneous flow and oscillatory dynamics in the flow-alignment regime, whereas spontaneous flow, oscillatory dynamics, and chaotic flow in the flow-tumbling regime. Although there is no interfacial effect in that work and they studied two different systems with different flow-alignment parameters, not activity, they observed chaotic flow only in one regime. It would be interesting as a future work to study nonlinear effects and the mechanism of chaotic flow.

Oscillatory dynamics such as traveling waves do appear as a secondary bifurcation, not as a linear instability through Hopf bifurcation. This is in contrast to the observations of an active nematic fluidGiomi et al. (2011, 2012) and active polar fluid with additional active stress in the form of $\sigma^{(a)}_{ij}\sim\zeta_{2}(\partial_{i}p_{j}+\partial_{j}p_{i})$ (see Giomi et al. (2008)). In these cases, the systems have additional time scales associated with the dynamics of the concentration field and/or velocity field (inertia term). Our system is in the low-Reynolds-number regime, and in the uniform concentration and thus lacks those time scales. Nevertheless, our system shows oscillatory dynamics such as zigzag motion after several transitions from the stationary state.

Recently, several experimental systems have been proposed to study self-propulsion of a drop containing liquid crystals with certain activityThakur et al. (2006); Krüger et al. (2016); Yamamoto and Sano (2017). Our model is probably too crude to explain the motion in these works, but we hope to convey a basic understanding of these phenomena.

Acknowledgements.

The authors are grateful to Rhoda Hawkins and Igor Aranson for helpful discussions. The authors acknowledge the support by JSPS KAKENHI Grant Nos. JP16H00793 and 17K05605.

Appendix A Zernike polynomials

Several lowest-order terms in Zernike polynomials are given by

[TABLE]

Orthogonality is expressed by

[TABLE]

The Laplacian acting on the Zernike polynomials results inJanssen (2014)

[TABLE]

This operator makes sense when the mapping is $N+2$ dimensions onto $N$ dimensions. Note that $R_{n}^{m}=0$ for odd $n-m$ and therefore the Laplacian operator induces the mapping from $[(N+2)/2]$ to $[N/2]$.

Appendix B Stability in Dynamics

The analysis in the previous section is associated with the stability of the steady state. To see this, the time evolution of the linearized equation is expressed, similar to (44), as

[TABLE]

Together with the boundary condition, we use the ansatz of $a_{n}^{1}(t)=a_{n}^{1}e^{\sigma t}$ and $b_{n}^{1}(t)=b_{n}^{1}e^{\sigma t}$ as

[TABLE]

where

[TABLE]

From $\det(L-\Lambda)=0$, the stability of the linearized equation is obtained from $\sigma$. When we neglect the second term on the right-hand side of (83), all the $\sigma$ are negative, and this suggests that the system relaxes to uniform orientation without the coupling to the shear and rotational flow.

Appendix C Velocity Field under Force

We solve the Stokes equation in two dimensions in the form of

[TABLE]

We will consider only the first mode, which is expressed by $m=1$ in (20). The pressure is given by

[TABLE]

With this expansion, incompressibility implies

[TABLE]

The Stokes equation is then rewritten as

[TABLE]

where

[TABLE]

The solution of the equation of the pressure is

[TABLE]

with the integral constant ${\bf A}$ to be determined by the boundary conditions. The solution of the velocity field is obtained as

[TABLE]

where ${\bf B}$ is another integral constant to be determined from the boundary condition, and $\bm{\psi}_{1}^{(1)}=(\tilde{\psi}_{1}^{(1)},\psi_{1}^{(1)})$. From the boundary condition, $({\bf v}_{r}(R)-{\bf u})\cdot{\bf n}=0$. Using the self-propulsive velocity ${\bf u}$, the tangential velocity at the surface is $({\bf v}_{t}-{\bf u})\cdot{\bf t}=v_{s}$, where the slip velocity, $v_{s}$ is associated with self-propulsive velocity ${\bf u}=-(1/2\pi)\int v_{s}{\bf t}d\theta$Yoshinaga and Liverpool (2018). The velocity field is then expressed as

[TABLE]

Once the tangential slip velocity is known, we may calculate the velocity field outside the drop, and, assuming same viscosity outside the drop, the boundary condition (29) results in

[TABLE]

In order to calculate (69)and (70), it is suffice to compute ${\bf v}^{\prime}_{n}+{\bf v}^{\prime}_{t}$. Using incompressibility and the result of the general solution, we obtain

[TABLE]

Using the recurrence relation of the Zernike polynomials

[TABLE]