Rayleigh–Plateau instability of anisotropic interfaces. Part 1. An analytical and numerical study of fluid interfaces

Abstract

Numerous experiments and theoretical calculations have shown that cylindrical vesicles can undergo a pearling instability similar to the Rayleigh–Plateau instability of a liquid jet when they are subjected to external tension. In a living cell, a Rayleigh–Plateau-like instability could be triggered by internal tension generated in the cell cortex. This mechanism has been suggested to play an essential role in biological processes such as cell morphogenesis. In contrast to the simple, passive and isotropic membrane of vesicles, the cortical tensions generated by biological cells are often strongly anisotropic. Here, we theoretically investigate how this anisotropy affects the Rayleigh–Plateau instability mechanism. We do so in the limit of both low and high Reynolds numbers and accordingly cover cell behaviour under anisotropic cortical tension as well as fast liquid jets with anisotropic surface tension. Combining analytical linear stability analysis with numerical simulations we report a strong influence of the anisotropy on the dominant wavelength of the instability: increasing azimuthal with respect to axial tension leads to destabilisation and to a shorter break-up wavelength. In addition, compared to the classical isotropic Rayleigh–Plateau situation, the range of unstable modes grows or shrinks when the azimuthal tension is higher or lower than the axial tension, respectively. We explore nonlinear effects like an altered break-up time and formation of satellite droplets under anisotropic tension. In Part 2 (Bächer et al. J. Fluid Mech., vol. xxx, 2021, Ax) of this series we continue our analysis by analytically investigating the influence of bending and shear elasticity, usually present in vesicles and cells, on this anisotropic Rayleigh–Plateau instability.

1. Introduction

The break-up of liquid jets into droplets, triggered by surface tension, was already investigated intensively by Plateau in the second half of the 19th century (Plateau 1873). Based on the concept of the fastest growing perturbation, Rayleigh derived a relation between the radius of the liquid jet and the dominant wavelength which determines the size of the droplets for an ideal fluid in the absence of an outer medium (Rayleigh 1878). Later, he extended the theoretical description of this Rayleigh–Plateau instability to viscous jets in the inertialess Stokes limit (Rayleigh 1892). Again in the Stokes limit, the presence of an outer medium with arbitrary viscosity ratio between inner and outer fluid has been investigated by Tomotika (1935). A simplified version for the important case of equal viscosities has been presented by Stone & Brenner (1996). The Rayleigh–Plateau instability is a prime example of the beauty of fluid mechanics and possesses great relevance in various applications. We refer to the review article by Eggers & Villermaux (2008) for further details.

However, pearling and break-up due to the Rayleigh–Plateau mechanism are not restricted to liquid jets. In 1994, Bar-Ziv & Moses (1994) reported a pearling instability for a tubular vesicle. A vesicle consists of a lipid bilayer membrane confining an interior fluid and is often considered as a model system for a biological cell (Seifert 1997). Under local application of laser tweezers the vesicle formed pearls (Bar-Ziv & Moses 1994). Using a hydrodynamic theory Nelson, Powers & Seifert (1995) and Goldstein et al. (1996) explained the pearling of the vesicle by a laser induced tension, which in turn triggers a Rayleigh–Plateau-like instability. Pearl formation starts at the site of application of the laser and the instability then propagates along the cylindrical vesicle (Goldstein et al. 1996; Bar-Ziv, Tlusty & Moses 1997). Later, several experimental studies demonstrated different ways to induce the tension which is required for the pearling instability (Powers 2010). Pulling on membrane tethers with optically trapped particles (Bar-Ziv, Moses & Nelson 1998; Powers, Huber & Goldstein 2002), protein mediated anchoring of membrane tethers to a substrate (Bar-Ziv et al. 1999), applying a magnetic field (Ménager et al. 2002), electric field (Sinha, Gadkari & Thaokar 2013) or osmotic pressure gradient (Yanagisawa, Imai & Taniguchi 2008; Sanborn et al. 2013) can all lead to pearling. Furthermore, Kantsler, Segre & Steinberg (2008) reported the transition of a finite-size, tubular vesicle to a pearling state due to stretching in an extensional flow and noted that the transition is reversible when the external flow stops. In shear flow the instability has also been observed (Pal & Khakhar 2019). Boedec, Jaeger & Leonetti (2014) derived theoretically the growth rate for the instability of a cylindrical vesicle under tension. They treat the fluid surrounding the vesicle in the limit of the Stokes equation and allow for variations of the tension along the vesicle. By means of boundary integral simulations Narsimhan, Spann & Shaqfeh (2015) showed that the initial shape of a closed vesicle in extensional flow influences the number of fragments after pearling.

In contrast to passive vesicles, where the tension triggering the Rayleigh–Plateau instability has to be imposed from the outside, living biological cells are able to internally create active stresses in their cytoskeletal network (Kruse et al. 2005; Marchetti et al. 2013; Prost, Jülicher & Joanny 2015; Salbreux & Jülicher 2017; Jülicher, Grill & Salbreux 2018). Such cytoskeletal networks can build a thin layer that underlies the plasma membrane, named the cell cortex (Alberts et al. 2007; Köster & Mayor 2016; Chugh & Paluch 2018), in which the action of motor proteins leads to active tension at the cell's interface (Chugh et al. 2017). A positive constant active tension caused by a homogeneous (Pleines et al. 2013) actomyosin distribution in the cortex describes the internal tendency of the cytoskeleton to contract (Needleman & Dogic 2017). Alternatively, proteins which anchor at the plasma membrane can trigger a pearling instability (Tsafrir et al. 2001) by bending the membrane and thus inducing a non-zero curvature (Campelo & Hernández-Machado 2007; Jelerčič & Gov 2015). For a viscous active surface Mietke et al. (2019a) and Mietke, Jülicher & Sbalzarini (2019b) report a Rayleigh–Plateau instability with mechano-chemical regulation. For a biological tissue composed of multiple cells Hannezo, Prost & Joanny (2012) provide an energy argument based on an effective surface tension. Berthoumieux et al. (2014) considered the Green's function for an elastic cell membrane subjected to active tension, which again leads to the prediction of a Rayleigh–Plateau instability. Bächer & Gekle (2019) confirmed the instability threshold predicted by Berthoumieux et al. (2014) and presented the shape of a membrane undergoing Rayleigh–Plateau instability in three-dimensional simulations of active membranes. For soft materials the Rayleigh–Plateau instability is an important mechanism in beaded object formation (Mora et al. 2010) and the production of synthetic vesicles (Anna 2016; Pal & Khakhar 2019). Especially in the biological context, Rayleigh–Plateau-like instabilities have been proposed to play an important role in microtubuli-driven cell deformation (Emsellem, Cardoso & Tabeling 1998), as a driving mechanism behind mitochondrial fission (Gonzalez-Rodriguez et al. 2015) as well as for pathological shapes of blood vessels during vasoconstriction (Alstrøm et al. 1999).

