Self-propulsion of active droplets without liquid-crystalline order

The swimming of cells, far from any boundary, can arise in the absence of long-range liquid-crystalline order within the cytoplasm, but simple models of this effect are lacking. Here we present a two-dimensional model of droplet self-propulsion involving two scalar fields, representing the cytoplasm and a contractile cortex. An active stress results from coupling between these fields; self-propulsion results when rotational symmetry is spontaneously broken. The swimming speed is predicted, and shown numerically, to vary linearly with the activity parameter and with the droplet area fraction. The model exhibits a Crowley-like instability for an array of active droplets.

Active fluids are an emerging class of non-equilibrium systems, where energy is injected into the system locally and continuously, by the constituent particles themselves [1].Many examples of active fluids are biological in nature, for example, actomyosin networks inside the cell cytoskeleton [2][3][4] and dense suspensions of microtubules and kinesins in vitro [5,6].In the case of actomyosin networks, each myosin motor can attach (and detach) to two actin filaments and pull the two filaments inwards, causing a net local contractile stress, which drives the system out-of equilibrium.In many cases, this local energy injection at the filament scale can be translated into a macroscopic motion.For example, actomyosin contraction at the rear of the cell cortex has been shown to play an important role in the swimming motility of cells in a bulk fluid environment [3,4], far from any boundary at which crawling can instead occur.
At the level of phase-field modeling and simulations, cell motility is often described as the spontaneous motion of a droplet of active fluid [7][8][9][10].The current field-theoretic understanding of cell swimming involves a scalar field φ, coupled to polar or nematic liquidcrystalline order, described by a vector or a second-rank tensor [7,8,11].The scalar field delineates the cell's interior (φ > 0) from its exterior (φ < 0) whereas the vectorial/tensorial field describes bulk internal alignment of a uniform or cortical 'cytoskeleton'.The propulsion mechanism then relies on a discrete broken symmetry along a pre-existing axis of orientational order, such as a spontaneous splay transition [11], or self-advection caused by net polymerization at the leading end of each polar filament [7,10,12].Thus, the vectorial/tensorial nature of the order parameter is crucial to obtain self-propulsion in such theories.On the other hand, experimental observations of cell swimming suggests direct rotational symmetry breaking of the actomyosin concentration delineating the cell cortex [3,4], implying that liquid-crystalline order is not a pre-requisite for self-propulsion in cells.
In this Letter, we present a field theory of selfpropelling active droplets, such as cells, in the absence of liquid-crystalline order.Our theory is given in terms of two scalar fields φ and ψ, which follow two separate (a) Flow driven by an active stress Σ A on the interface of the droplet phase field φ.A positive (negative) Σ A implies an extensile (contractile) fluid flow tangential to the φ-interface.(b) A cell is represented by φ-droplet, while concentration of actomyosin in the cell cortex given by φ−ψ (red ring).Under a small random initial perturbation, the ψ > 0 domain is displaced to the front (right), creating excess contractile stress at the back.This leads to an active flow which squashes the ψ-droplet in the front, causing sustained self-propulsion along a direction selected by spontaneous symmetry breaking.
Cahn-Hiliard equations, coupled by a momentum conservation or Stokes equation.The φ-field delineates as usual the interior/exterior of the droplet, whereas ρ = φ−ψ describes the presence of a cortex surrounding the droplet (a ring in 2D), rendered contractile by, e.g., the action of myosin motors.Self-propulsion of the active droplet emerges generically by the direct breakdown of rotational symmetry in the ρ field.Once this is broken, the droplets moves with a speed that increases linearly with the activity strength and decreases linearly with the droplet area fraction within a periodic domain.(We rationalize both scalings analytically.)We then extend our theory to study an array of active droplets and identify an active Crowley-like instability.
Eq. ( 1) is adopted as the simplest way to stabilize two concentric phase-field domains (in φ and ψ) of unequal size.On identification of ρ = φ − ψ, this becomes a droplet surrounded by a cortex, as required.Parameters are chosen so that each phase field approaches ±1 in the interior/exterior bulk phases.The coupling term, with β > 0, is chosen to make the φ and ψ domains sit on top of each other.Since both fields are conserved, the volumes of these φ, ψ droplets (equivalently, the φ droplet and its φ − ψ cortex) are separately constant in time.We choose these volumes so that in equilibrium the ψ-droplet resides within the φ-droplet giving a cortex in between, see Fig. 1(b).The φ-droplet then has interfacial width ξ 0 −κ/2a and tension γ 0 −8κa 3 /9b 2 [13], with similar expressions for the ψ-droplet.These tensions govern respectively the cortex/exterior and cortex/interior interfaces.Note that the β-term also renormalizes the interfacial width ξ 0 and tension γ 0 ; however, we choose β small so that this difference is not appreciable.
The only active term in our model is a contractile stress, which lives, for simplicity, at the outer interface where φ passes through zero.This active stress is, however, modulated by the local concentration of cortical material which (in some units) is ρ(r, t) = φ − ψ −ψ > 0 (see Fig. 2(a)).If rotational symmetry is maintained, the cortex is concentric with the droplet and the active stress is likewise symmetric (Fig. 1(b), left).In the broken symmetry state, with more cortical material at the rear of the droplet, the active stress is larger there, creating a fluid flow that sustains the broken symmetry, by sweeping the cortex of actomyosin towards the rear so that −ψ is larger there than at the front (Fig. 1

(b) right).
This flow is governed by the hydrodynamic velocity v(r, t), which describes the average velocity of the cellular materials plus the solvent.The conserved dynamics of φ and ψ is then as follows: such that the total φ dr and ψ dr are constant in time.
The first term inside the brackets describes advection of φ and ψ by the fluid velocity v.The second term describes diffusion of φ and ψ along the negative gradient of the chemical potential δF/δφ and δF/δψ respectively.M φ,ψ > 0 are the mobilities for each field.
In the limit of low Reynolds number, which is appropriate for sub-cellular materials, the fluid flow v(r, t) is obtained from solving the Stokes equation: where σ = −pI + η[(∇v) + (∇v) T ] is the Cauchy stress tensor in a fluid of viscosity η, I is the identity matrix, and p is the isotropic pressure, which enforces the incompressibility condition ∇ • v = 0. Σ E in (3) is the equilibrium interfacial stress, which is derived from the free energy functional (1) [13]: Physically, Σ E is the elastic response to a deformation in the interface of φ-and ψ-droplet.Σ A in (3) is the active stress, which drives the system out of equilibrium.The form of Σ A is adapted from active model H [14,15]: where α ≥ 0 is a cortical contractile activity parameter, such that the equilibrium limit is recovered when α → 0. From ( 5), the activity is always localized at the interface of the φ-droplet.This differs from other models of active fluid droplets such as active nematics [6,11,16,17], where the active stress affects the bulk of the interior.The physical significance of Σ A is illustrated in Fig. 1(a).Consider an interface separating φ = +1 and φ = −1 regions.If αψ > 0 so that Σ A ∝ (∇φ)(∇φ), the active stress creates an extensile fluid flow in the tangential direction of the interface.On the other hand if αψ < 0, as will hold here, Σ A ∝ −(∇φ)(∇φ), and the active stress creates a contractile flow tangential to the interface, corresponding to actomyosin contractility α > 0 in the cell cortex.(Recall that the density of this cortex is ρ ∼ −ψ > 0 at the outer droplet interface.)Note that in the literature on active phase separation [14,15], contractile and extensile are defined with respect to the microswimmer orientation, normal to the interface.The opposite convention, chosen here, refers to the tangential cortex layer and is used in the cellular literature [3,4].
Our numerical system consists of a 2-dimensional square box with linear size L and periodic boundary conditions.We initialize a φ-droplet with radius R at the centre of the box, and similarly a ψ-droplet with a smaller radius 0.9R; see Fig. 1(b) left.We then solve Eqns.(2a-3) numerically using apseudo-spectral method, as detailed in the SI [18].
Mechanism of self-propulsion: Now we will illustrate the mechanism for self-propulsion.We fix the activity α to be finite and positive.First, let us consider what happens when both φ-and ψ-droplet are concentric as shown in Fig. 1(b) left.At the interface of the φ-droplet, the value of ψ is negative and the same everywhere along the interface (ψ −1).Thus, the active stress Σ A (5) is contractile and its magnitude is the same everywhere along the interface.This represents an isotropic distribution of active cortex material, and thus, by symmetry, we should not expect to see any motion.Now imagine that we give a small displacement (which can come either from thermal fluctuations or an initial perturbation -we consider the latter case here) to the ψ-droplet as shown in Fig. 1(b) right.Suppose that at the front, the interface of the φ-droplet coincides with that of the ψ-droplet.Thus the value of ψ is zero at the front of φ-droplet while ψ −1 at the back of φ-droplet.Therefore at the front, the active stress Σ A (5) is approximately zero, whereas at the back, the active stress is finite and contractile.This excess contractile stress at the back pulls the fluid from the front to the back along the interface of φ, as indicated by black arrows in Fig. 1(b) right.This translates into persistent motion.Fig. 2(a) shows the values of φ and ψ measured at the cross-section, which supports this mechanism.Physically, we have accumulation of actomyosin at the back of the cell which creates excess contractility at the rear cell cortex.This is consistent with experiments on cellular swimming motility in the absence of any boundary on which to crawl [4,19].Self-propulsion speed, U : Having described the mechanism, we now obtain an analytical expression for U using appropriate boundary conditions corresponding to our model, as given below: The above equations represent no-flux (6a), continuity of tangential slip velocity (6b) at the interface, and the fact that the discontinuity of the Cauchy stress, σ, at the interface is related to the surface tension γ(θ) of the droplet, via (6c) [24,25].We solve the Stokes equation both inside and outside the droplet using the above boundary conditions to obtain the the speed U 0 of the droplet in an infinite domain (see [18]): where ∆γ is the change in surface tension between the front and rear of the droplet.Our simulations are done for a given area fraction ϕ = πR 2 /L 2 using periodic boundary conditions inside a box of size L. To account for this, we need to sum the flow due to periodic images of the droplet, which lie on a two-dimensional square lattice [18].The final expression for the speed U is then Fig. 3(a) shows that U increases linearly with the activity parameter α and decreases linearly with the area fraction ϕ in accord with Eq. ( 8).However, the coefficients are different in theoretical calculation (8) and simulations (see Fig. 3).As described above, in presence of the active stress Σ A , ∆γ = α∆ψγ 0 /κ, and thus, we obtain the linear scaling with α using Eq.(7)(8).This implies that the speed increases linearly with contractile activity α and actomyosin concentration in the cortex ρ.In simulations, ψ-droplet is squashed inside the φ-droplet, and thus, there will also be a dipolar fluid flow [13], while we only consider a quadrupolar flow in our analysis [18].An absolutely convergent expression of the fluid flow [26,27] should then be used to obtain contributions from other multipoles, which we postpone to a future work.
Crowley-like instability of an array of active droplets: The model and simulations above can be easily extended to the case of multiple droplets.This is operationally done by generalizing the free energy functional of (1) for N active droplets as: Here the term proportional to β > 0 effectively leads to repulsion between the droplet φ i and φ j [28], and thus, precludes any overlap.The dynamics for each φ i , ψ i and v is the same as in (2a-3) and the equilibrium interfacial stress Σ E and active stress Σ A are now the sum of each contribution from φ i and ψ i given in ( 4) and (5).We use the above to study the instability in a linear array of three droplets as shown in Fig. 4. We see an active-Crowley-like instability when there is a well-defined central particle.Unlike in the case of the sedimentation [29], here the central particle lags.This can be understood from the fluid flow in Fig. 2(c), which has the effect of the central particle being pushed backwards by the neighboring particles.On the other hand if the particles are equally separated as in Fig. 4(b), there is no instability as there is no well-defined central particle, when taking into account the periodic boundary condition used.