All the above mentioned studies on the Rayleigh–Plateau instability in different contexts have in common that they consider an isotropic tension. However, in reality, cytoskeletal systems often exhibit strong anisotropy (Reymann et al. 2012; Murrell et al. 2015; Blackwell et al. 2016; Zhang et al. 2018), e.g. due to the formation of stress fibres (Tojkander, Gateva & Lappalainen 2012). Accordingly, the tension in the cell cortex can be anisotropic (Rauzi et al. 2008; Mayer et al. 2010; Behrndt et al. 2012; Callan-Jones et al. 2016), which is important for many biological phenomena such as cell-shape regulation (Callan-Jones et al. 2016), cell polarisation (Mayer et al. 2010), ingression formation (Reymann et al. 2016), the formation of a furrow during cell division (White & Borisy 1983; Salbreux, Prost & Joanny 2009) and the production of blood platelets (Bächer, Bender & Gekle 2020). For a solid rod in the absence of any kind of fluid, Gurski & McFadden (2003) proposed an instability mechanism based on the bulk anisotropy of the underlying crystal lattice for the growing of nanowires. How anisotropic surface tension affects the Rayleigh–Plateau instability of vesicles, cells or even liquid jets, however, remains an open question.

In this work, we analytically extend the framework of the Rayleigh–Plateau instability to include anisotropic interfacial tension for low (Stokes fluid) and high (ideal fluid) Reynolds numbers. In both situations, we derive the dispersion relation depending on the tension anisotropy and report a striking influence on the dominant wavelength and maximum growth rate of the instability. Compared to the classical Rayleigh–Plateau criterion for isotropic surface tension, we observe a decrease in wavelength for dominating azimuthal tension and an increase for dominating axial tension. The analytical predictions agree very well with numerical simulations using a boundary integral method (BIM) and a lattice-Boltzmann/immersed boundary method (LBM/IBM). From these simulations we also compute the nonlinear correction to the linear break-up time. Including interface viscosity in the stability analysis for the Stokes regime also influences the dominant wavelength and growth rate of the instability albeit less pronounced than the tension anisotropy. Finally, we use a long-wavelength expansion to investigate the formation of satellite droplets (Ashgriz & Mashayek 1995; Martínez-Calvo et al. 2020) under anisotropic interfacial tension. In Part 2 (Bächer, Graessel & Gekle 2021) we consider the anisotropic Rayleigh–Plateau instability of vesicles or capsules endowed with bending and shear elasticity.

We start by introducing our theoretical model for an anisotropic interface, the coupling to the surrounding fluid as well as the numerical methods used in the simulations in § 2. We then present the dispersion relation for the Rayleigh–Plateau instability of an anisotropic interface obtained by analytical linear stability analysis in § 3 first for a Stokes fluid in § 3.1 and then for an ideal fluid in § 3.2. In § 4.1 we discuss the effect of anisotropic interfacial tension on the dominant wavelength of the instability by comparing analytical and simulation results and show a transition between Stokes fluid and ideal fluid in § 4.2. In § 4.3 we discuss the dominant growth rate in comparison to numerical results. Nonlinear corrections of the linear break-up time are investigated in § 4.4. In § 5 we discuss the combination of anisotropic tension and interface viscosity and finally present the formation of satellite droplets under the influence of tension anisotropy for an ideal fluid jet without ambient fluid in § 6. We conclude in § 7.

2. Description of an anisotropic interface

2.1. Problem illustration

We consider a general complex interface, as sketched in figure 1, which is surrounded by a fluid on both sides. This can either represent the interface of a liquid jet in the co-moving frame or the membrane of a vesicle or biological cell. As usual (Eggers & Villermaux 2008; Boedec et al. 2014), we assume that the interface is infinitely long. In the analytical stability analysis we consider an axisymmetric interface, which is parametrised by the axial position $z$ and the local radius $R(z,t)$. Initially, the interface is cylindrical with radius $R(z,0)=R_0$. At arbitrary time $t$ the interface is subjected to a perturbation $\delta R(z,t)$, such that the local radius is given by $R(z,t) = R_0 + \delta R(z,t)$.

Figure 1. Illustration of the set-up. We consider a complex interface which can be either a liquid jet of Newtonian fluid in the limit of vanishing viscosity $\eta$ or the membrane of a vesicle or cell immersed in a fluid in the limit of the Stokes equation, i.e. density $\rho =0$. The fluid jet is immersed in an ambient fluid with $\eta ^0, \rho ^0$. The cylindrical interface of initial radius $R_0$ (dashed line) is subjected to a periodic perturbation with amplitude $\epsilon$ (solid blue line). The interface is parametrised by the position along the cylinder axis $z$ and the radius $R(z,t)$. We consider the interfacial tension in the axial direction $\gamma^z$ (orange) different from that in the azimuthal direction $\gamma^\phi$ (green), both of which contribute to the membrane force acting onto the fluid with different curvature components (grey circles).

In order to perform a linear stability analysis of the complex interface in the presence of anisotropic interfacial tension, we apply a periodic perturbation to its shape (Drazin & Reid 2004). The perturbation of the interface is illustrated in figure 1: it modulates the radius in $z$-direction along the cylinder axis with amplitude $\epsilon (t)=\epsilon _0\,{\textrm {e}}^{\omega t}$, a wavelength $\lambda$ and a corresponding wavenumber $k={2{\rm \pi} }/{\lambda }$ of the wave vector pointing along the cylinder axis. The perturbation with initial amplitude $\epsilon _0$ grows in time with growth rate $\omega$. Accordingly, the interface of the jet can be described by its radius as

(2.1)\begin{equation} R(z,t)=R_0+\epsilon_0 \exp({\omega t +{\textrm{i}} kz}). \end{equation}

Throughout this work, we consider an anisotropic interfacial tension, i.e. the value of the axial tension differs from the value of the azimuthal tension. This anisotropic tension accounts for two fundamentally different situations. First, in a liquid jet anisotropy can arise, e.g. from covering the interface with passive anisotropic surfactant molecules, thus extending the classical concept of liquid–liquid surface tension to an anisotropic situation. Second, in biological cells or tissue, an active biological machinery, cytoskeletal filaments with motor proteins underlining the plasma membrane, can produce anisotropic tensions at the interface as described in more detail in the Introduction. Due to their usually contractile nature, these active tensions enter the physical equations in the same way as the classical surface tension, despite their fundamentally different origin. In the following, we therefore use the same symbol and refer to both scenarios with the general term interfacial tension.

2.2. Interface coupled to a surrounding fluid

Interfacial tension leads to internal forces being transmitted from the interface to the fluid (Green & Zerna 1954; Barthès-Biesel 2016; Salbreux & Jülicher 2017). In contrast to the classical isotropic Rayleigh–Plateau scenario, we assign anisotropic tension to the interface, i.e. we distinguish between the azimuthal $\gamma ^{\phi }$ and the axial interfacial tension $\gamma ^{z}$. As sketched in figure 1 the periodic perturbation along the axis changes the curvature of the interface both in the azimuthal and in the axial direction. In azimuthal direction the curvature ${1}/{R_\phi }$ is the inverse of the local radius of the interface $R_\phi =R(z,t)$, where we follow the convention that the curvature of a cylinder is positive. Accordingly, the curvature along the axis is given by the negative second derivative of the radius, ${1}/{R_z} = - R^{\prime \prime }$ such that the curvature is negative at a neck and positive at a bulge (compare figure 1). The derivative $R^{\prime \prime }$ follows directly from (2.1).

The anisotropic tension does not depend on the position along the interface, therefore its derivative vanishes and for a liquid–liquid interface no internal forces tangential to the interface arise (Green & Zerna 1954; Salbreux & Jülicher 2017). The internal force normal to the interface is given by the interfacial tension components weighted by the corresponding principal curvature. Balance of forces requires that this normal force is in equilibrium with the difference in tractions exerted by the fluids on either side of the interface. Thus, the normal traction jump across the interface ${\rm \Delta} f^n$ reads

(2.2)\begin{equation} \frac{ \gamma^\phi}{R_\phi} + \frac{ \gamma^z}{R_z} = {\rm \Delta} f^n. \end{equation}

The traction jump is given by the projection of the three-dimensional viscous stress tensor of the outer and inner fluid onto the interface normal vector (Chandrasekharaiah & Debnath 1994). For an incompressible interface or for negligible viscous effects, i.e. for an ideal fluid, the traction jump is determined by the pressure $p$ of the fluid. With the normal vector pointing outwards from the interface and considering the outer and inner fluid as incompressible, the traction jump in normal direction is thus given by

(2.3)\begin{equation} {\rm \Delta} f^n = -p^{{out}} + p^{{in}} = p (r=R), \end{equation}

with pressures $p^{{out}}$ and $p^{{in}}$ of the outer and inner fluid, respectively, and $p (r=R)$ denoting the pressure difference at the interface. Together (2.2) and (2.3) lead to

(2.4)\begin{equation} \frac{\gamma^\phi}{R_\phi} + \frac{\gamma^z}{R_z}= p (r=R), \end{equation}

for an incompressible interface or an ideal fluid interface. The relation in (2.4) reduces to the classical Young–Laplace equation ${\gamma }/{R} = p$ in the limit of isotropic surface tension $\gamma =\gamma ^\phi =\gamma ^z$ and vanishing curvature along $z$, i.e. $R_z \rightarrow \infty$. The anisotropic interfacial tension thus leads to a pressure disturbance of the inner fluid, where the two contributions of the interfacial tension are weighted with their respective radii of curvature.

2.3. Fluid dynamics

The motion of the fluid inside and outside the jet is in general described by the Navier–Stokes equation