Conclusion:
We have presented a minimal (and scalable) hydrodynamic model of active droplets without any liquid-crystalline order parameter [7,11].Nor do we require explicit chemical reactions in the form of source and sink terms [30][31][32].The self-propulsion of droplets in our model relies on the fact that activity sustains a finite surface tension gradient, ∆γ = 0, on the surface of the droplet.Thus, although motivated by a simple description of swimming cells, our model can also capture the self-propulsion due to Marangoni stresses on the surface of active emulsion droplets [21][22][23]33].Our theory might possibly be extended to address a passive liquid crystal inside an active scalar droplet, to mimic experimental systems of [21,33], or, in the chiral case, the helical trajectories seen in experiments of [34].
Our analytical results account for the linear scalings for the self-propulsion speed, as a function of activity parameter α and area fraction ϕ, seen in the simulations; however this belies the true complexity of the problem which exhibits linear scaling far beyond any perturbative regime (and indeed with a different coefficient).We also showed the feasibility of our model for the study of many droplets by studying an active Crowley-like instability in a linear array of active droplets.A more detailed study of an active droplet suspension using our theory and its comparison to the particulate theories [35] of active matter will be presented in future work.
RS is funded by a Royal Society-SERB Newton International Fellowship.MEC is funded by the Royal Society.Numerical work was performed on the Fawcett HPC system at the Centre for Mathematical Sciences.Work funded in part by the European Research Council under the Horizon 2020 Programme, ERC grant agreement number 740269.* rs2004@cam.ac.uk In this Section, we will derive the analytic solution for the fluid velocity in the limit of infinite boundary (L → ∞) and sharp φ-interface (ξ 0 → 0).We consider a two-dimensional active droplet swimming with velocity −U 0 x in the lab frame.In the co-moving frame, the droplet will be stationary, and the fluid velocity at far field r → ∞ is U 0 x.We can then define v(r) and ṽ(r) to be the fluid velocity inside and outside the droplet respectively; v(r) and ṽ(r) are obtained from solving two independent Stokes equation.We then match the two solutions at the interface of the droplet r = R via the boundary conditions in Eq.( 6) of the main text.We can assume the droplet to be centered at the origin.It is convenient to obtain the solution of the fluid flow using stream-function ϑ, defined as v = ∇ × (ϑẑ) [36].The components of the fluid velocity are given in terms of the stream-function in Cartesian (x, y) and plane polar (r, θ) coordinates as Now using the incompressibility condition ∇ • v = 0 and Stokes equation, it can be shown by standard arguments that the stream-function satisfies the biharmonic equation: subject to the boundary conditions for the corresponding velocity field given in the main text.Now let us denote v = ∇ × (ϑẑ) and ṽ = ∇ × ( θẑ) to be the fluid velocity inside and outside the droplet respectively.
Exterior flow v: First, we will solve the fluid velocity outside the droplet.The boundary condition for the fluid velocity at infinity r → ∞ is: This implies that the stream-function has the following form at r → ∞: Thus we can use separation of variables: Substituting ( 14) into (11), f (r) satisfies the following equation: Using trial solution f (r) ∼ r n , we can obtain the general solution to (15): From the boundary condition ( 13), we get C = D = 0 and B = U 0 .Now the no-flux boundary condition for the normal component of the velocity at the interface r = R reads: This implies A = −U 0 R 2 .Thus the exterior fluid velocity in polar coordinates and co-moving frame is: Interior flow ṽ: Again, defining the stream-function θ(r, θ) and using separation of variables, the solution to the biharmonic equation ( 11) is Now since the fluid velocity ṽ = ∇ × ( θẑ) has to remain finite at r = 0, we then require Ã = C = 0. We then impose no-flux boundary condition: and continuity of the tangential velocity at the interface: to finally find: B = −U 0 and D = U 0 /R 2 .Thus the interior fluid velocity in polar coordinates is: This fluid flow, derived from above expressions for the interior and exterior region, has been plotted in Fig. (5).This can be compared with direct numerical simulations in the main text.It is worthwhile to note that there is no solution for Stokes flow around a disc, the so-called Stokes paradox [36], if a no-slip boundary condition at r = R is used.The solution described above relies on the fact that there is a free-slip boundary condition on the surface of the droplet.
Self-propulsion speed: The hydrodynamic stress is discontinuous at the interface r = R due to the surface tension γ(θ): The stress in polar coordinates is: First, let us look at the r-component of ( 25): which is just the Laplace pressure difference between the interior and exterior of the droplet.Now we look at the θ-component of ( 25): where η and η is the fluid viscosity inside and outside the droplet respectively.Now substituting ( 19) and ( 24) into (31), we get: Thus, the velocity of the droplet is non-zero only if we have a surface tension gradient along the tangential direction.Now we can average (32) over the angle θ ∈ [0, π] to obtain: where ∆γ is the surface tension difference between the back and the front.Note that θ = π corresponds to the front the droplet and θ = 0 corresponds to the back of the droplet.For active model H [14], the effective surface tension is γ = γ 0 (1 − αψ/κ), and thus we get: since ψ = 0 in the front and ψ = −1 in the back (see Fig. 1 in the main text).Thus, we obtain linear scaling with activity α, consistent with the numerical result in the main text.

II. ACTIVE DROPLETS ON A LATTICE
To compute this finite size scaling theoretically, we consider an infinite array of active droplets, whose centre of mass are located on a square lattice.The renormalized fluid velocity in the region around the central droplet, located at the origin, is the sum of all the fluid velocities generated by each droplet on the lattice.The numerical simulations presented in the main text assume periodic boundary condition on each side of the box: (x, ±L/2), (±L/2, y).This is equivalent to having an infinite number of active droplets located on a square lattice (mL, nL), where m, n ∈ Z.In this Section, we will derive the approximate fluid velocity generated by these droplets.
First, let us consider the dilute limit R/L → 0 (Section I).Let us consider an active droplet swimming with velocity U 0 x and located at the origin.The fluid flow generated by this single droplet in the lab frame is (see Fig. 5(b) and Eqns.(18)(19)): In particular, the fluid flow at the leading edge of a single droplet (x, y) = (R, 0) is equal to the velocity of the droplet itself: Now, the fluid flow generated by infinite droplets in a lattice is given by (assuming dilute flow solution for each droplet): The first term is the fluid flow created by the active droplet at the origin.In particular, the fluid flow at the edge of the centre droplet is given by: where G(R/L) is the infinite sum contained in (39).It can be shown that the infinite sum converges and G(R/L) as a function of R/L is shown in the plot of Fig. 6.From the plot, G(R/L) constant for small R/L and the constant value is −π.Therefore, the fluid velocity at (x, y) = (R, 0) is Now we can compare (41) to (37) to deduce that the velocity of the central droplet at the origin is renormalized Thus we find linear scaling with decreasing ϕ and increasing α or ∆γ, consistent with the numerical result in the main text, although the pre-factors are not exactly matched.The reason is because the neighboring droplets also renormalize the fluid velocity at the interface of the central droplet r = R. Thus the stress boundary condition (25), no-flux condition, and continuity of tangential velocity at the interface has to be reevaluated.Secondly in simulations, inside the φ-droplet, there is also a ψdroplet, which gets squashed.This contributes to a dipolar fluid flow [15], while we only consider the quadrupolar flow in our analysis.A more detailed study, accounting for all these effects and the role of periodic images in evaluation the fluid flow at all orders using a absolutely convergent expression of the fluid flow [26,27], will be presented in a future work.