(2.5)\begin{equation} \frac{\partial \boldsymbol{v}}{\partial t}+ \left(\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{v} = -\frac{1}{\rho} \boldsymbol{\nabla} p + \nu {\rm \Delta} \boldsymbol{v}, \end{equation}

with the velocity field $\boldsymbol {v}$, fluid density $\rho$, kinematic viscosity $\nu = {\eta }/{\rho }$ and shear viscosity $\eta$. Density and viscosity of the outer fluid are denoted by $\rho ^o$ and $\eta ^o$, respectively. The fact that the liquid is incompressible, which is true for both the liquid of a jet and the liquid encapsulated in a vesicle, leads furthermore to the continuity equation for the incompressible liquid

(2.6)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0. \end{equation}

The Navier–Stokes and continuity equations together govern the motion of the fluid.

Besides the well-known Reynolds number $\textit {Re} = {\rho R_{0} V_{0}}/{\eta }$ with a typical velocity $V_{0}$ and the unperturbed interface radius $R_0$ as the typical length, we use the Ohnesorge number (Eggers & Villermaux 2008), which relates characteristic scales of the interface and the surrounding fluid

(2.7)\begin{equation} \textit{Oh} = \frac{\eta}{\sqrt{\rho R_0 \gamma}}. \end{equation}

In our anisotropic scenario with $\gamma ^{z} \neq \gamma ^{\phi }$ we define two distinct Ohnesorge numbers $\textit {Oh}_z$ and $\textit {Oh}_{\phi }$ for the respective interfacial tensions.

In the limit of small velocities or large viscosity, i.e. the Stokes regime, the Reynolds number approaches zero while the Ohnesorge number becomes large. In this regime, the Navier–Stokes equation can be replaced by the linear Stokes equation

(2.8)\begin{equation} 0 = -\boldsymbol{\nabla} p + \eta {\rm \Delta} \boldsymbol{v}. \end{equation}

In the limit of large velocities or small viscosity, i.e. for an ideal fluid, the Reynolds number is larger than one and the Ohnesorge number becomes small. Here, the Euler equation applies

(2.9)\begin{equation} \frac{\partial \boldsymbol{v}}{\partial t}+ \left(\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{v} = -\frac{1}{\rho} \boldsymbol{\nabla}p. \end{equation}

2.4. Numerical simulations

We aim for a comparison of our main analytical results, the dominant wavelength of the instability and its growth rate presented in §§ 4.1 and 4.3, respectively, with numerical simulations solving the coupled fluid and interface dynamics. The simulations further provide us a glimpse on the nonlinear aspects of the instability dynamics. For a Stokes fluid and an ideal fluid with an outer fluid with the same properties, we perform three-dimensional boundary integral method and lattice-Boltzmann/immersed boundary method simulations, respectively.

2.4.1. Three-dimensional numerical investigation of the instability

We consider a fluid column, the liquid jet, immersed in an ambient fluid. We use fully three-dimensional simulations, thus testing also for non-axisymmetric instabilities triggered by anisotropic interfacial tension (which we did not observe).

The interface encapsulates a Newtonian fluid and is surrounded by another Newtonian fluid of the same density $\rho ^0 = \rho$ or viscosity $\eta ^0 = \eta$. For the fluid, periodic boundary conditions are chosen in each of the three spatial directions together with a kinematic boundary condition at the interface. Fluid dynamics is either solved by the BIM or the LBM/IBM, as detailed below.

Following the set-up sketched in figure 1, we consider an initially cylindrical interface which is modelled as a thin shell and which we discretise by nodes connected to triangles. Anisotropic tension of the interface is realised using the recently developed and validated computational method for active membranes in flows (Bächer & Gekle 2019). This approach can treat both the anisotropic surface tension of a liquid jet in the co-moving frame and the anisotropic active tension of a biological cell cortex on the same footing. To investigate the instability dynamics, we initially apply a small periodic perturbation to the cylindrical interface, as shown in figure 1. From this initial configuration the temporal evolution of the interface and the suspending fluid is solved including a dynamical two way coupling of interface and fluid. The method to determine the dominant wavelength and growth rate is described in appendix A.

2.4.2. Boundary integral method

As the simulation method at zero Reynolds number we use the BIM to solve the fluid dynamics (Pozrikidis 2001; Zhao et al. 2010). The BIM solves the Stokes equation in the presence of discretised boundaries based on hydrodynamic Green's functions. It directly solves for the fluid velocity at the nodes of the discretised interface for given interface shape and interfacial force density. As a consequence of using the Stokes equation, in BIM simulations inertial effects are excluded. Neglecting inertial effects corresponds to $\textit {Re} = 0$ and an Ohnesorge number $\textit {Oh} \rightarrow \infty$. Membrane forces due to interfacial tension are calculated as detailed in Bächer & Gekle (2019). As the interfacial tension in the simulation we use $\gamma ^\phi \approx 10\ \textrm {pN}\ \mathrm {\mu }\textrm {m}^{-1}$, a typical tension expected for blood cells (Dmitrieff et al. 2017). For details on the implementation of the BIM we refer to Guckenberger & Gekle (2018).

As the box size, we consider the length of the cylindrical tube, which is typically 80 times the tube radius, along the axis and ten times the tube radius in lateral directions. A typical discretisation of the interface consists of approximately 17 000 nodes and 32 500 triangles. We do not use periodic boundary conditions for the membrane due to technical issues of the implementation used for BIM simulations. We rather place the outer rings of nodes exactly at the beginning and end of the box, respectively, and fix the nodes by elastic springs. Due to an insufficient number of neighbouring nodes, at the boundary nodes the force from the interfacial tension is not calculated. We note that using this set-up, the fluid encapsulated by the membrane remains inside. For the initial small deformation of the interface we choose the amplitude $\epsilon _0=0.02$. The fluid viscosity is chosen as $\eta \approx 1.2 \times 10^{-3}$ Pa s. We simulate for approximately 100 time steps with adaptive step size and a total simulation time of approximately 20 ms.

2.4.3. Lattice-Boltzmann/immersed boundary method

As a method to solve fluid dynamics at large/finite Reynolds number, we use the LBM (Aidun & Clausen 2010; Krüger et al. 2016) together with the IBM (Peskin 2002; Mittal & Iaccarino 2005; Bächer, Schrack & Gekle 2017; Mountrakis, Lorenz & Hoekstra 2017; Bächer et al. 2018). Our LBM/IBM is implemented in the software package ESPResSo (Limbach et al. 2006; Roehm & Arnold 2012; Arnold et al. 2013; Weik et al. 2019) and has been extensively validated (Gekle 2016; Guckenberger et al. 2016; Bächer et al. 2018; Bächer & Gekle 2019).

The LBM solves the fluid dynamics on the basis of the mesoscopic Boltzmann equation and accounts for the fluid dynamics according to the full Navier–Stokes equation. The fluid thus has a finite density, a finite viscosity and, therefore, a finite Ohnesorge number. The fluid is discretised by an Eulerian grid and populations representing the distribution functions for the different velocities are assigned to each fluid node. We here use the D3Q19 velocity set and a typical fluid mesh with dimensions of approximately $650\times 40\times 40$ with some simulation lattices extending up to $800\times 40\times 40$. A typical simulation runs for 500 000 steps. In the limit of an ideal fluid we choose Ohnesorge numbers in the range of $10^{-3}\text {--}10^{-4}$. Initially, the fluid has zero velocity.

The discretised interface is coupled to the background fluid using the IBM. A typical interface contains 18 240 nodes and 36 480 triangles, has a radius of 6 LBM grid cells and is periodic along the axial direction with an initial perturbation amplitude $\epsilon _0 = 0.02$ in simulations to determine the dominant wavelength. An additional refined simulation set-up is used for determination of the dominant growth rate, where we simulate one period of the dominant mode with initial perturbation $\epsilon _0 = 0.002$ and increased resolution with a radius of 13 LBM grid cells. Here, a typical fluid lattice consists of $180\times 70\times 70$ nodes and a membrane mesh of 15 416 nodes and 30 832 triangles. We again note that axisymmetry is not imposed and that the simulations are fully three-dimensional. The average distance between two interface nodes is approximately the length of one LBM grid cell. An interface node moves with the local fluid velocity which is interpolated at the node position from the surrounding fluid nodes by an eight-point stencil. The force stemming from the interfacial tension and acting from the membrane onto the fluid at the site of each interface node is transmitted to the fluid by the same eight-point stencil interpolation scheme. Thus, the IBM provides a dynamic two way coupling of membrane and fluid.

3. Dispersion relation for anisotropic interfacial tension

3.1. Anisotropic Rayleigh–Plateau instability for a Stokes fluid

Biological cells as well as their synthetic counterpart (vesicles) are typically a few tens micrometres in size. Therefore, we consider the limit of small Reynolds numbers, i.e. $\textit {Re} \ll 1$, where the Navier–Stokes equation (2.5), reduces to the linear Stokes equation (2.8) which together with the continuity equation (2.6) describes the fluid behaviour. In the following, we consider identical fluid viscosity inside and outside the vesicle, i.e. $\eta ^o = \eta$. Our aim is to obtain the dispersion relation in the case of anisotropic interfacial tension, which gives the growth rate depending on the wavenumber of the perturbation (2.1). As detailed in appendix B.1 we perform a linear stability analysis and obtain the dispersion relation for a Stokes fluid

(3.1)\begin{align} \omega (k) &= \frac{\gamma^\phi}{R_0 \eta} \left( 1 - \frac{\gamma^z}{\gamma^\phi}(R_0k)^2 \right) \left[\vphantom{\frac{kR_0}{2}} \textrm{I}_1(kR_0)\textrm{K}_1(kR_0)\right. \nonumber\\ &\quad \left. + \frac{kR_0}{2} \left( \textrm{I}_1(kR_0)\textrm{K}_0(kR_0) - \textrm{I}_0(kR_0) \textrm{K}_1(kR_0) \right) \right], \end{align}

with $\textrm {I}_\nu (x), \textrm {K}_\nu (x)$ being the modified Bessel functions of the first and second kind, respectively, and of order $\nu$.

Positive values of $\omega$ correspond to growing perturbations (2.1), whereas perturbation modes with negative growth rate are dampened. Because of the positive prefactor $\omega _0^{S}={\gamma ^{\phi }}/({R_0\eta })$, which is the inverse of the viscocapillary time based on the azimuthal tension $\gamma ^{\phi }$, and the positive modified Bessel functions for positive $kR_0$, the tension anisotropy determines the range of growing, i.e. unstable, modes. We obtain from (3.1) the range of growing modes for values of $kR_0$ between $-\sqrt {{\gamma ^{\phi }}/{\gamma ^{z}}}$ and $\sqrt {{\gamma ^{\phi }}/{\gamma ^{z}}}$. This range depends on the square root of the anisotropy of the interfacial tension. We recover for isotropic tension $\gamma ^z=\gamma ^\phi =\gamma$ the range of growing wavelengths between $-\sqrt {{\gamma ^{\phi }}/{\gamma ^{z}}} = -1$ and $\sqrt {{\gamma ^{\phi }}/{\gamma ^{z}}} = 1$ as found by Plateau (1873). So in the case of the classical Rayleigh–Plateau instability, the growing wavelengths do not depend on the interfacial tension $\gamma$ but only on the undisturbed radius $R_0$ of the jet (Drazin & Reid 2004). Here, in addition the anisotropy of interfacial tension enters as a factor.

We show in figure 2(a) the dispersion relation of the classical Rayleigh–Plateau instability, i.e. for isotropic interfacial tension. The growth rate $\omega$ is plotted only against positive $kR_0$ due to symmetry. We further distinguish the individual contributions from $\gamma ^\phi$, the first term in (3.1), and from $\gamma ^z$, the second term in (3.1), and plot them in orange and green, respectively, together with the total dispersion relation in blue. It is the interplay of the two contributions of the interfacial tension $\gamma ^z$ and $\gamma ^\phi$ that determines the dispersion relation. The azimuthal tension $\gamma ^{\phi }$, i.e. the first term in (3.1), is positive and thus the system would be unstable against any perturbation with arbitrary wavenumber. However, this is not the case, because this term is balanced by the damping contribution from $\gamma ^z$. Both contributions together determine a finite maximum of the dispersion relation, which corresponds to the dominant mode that grows fastest.

Figure 2. Dispersion relation in the Stokes regime for $\eta =\eta ^o$. Curves are shown for (a) isotropic interfacial tension ${\gamma ^{z}}/{\gamma ^{\phi }} = 1.0$ and for anisotropic interfacial tension with (b) ${\gamma ^{z}}/{\gamma ^{\phi }} = 0.5$ and (c) ${\gamma ^{z}}/{\gamma ^{\phi }} = 2.0$. We distinguish the contributions from $\gamma ^\phi$ (green) and $\gamma ^z$ (orange). An anisotropic tension strongly alters the range of growing modes and shifts the maximum towards larger $kR_0$ in (b) or smaller $kR_0$ in (c). (d) Dispersion relation for vanishing axial interfacial tension, i.e. $\gamma ^z=0$. The $\gamma ^\phi$ contribution (green) has its maximum at $kR_0=1.59$ in each of the panels, because $\gamma ^\phi$ is kept constant. Thus, although all modes are unstable in (d), in principle, there still exists a well-defined finite dominant wavelength for a Stokes fluid due to fluid stresses.

We now consider an anisotropic interface were the contributions $\gamma ^\phi$ and $\gamma ^z$ are no longer identical. Their changing ratio leads to a different weighting of the contributions to the dispersion relation (3.1). If the destabilising contribution from $\gamma ^\phi$ rises compared to $\gamma ^z$ as shown in figure 2(b), the range of growing wavelengths increases. On the other hand, if $\gamma ^\phi$ decreases relative to $\gamma ^z$, the range becomes smaller (see figure 2c). Because the destabilising $\gamma ^{\phi }$ contribution reaches a maximum at $kR_0=1.59$ and tends to zero for $kR_0\to \infty$, independent of the anisotropy ratio, also for vanishing $\gamma ^z \rightarrow 0$ a well-defined mode at finite $kR_0$ has the largest growth rate. This is shown by figure 2(d) in the case of $\gamma ^z=0$ with an extended $kR_0$-range on the horizontal axis. In the limit $\gamma ^\phi =0$ all modes are stable. Furthermore, changes in the anisotropy ratio shift the position of the maximum of the dispersion relation. If $\gamma ^\phi >\gamma ^z$ (figure 2b), the position of the maximum of $\omega$ shifts to larger values of $kR_0$, if $\gamma ^\phi <\gamma ^z$ the maximum is found at smaller $kR_0$ as shown in figure 2(c). This means for dominating axial tension $\gamma ^z$ the instability wavelength increases.

3.2. Anisotropic Rayleigh–Plateau instability for an ideal fluid

Performing again a linear stability analysis using the same solution procedure as before, we calculate the dispersion relation for an ideal fluid jet with same density inside and outside the jet as detailed in appendix B.2. We obtain the dispersion relation

(3.2)\begin{equation} \omega^2 = \frac{\gamma^\phi}{\rho R_0^3} (kR_0)^2 \left( 1 - \frac{\gamma^z}{\gamma^\phi} (kR_0)^2 \right) \textrm{I}_1(kR_0) \textrm{K}_1(kR_0). \end{equation}

Compared to the dispersion relation for a Stokes fluid in (3.1), we here obtain an equation for the squared growth rate $\omega ^2$. According to the ansatz for the perturbed interface in (2.1), perturbations with real and positive $\omega$ will grow. Imaginary $\omega$ describe oscillatory perturbations of the surface which do not grow in time. Imaginary $\omega$ correspond to negative values of $\omega ^2$. Each positive $\omega ^2$ has a positive and negative solution $\omega$. The positive solution will grow while the negative is damped. Thus, we are interested in non-negative values of $\omega ^2$ which are obtained from (3.2). This leads to the same expression for the range of growing wavelengths as for the Stokes fluid in § 3.1 because the relevant factor in the dispersion relation is identical, i.e. the anisotropy ratio ${\gamma ^{z}}/{\gamma ^{\phi }}$ enters the equation in the same way. However, the prefactor of the growth rate changes $\omega _0^2 = {\gamma ^{\phi }}/({\rho R_0^3})$, it now depends on the density rather than on the viscosity and is the squared inverse of the capillary time based on the azimuthal tension $\gamma ^\phi$. Furthermore, the geometrical factor containing the Bessel functions is remarkably different.

The dispersion relation (3.2) for the ideal fluid is shown in figure 3(a–c) for same values of the anisotropy ratio as in figure 2(a–c). We again observe a strong variation of the maximum position and range of unstable modes with changing anisotropy ratio. A remarkable difference to the Stokes fluid is the shape of the $\gamma ^{\phi }$ contribution. Here for an ideal fluid the destabilising $\gamma ^{\phi }$ contribution no longer reaches a maximum at finite $kR_0$ but instead increases indefinitely. Thus, in the limit $\gamma ^z=0$ all modes are unstable with steadily increasing growth rate. Compared to the Stokes fluid, the total dispersion relation is furthermore more asymmetric between $kR_0 = 0$ and $\sqrt {{\gamma ^{\phi }}/{\gamma ^{z}}}$ and the maximum shifts towards larger wavenumbers (compare e.g. figure 2b to figure 3b).

Figure 3. Dispersion relation for an ideal fluid with $\rho =\rho ^o$. Curves are shown for (a) isotropic interfacial tension ${\gamma ^z}/{\gamma ^\phi } = 1.0$ and for anisotropic interfacial tension with (b) ${\gamma ^z}/{\gamma ^\phi } = 0.5$ and (c) ${\gamma ^z}/{\gamma ^\phi } = 2.0$. We distinguish the contributions from $\gamma ^\phi$ (green) and $\gamma ^z$ (orange). While $\gamma ^z$ is purely damping, $\gamma ^\phi$ is destabilising. An anisotropic tension strongly alters the range of growing modes and shifts the maximum towards larger $kR_0$ in (b) or smaller $kR_0$ in (c).

In appendix E we further generalise our results to the dispersion relation including a general density and viscosity contrast as derived by Tomotika (1935).

4. Quantitative analysis of the effects due to tension anisotropy

4.1. Dominant wavelength

Having determined the range of (un)stable wavenumbers in the previous section, we now explicitly investigate how tension anisotropy affects the value of the dominant, i.e. fastest growing, wavelength $\lambda _{m}$. This quantity is of practical interest as it determines the size of the fragmented vesicles/droplets and provides an intrinsic length scale of the instability. For arbitrary ratios ${\gamma ^{z}}/{\gamma ^{\phi }}$, we use Mathematica to determine numerically the maximum of the dispersion relation (3.1) and (3.2). We further perform fully three-dimensional simulations of a membrane endowed with interfacial tension using BIM and LBM/IBM, as detailed in § 2.4. While BIM intrinsically solves the fluid dynamics in the Stokes limit, LBM/IBM simulations are run for $\textit {Oh} \approx 0.00025$, i.e. very close to the ideal fluid limit.

The result of the linear stability analysis is compared to the simulation results in figure 4. The solid lines show how the position of the maximum of the dispersion relations (3.1) (orange line in figure 4a) and (3.2) (red line in figure 4b), i.e. the dominant wavelength $\lambda _{m}$, changes with the ratio ${\gamma ^{z}}/{\gamma ^{\phi }}$. The simulation results for the Stokes fluid using BIM are drawn as triangles, those for the ideal fluid using LBM/IBM as squares. Both are in very good agreement with the respective theoretical predictions. For the Stokes fluid the obtained value $k_m R_{0} \approx 0.562$ (i.e. ${\lambda _m}/{R_0} \approx 11.18$) for isotropic tension ${\gamma ^{z}}/{\gamma ^{\phi }}=1$ is in good agreement with Tomotika (1935) and Stone & Brenner (1996). For the ideal fluid the dominant wavelength is smaller compared to the Stokes limit, which is true for all values of $kR_0$. For both Stokes fluid and ideal fluid increasing the ratio ${\gamma ^{z}}/{\gamma ^{\phi }}$ leads to an increase in the wavelength compared to its value at isotropic tension (illustrated by the insets from simulations on the right-hand sides of figures 4a and 4b). For decreasing ratio the opposite happens, the wavelength decreases (see inset from simulations at the top of figures 4a and 4b). Over the entire range of interfacial tension ratios we observe a nonlinear dependence of the wavelength on the anisotropy ratio. In the Stokes regime at an anisotropy ratio of zero a finite wavelength dominates. This is due to the fact that the $\gamma ^{\phi }$-contribution to the dispersion relation, as shown in figure 2(d), does not diverge for large wavenumbers but rather has a maximum at $k_m R_0=1.59$ which corresponds to a wavelength of $\lambda _m = 3.96 R_0$. This value matches the $y$-axis intercept of the dominant wavelength in figure 4(b). For the ideal fluid, however, for vanishing anisotropy ratio the wavelength goes to zero. In the limit of infinite ratio the wavelengths tend to infinity for both Stokes and ideal fluid.

Figure 4. Dominant wavelength as function of the anisotropy in interfacial tension. (a) Simulation results from BIM for the Stokes fluid are in very good agreement with the analytical results obtained from the dispersion relation (3.1). (b) Results for the ideal fluid from LBM/IBM agree very well with dominant wavelength obtained from the analytical dispersion relation (3.2). While the whole curve is at larger values in the Stokes limit, in both cases the dominant wavelength increases steadily with increasing ${\gamma ^z}/{\gamma ^\phi }$. Simulation snapshots of the interface are shown for different ratios ${\gamma ^{z}}/{\gamma ^{\phi }}$ over a length of about $55 R_0$ as insets.

In order to explain the effect of anisotropic interfacial tension, we first recall the classical Rayleigh–Plateau mechanism where two opposing effects influence the break-up. Since the radius in the region of a constriction is smaller than in a peak region, a pressure gradient develops pushing fluid out of the constriction and thus amplifying the disturbance. At the same time, however, due to the perturbation of the surface, the radius of curvature along the $z$-direction is negative in the region of the constriction and positive in the region of a peak. As can be seen from the Young–Laplace equation (2.4) this introduces another pressure gradient dragging the liquid back from the peak regions thus counteracting the pressure difference due to variations of the radius. The instability is a result of the interplay of both effects. An anisotropic interfacial tension weights these effects by either the azimuthal tension $\gamma ^\phi$ or axial tension $\gamma ^z$. Thus, a change in the ratio of the interfacial tension leads to a change in the weighting, shifting the region of growing wavelength and also altering the most unstable wavelength. This argument is illustrated by the three cases of the dispersion relation with its different contributions shown in figures 2 and 3.

The limit $\lambda _{m}\rightarrow \infty$ for $\gamma ^z \rightarrow \infty$ can be understood on the basis of the Young–Laplace equation for anisotropic interfacial tension in (2.4), as well. For infinite $\gamma ^z$ a finite curvature along $z$ would result in an infinite pressure difference. Thus, the interfacial tension must be balanced by a vanishing curvature, i.e. by an infinite curvature radius, which is equivalent to an infinite wavelength.

Finally, we discuss the limit $\gamma ^z\rightarrow 0$. We start with considering the ideal fluid. Due to finite $\gamma ^\phi$ every circular segment of the interface along the cylinder axis tends to contract. The incompressibility of the liquid inside prevents this homogeneous contraction. This means that a volume conserving neck–tail perturbation between two neighbouring thin circular segments, thus with very small wavelength and very large curvature in the $z$-direction, can in principle be established. Since $\gamma ^z\rightarrow 0$ there is no counteracting contribution which balances this tendency. The monotonic increase of the growth rate for the $\gamma ^\phi$-contribution with increasing $k$, as shown by the course of the dispersion relation in figure 3, suggests that such a perturbation grows fastest. This, in total, results in $\lambda _{m} \rightarrow 0$ for the ideal fluid. In the Stokes limit, however, the $\gamma ^{\phi }$-contribution is finite for large $kR_0$, possibly due to viscous stresses from the fluid which have a damping effect on the perturbation.

4.2. Transition between the two regimes

In the following, we will show that the border between the Stokes regime and the ideal fluid limit is not necessarily clear cut. In fact, by varying nothing more than the anisotropy ratio, the system can undergo a transition from one regime to the other. To demonstrate this transition, we present simulations using the LBM/IBM for typical parameters of biological cells/vesicles. We choose an interfacial tension of $\gamma ^{\phi } = 10^{-4}\ \textrm {N}\ \textrm {m}^{-1}$, which is the cortical tension reported for neutrophils (Tinevez et al. 2009) and in the middle of the range of typical tensions reported by Winklbauer (2015). As typical diameter of the tubular vesicle we choose $2R_0 = 1\ \mathrm {\mu }\textrm {m}$ and for the surrounding fluid density $\rho =1000\ \textrm {kg}\ \textrm {m}^{-3}$ and viscosity $\eta =1.2\times 10^{-3}\ \textrm {Pa}\ \textrm {s}$.

The dominant wavelength is shown in figure 5(a) as a function of the anisotropy ratio with green dots. For a clear comparison we also show the Stokes fluid dispersion relation in orange and the relation for the ideal fluid in red. The numerical simulations exhibit a transition between the two curves. For a small anisotropy ratio we obtain a finite wavelength, which nearly matches the result in the Stokes regime. In the case of isotropic tension we obtain a dominant wavelength of $k_m R_{0} \approx 0.628$ (i.e. ${\lambda _{m}}/{R_0} \approx 10.005$), a value in between the one for the ideal fluid and the Stokes regime. For larger values of the anisotropy ratio we end up close to the curve for an ideal fluid.

Figure 5. Transition between both regimes. (a) LBM/IBM simulations with typical vesicle and cell parameters (green dots) show dominant wavelengths between the two curves obtained in the limit of a Stokes fluid (orange) and an ideal fluid (red). (b) The transition between the two regimes in the wavelength is accompanied by a strong variation in the Ohnesorge number with respect to the tension along the axis $z$, i.e. $\textit {Oh}_z$.

We explain the transition between both regimes by the varying Ohnesorge number, which is shown in figure 5(b). By varying the anisotropy ratio, either the Ohnesorge number $\textit {Oh}_{\phi }$ or $\textit {Oh}_{z}$ varies, while the other can be kept constant. Here, we keep $\textit {Oh}_{\phi } \approx 5.5$ shown by the dark-green downwards-pointing triangles in figure 5(b). Consequently, $\textit {Oh}_{z}$ changes from 30 to approximately 2.5, as shown by the light-green upwards-pointing triangles. The transition in the Ohnesorge number $\textit {Oh}_{z}$ is matched by the transition in the wavelength. At small anisotropy ratios with large Ohnesorge number, the wavelength is close to the analytical prediction for the Stokes equation. This is in good agreement with the Stokes equation having $\textit {Oh} \rightarrow \infty$. Towards larger anisotropy ratios the Ohnesorge number becomes smaller and the wavelength approaches the predictions for an ideal fluid. We thus conclude that finite inertia effects trigger the transition, even though the Ohnesorge number is still larger than one. Our results clearly show that finite inertia effects can alter the Rayleigh–Plateau instability of tubular vesicles, even though their micrometric dimensions may at first sight suggest the opposite.

4.3. Dominant growth rate

We now investigate the growth rate of the most unstable mode $\omega _{m}$, i.e. the value of the maximum of the dispersion relation, in a quantitative manner. Using Mathematica we determine the maximum growth rate for varying tension anisotropy from the analytical dispersion relation. In addition, we perform simulations as described in § 2.4, where we extract the growth rate as described in appendix A.

The results in the limit of a Stokes fluid are shown in figure 6(a) and the limit of an ideal fluid is shown in figure 6(b). With increasing tension anisotropy the dominant growth rate decreases strongly. This can again be explained by the stabilising nature of the axial tension $\gamma ^z$ which slows down the instability. For tension anisotropy approaching infinity the growth rate approaches zero. At tension anisotropy equal zero we observe a finite growth rate of $\omega _{m} \approx 0.087$ for the Stokes fluid in (a), which results from viscous stresses in the Stokes fluid. In stark contrast, for an ideal fluid in (b) the maximum growth rate increases more strongly and even diverges for tension anisotropy to zero due to the destabilising nature of the azimuthal tension $\gamma ^{\phi }$. Simulation results for the Stokes fluid with BIM (triangles) and for the ideal fluid with LBM/IBM (squares) are in very good agreement with the corresponding analytical results. From the inverse of the growth rate the linear break-up time can be estimated, which is the time it takes until a droplet or vesicle pinches off. From figure 6 we can conclude that with decreasing tension anisotropy the break-up of the interface is strongly accelerated.

Figure 6. Growth rate of the dominant mode as a function of the anisotropy in interfacial tension. The dominant growth rate according to the dispersion relation (a) for a Stokes fluid (3.1) (orange line) and (b) for an ideal fluid (3.2) (red line) is shown with corresponding BIM simulations (triangles) and LBM/IBM simulations (squares), respectively, depending on the tension anisotropy ${\gamma ^z}/{\gamma ^\phi }$. The dominant growth rate decreases steadily and strongly with increasing tension anisotropy. While the decrease with increasing anisotropy is similar, the growth rate is one order of magnitude larger for the ideal fluid and it does not remain finite at zero anisotropy in contrast to the Stokes fluid in (a). In both cases simulation results are in perfect agreement with the theory.

4.4. Nonlinear correction to the linear break-up time

After discussing dispersion relation, growth rate and dominant wavelength obtained by linear stability analysis, we now proceed to investigate the nonlinear behaviour of the Rayleigh–Plateau instability. This is covered by the simulations presented above using BIM for a Stokes fluid and LBM/IBM for an ideal fluid. In the following, we extract the nonlinear correction of the linear break-up time (Ashgriz & Mashayek 1995; Martínez-Calvo et al. 2020), which we define based on (A 1) by

(4.1)\begin{equation} {\rm \Delta} t_{nl} = t_{b} - \frac{\ln(\epsilon_0^{-1})}{\omega_{m}}, \end{equation}

from simulations with initial perturbation amplitude $\epsilon _0$ as described in appendix A. This correction compares the break-up time obtained from simulations $t_{b}$ to the linear break-up time obtained from the maximum growth rate of the dispersion relation $\omega _{m}$ and is shown in figure 7 in relation to $t_{b}$. In the limit of a Stokes fluid the nonlinear correction varies strongly: we observe a change of sign above an anisotropy ratio of about ${\gamma ^z}/{\gamma ^\phi }=0.5$ where the correction becomes strongly negative. The nonlinear correction for the ideal fluid for the LBM/IBM is smaller, positive, and slightly increases with increasing tension anisotropy. For isotropic tension but varying Marangoni number of a surfactant-covered fluid jet Martínez-Calvo et al. (2020) report effects going in the same direction with a more pronounced variation of the nonlinear correction of the linear break-up time towards the Stokes limit and less variation towards the ideal fluid limit.

Figure 7. Nonlinear correction of the linear break-up time for varying tension anisotropy. The nonlinear correction of the linear break-up time is shown relative to the break-up time $t_{b}$ obtained from simulations. In the limit of an ideal fluid the LBM/IBM simulations show a slightly increasing nonlinear correction to the linear break-up time with increasing tension anisotropy. In contrast, BIM simulations show the reversed behaviour for a Stokes fluid, where in addition the sign changes and the amplitude variations are more pronounced.

5. Influence of interface viscosity

We now investigate how anisotropic interfacial tension influences the instability wavelength and growth rate in the Stokes regime if the interface in addition possesses interface viscosity (Boussinesq 1913; Scriven 1960; Whitaker 1976; Hajiloo, Ramamohan & Slattery 1987; Powers 2010; Yazdani & Bagchi 2013; Narsimhan et al. 2015; Martínez-Calvo & Sevilla 2018; Guglietta et al. 2020). The dispersion relation in presence of interface viscosity and tension anisotropy is derived in appendix B.3 and reads

(5.1)\begin{align} \omega &= \frac{\gamma^\phi}{\eta R_0} \left(1 - \frac{\gamma^z}{\gamma^\phi} (kR_0)^2 \right) \left( 1+ \frac{2\eta_{S}}{\eta R_0} kR_0 \xi \right)^{-1} \left[\vphantom{\frac{kR_0}{2}} {\rm I}_1(kR_0){\rm K}_1(kR_0) \right. \nonumber\\ &\left.\quad + \frac{kR_0}{2}\left( {\rm I}_1(kR_0) {\rm K}_0(kR_0) - {\rm I}_0(kR_0) {\rm K}_1(kR_0) \right)\right] \end{align}

with the abbreviation

(5.2)\begin{equation} \xi = - \tfrac{1}{2} \left[ kR_0 {\rm I}_1 (kR_0){\rm K}_1(kR_0) - kR_0 {\rm I}_0(kR_0) {\rm K}_0(kR_0) \right]. \end{equation}

For vanishing interface viscosity $\eta _{S}$ the dispersion relation correctly reduces to (3.1).

From the dispersion relation we numerically calculate the dominant wavelength $\lambda _{m}$ and the corresponding growth rate $\omega _{m}$ for varying anisotropy ratio ${\gamma ^z}/{\gamma ^\phi }$ and interface viscosity $\eta _{S}$, which we measure relative to fluid viscosity $\eta$. Phase diagrams of the dominant wavelength and the corresponding growth rate are shown in figures 8(a) and 8(b), respectively. Typical values for the interface viscosity of vesicles and red blood cells reported in the literature vary from $10^{-10}$ Pa s m to $10^{-7}$ Pa s m (Dimova et al. 2006; den Otter & Shkulipa 2007; Guglietta et al. 2020). In combination with typical fluid viscosities (§ 2.4) and sizes of the order of $R_0=1\ \mathrm {\mu }\mathrm {m}$ for vesicles and $R_0=4\ \mathrm {\mu }\mathrm {m}$ for red blood cells this leads to typical viscosity ratios ${2\eta _{{S}}}/({\eta R_0}) = 0.1 \text {--} 40$. In figure 8 these values correspond to the range from approximately $-1$ to approximately 1.6 on the ordinate. For fixed values of the viscosity ratio, an increase in the anisotropy ratio leads to a larger dominant wavelength and smaller maximum growth rate. This is in line with the results discussed above for zero interface viscosity in §§ 4.1 and 4.3. Especially in the region of small interface viscosities, changes in $\eta _{S}$ have little effect on dominant wavelength and growth rate, but also around a fraction of 1 the anisotropy dominates. If the interface viscosity $\eta _{S}$ becomes very large compared to the fluid viscosity, the growth rate tends to zero, as a more viscous interface leads to a slower dynamics of the perturbation. The latter is consistent with the results of Narsimhan et al. (2015) for isotropic interfaces. We also note that the maximum wavenumber of the instability, beyond which perturbations do not grow, is influenced only by the anisotropy ratio and not the viscosity ratio, because the root of the dispersion relation (5.1) is determined by the first term in brackets.

Figure 8. Influence of interface viscosity on the anisotropic Rayleigh–Plateau instability. Phase diagrams for (a) the dominant wavelength and (b) the corresponding growth rate. The dominant wavelength increases both with increasing anisotropy ratio and increasing interface viscosity $\eta _{S}$. The growth rate of the dominant perturbation decreases with increasing anisotropy ratio and interface viscosity. Despite the increase in the wavelength and the slowing down of the instability for very large values of the interface viscosity, the tension anisotropy is the dominating parameter.

6. Formation of satellite droplets for an ideal fluid jet without ambient fluid under influence of tension anisotropy

We eventually consider an ideal fluid jet without ambient fluid, i.e. $\eta ^0 = 0$ and $\rho ^0 \ll \rho$. In appendix C we derive the dispersion relation including anisotropic interfacial tension and show that for the ideal fluid jet without ambient fluid similar results hold as for the ideal fluid discussed in §§ 3.2 and 4. To compare our analytical results to simulations of an ideal fluid jet without ambient fluid we develop an axisymmetric simulation method based on a long-wavelength (or small-$k$) approximation as detailed in appendix D. We consider the Navier–Stokes equation (D 9) together with the kinematic boundary condition (D 10) in the small-$k$ approximation and solve them numerically. From the dynamic evolution the dominant growth rate and wavelength can be calculated, similarly to the procedure for the three-dimensional simulations detailed in appendix A. They are in good agreement with the theory as detailed in appendix C. We further obtain the jet shape over time.

With this, we are able to study the formation of satellite droplets under anisotropic interfacial tension. Satellite droplets are typically much smaller than and form in between the main drops during break-up of a liquid jet for certain parameter combinations (Eggers & Dupont 1994; Eggers 1997; Martínez-Calvo et al. 2020). In the presence of an ambient fluid (BIM and LBM/IBM simulations) we do not observe satellite droplet formation, which may be related to limitations in the resolution of our simulation and the fact that we do not simulate the actual pinch-off. In the following, we therefore investigate satellite droplet formation using the small-$k$ simulations for a liquid jet without ambient fluid at varying tension anisotropy and viscosity. Starting with the zero-viscosity, ideal fluid jet ($Oh=0$), we observe satellite droplets as shown at the bottom of figure 9. As the dominant wavelength increases with increasing anisotropy, the satellite droplet at the bottom right of figure 9 is more elongated compared to the left one. However, both satellites possess the same relative volume of approximately $\varXi =3$ % (Rutland & Jameson 1970; Lafrance 1975; Mansour & Lundgren 1990; Ashgriz & Mashayek 1995; Eggers 1997), which is defined as the volume integral over the satellite droplet divided by the volume integral over one period of the perturbation (Martínez-Calvo et al. 2020). For the ideal fluid jet, $\varXi$ is independent of the tension anisotropy.

Figure 9. Formation of satellite droplets under the influence of tension anisotropy in the absence of an ambient fluid $\rho ^o=0, \eta ^o=0$. Relative volume $\varXi$ of the satellite droplet for varying tension anisotropy ${\gamma ^z}/{\gamma ^\phi }$ and varying Ohnesorge number $Oh$. Parameter combinations for the shapes shown around the colour map are indicated in the phase diagram by black crosses. For the ideal fluid jet without ambient fluid with $\rho ^o=0$ ($Oh=0$) the relative volume of 3 % remains constant while for the Stokes fluid without ambient fluid $\eta ^o = 0$ (large $Oh$) no satellites appear. In the intermediate range, a significant influence of tension anisotropy on the relative volume is observed. Especially, the satellite droplet becomes cylindrical for large tension anisotropy, as shown in the top right image.

As our small-$k$ simulations are based on the full Navier–Stokes equation, they allow us to explore the effect of fluid viscosity measured in terms of the Ohnesorge number $Oh$ on the satellite droplet. In the centre of figure 9 we show the relative volume $\varXi$ in a colour map. For isotropic interfacial tension, we observe a decrease in the satellite volume, which is in agreement with results reported by Martínez-Calvo et al. (2020) in the absence of surfactants. With increasing Ohnesorge number the phase diagram in figure 9 shows a growing influence of the anisotropy. For a small but non-zero $Oh\approx 0.2$, the volume assumes its maximum at ${\gamma ^z}/{\gamma ^\phi }=0$ and decreases with increasing anisotropy. At higher $Oh$, the effect is reversed and the volume strongly increases with increasing tension anisotropy. In contrast to $Oh=0$, at larger Ohnesorge number $Oh=0.9$ the shape of the satellite droplet is strongly influenced by an increase in tension anisotropy from ${\gamma ^z}/{\gamma ^\phi }=0.5$ to ${\gamma ^z}/{\gamma ^\phi }=3.5$ as illustrated by the images at the top of figure 9. For small anisotropy (upper left shape) the satellite droplet is small and spherical but most remarkably for larger tension anisotropy (upper right) it develops a more cylindrical shape, it is longer and larger. In both cases the satellite droplet is connected to the main droplets by a thin fluid string. At very large Ohnesorge number beyond one, i.e. towards the Stokes regime, the satellite volume decreases to zero over the whole range of anisotropy. All in all, over a broad range of intermediate Ohnesorge numbers we observe a striking influence of tension anisotropy on the satellite droplet, where at larger $Oh$ larger tension anisotropy stabilises the satellite droplet.

7. Conclusion

Using linear stability analysis supported by numerical simulations we generalised the Rayleigh–Plateau mechanism for the break-up of a liquid cylinder to situations where the interfacial tension is anisotropic. Two physically relevant situations were studied: a vesicle/biological cell in the limit of small Reynolds number (Stokes equation) and an ideal fluid in the limit of large Reynolds number (Euler equation). We found that anisotropic interfacial tension alters not only the range of growing perturbations but also strongly affects the dominant wavelength of the instability. If the axial tension is inferior/superior to the azimuthal tension, the dominant wavelength becomes smaller/larger than the wavelength of the classical isotropic Rayleigh–Plateau instability. For strong anisotropy, the dominant wavelength can even lie outside the instability range of the classical isotropic Rayleigh–Plateau instability. The predictions of our linear stability analysis were found to be in excellent agreement with numerical simulations using a boundary integral method at low and a lattice-Boltzmann/immersed boundary method at high Reynolds number. LBM/IBM simulations with typical vesicle/cell parameters, which include finite inertia effects, show a transition in the wavelength for decreasing Ohnesorge number from the Stokes regime to the ideal fluid. Using simulations, we found the nonlinear correction to the linear break-up time to decrease/increase with the anisotropy ratio in the Stokes/Euler limits, respectively. We found that including viscosity of the interface surrounded by Stokes fluid leads to an increase in the dominant wavelength and a slower dynamics of the perturbation. However, except at very large interface viscosity, changes are small compared to the effects of varying tension anisotropy. Finally, we showed that the satellite droplet volume of 3 % is not affected by varying tension anisotropy in the limit of an ideal fluid. For Ohnesorge numbers just below one, the satellite volume decreases with increasing anisotropy.

In Part 2 we investigate the alteration of the anisotropic Rayleigh–Plateau instability due to the bending and shear elasticity present in a typical cell. Besides biological cells, our results apply to synthetic interfaces with anisotropic properties such as nematic liquid crystals confined to interfaces (Keber et al. 2014) or cell-laden hydrogels extruded from a nozzle during bioprinting applications (Snyder et al. 2015; Mandrycky et al. 2016).