III. SIMULATION DETAILS
The fluid flow satisfies the Stokes equation where, is the force density on the fluid.The solution is obtained using Fourier transforms as we now describe.Defining the Fourier transform of a function φ(r) as φ k = φ(r)e −ik•r dr (46) φ(r) = 1 (2π) 3  φ k e ik•r dk, The simulations reported are in two space dimensions, while the analytic predictions would be maintained in higher dimensions.The parameter β is fixed to be small and positive β = 0.02 so that the ψ-droplet is confined inside the larger φ-droplet.We keep the ratio of the radius of the two droplets fixed at 0.9, such that φ droplet is bigger than the ψ droplet.We fix the mobilities to be unity: M φ = M ψ = 1 and the viscosity to be η = 0.1.The values of L, R, and α vary.We also define the area fraction to be ϕ = πR 2 /L 2 .The time discretization is fixed to be ∆t = 0.003 and spatial discretization is fixed to be ∆x = 1.

Figure 1 .
Figure 1.(a) Flow driven by an active stress Σ A on the interface of the droplet phase field φ.A positive (negative) Σ A implies an extensile (contractile) fluid flow tangential to the φ-interface.(b) A cell is represented by φ-droplet, while concentration of actomyosin in the cell cortex given by φ−ψ (red ring).Under a small random initial perturbation, the ψ > 0 domain is displaced to the front (right), creating excess contractile stress at the back.This leads to an active flow which squashes the ψ-droplet in the front, causing sustained self-propulsion along a direction selected by spontaneous symmetry breaking.

Figure 2 .
Figure 2. (a) Inset: a snapshot of an active droplet steadily moving in the x-direction.Main: corresponding plots of φ and ψ as a function of x at the cross-section passing through the centre of the droplet.(b,c) show the streamlines of the fluid flow of the same active droplet in (b) comoving frame and (c) lab frame.The streamlines are superimposed on the pseudo-color plot of the φ-field.(Parameter used: α = 0.11, ϕ = πR 2 /L 2 = 0.126 and L = 100.)

Fig. 2 (
b,c) show the steady state fluid velocity from the full hydrodynamic simulations in the co-moving (b)

Figure 3 .
Figure 3. (a) The self-propulsion speed of the droplet U increases linearly with the activity parameter α at a fixed area fraction ϕ = πR 2 /L 2 = 0.049.(b) The self-propulsion speed of the droplet U decreases linearly with the area fraction ϕ for a fixed α = 0.2.From the best linear fit of the simulation data we obtain U = α (0.128 − 0.422 ϕ).System size L = 100.

Figure 4 .
Figure 4. (a) Active Crowley-like instability: a linear array of active droplets is unstable.This instability is opposite to that in an array of sedimenting particles.System size is L = 120.(b) There is no instability if there is no well-defined central particle.System size is L = 90.

Figure 5 .
Figure 5. Streamlines of the theoretical fluid flow around an active droplet moving to the right in the co-moving (a) and lab-frame (b) in an infinite domain.The red circle denotes the interface of the droplet.(Note that the calculation in Section I is done with the droplet moving to the left, thus the fluid velocity above is flipped around the y-axis.)

Figure 6 .
Figure 6.Shows the plot of the infinite sum in (39) as a function of R/L,