## 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 Reference Plateau1873). 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 Reference Rayleigh1878). Later, he extended the theoretical description of this Rayleigh–Plateau instability to viscous jets in the inertialess Stokes limit (Rayleigh Reference Rayleigh1892). 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 (Reference Tomotika1935). A simplified version for the important case of equal viscosities has been presented by Stone & Brenner (Reference Stone and Brenner1996). 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 (Reference Eggers and Villermaux2008) for further details.

However, pearling and break-up due to the Rayleigh–Plateau mechanism are not restricted to liquid jets. In Reference Bar-Ziv and Moses1994, Bar-Ziv & Moses (Reference Bar-Ziv and Moses1994) 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 Reference Seifert1997). Under local application of laser tweezers the vesicle formed pearls (Bar-Ziv & Moses Reference Bar-Ziv and Moses1994). Using a hydrodynamic theory Nelson, Powers & Seifert (Reference Nelson, Powers and Seifert1995) and Goldstein *et al.* (Reference Goldstein, Nelson, Powers and Seifert1996) 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.* Reference Goldstein, Nelson, Powers and Seifert1996; Bar-Ziv, Tlusty & Moses Reference Bar-Ziv, Tlusty and Moses1997). Later, several experimental studies demonstrated different ways to induce the tension which is required for the pearling instability (Powers Reference Powers2010). Pulling on membrane tethers with optically trapped particles (Bar-Ziv, Moses & Nelson Reference Bar-Ziv, Moses and Nelson1998; Powers, Huber & Goldstein Reference Powers, Huber and Goldstein2002), protein mediated anchoring of membrane tethers to a substrate (Bar-Ziv *et al.* Reference Bar-Ziv, Tlusty, Moses, Safran and Bershadsky1999), applying a magnetic field (Ménager *et al.* Reference Ménager, Meyer, Cabuil, Cebers, Bacri and Perzynski2002), electric field (Sinha, Gadkari & Thaokar Reference Sinha, Gadkari and Thaokar2013) or osmotic pressure gradient (Yanagisawa, Imai & Taniguchi Reference Yanagisawa, Imai and Taniguchi2008; Sanborn *et al.* Reference Sanborn, Oglecka, Kraut and Parikh2013) can all lead to pearling. Furthermore, Kantsler, Segre & Steinberg (Reference Kantsler, Segre and Steinberg2008) 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 Reference Pal and Khakhar2019). Boedec, Jaeger & Leonetti (Reference Boedec, Jaeger and Leonetti2014) 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 (Reference Narsimhan, Spann and Shaqfeh2015) 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.* Reference Kruse, Joanny, Jülicher, Prost and Sekimoto2005; Marchetti *et al.* Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013; Prost, Jülicher & Joanny Reference Prost, Jülicher and Joanny2015; Salbreux & Jülicher Reference Salbreux and Jülicher2017; Jülicher, Grill & Salbreux Reference Jülicher, Grill and Salbreux2018). Such cytoskeletal networks can build a thin layer that underlies the plasma membrane, named the cell cortex (Alberts *et al.* Reference Alberts, Johnson, Lewis, Raff, Roberts and Walter2007; Köster & Mayor Reference Köster and Mayor2016; Chugh & Paluch Reference Chugh and Paluch2018), in which the action of motor proteins leads to active tension at the cell's interface (Chugh *et al.* Reference Chugh, Clark, Smith, Cassani, Dierkes, Ragab, Roux, Charras, Salbreux and Paluch2017). A positive constant active tension caused by a homogeneous (Pleines *et al.* Reference Pleines, Dutting, Cherpokova, Eckly, Meyer, Morowski, Krohne, Schulze, Gachet and Debili2013) actomyosin distribution in the cortex describes the internal tendency of the cytoskeleton to contract (Needleman & Dogic Reference Needleman and Dogic2017). Alternatively, proteins which anchor at the plasma membrane can trigger a pearling instability (Tsafrir *et al.* Reference Tsafrir, Sagi, Arzi, Guedeau-Boudeville, Frette, Kandel and Stavans2001) by bending the membrane and thus inducing a non-zero curvature (Campelo & Hernández-Machado Reference Campelo and Hernández-Machado2007; Jelerčič & Gov Reference Jelerčič and Gov2015). For a viscous active surface Mietke *et al.* (Reference Mietke, Jemseena, Kumar, Sbalzarini and Jülicher2019*a*) and Mietke, Jülicher & Sbalzarini (Reference Mietke, Jülicher and Sbalzarini2019*b*) report a Rayleigh–Plateau instability with mechano-chemical regulation. For a biological tissue composed of multiple cells Hannezo, Prost & Joanny (Reference Hannezo, Prost and Joanny2012) provide an energy argument based on an effective surface tension. Berthoumieux *et al.* (Reference Berthoumieux, Maître, Heisenberg, Paluch, Jülicher and Salbreux2014) 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 (Reference Bächer and Gekle2019) confirmed the instability threshold predicted by Berthoumieux *et al.* (Reference Berthoumieux, Maître, Heisenberg, Paluch, Jülicher and Salbreux2014) 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.* Reference Mora, Phou, Fromental, Pismen and Pomeau2010) and the production of synthetic vesicles (Anna Reference Anna2016; Pal & Khakhar Reference Pal and Khakhar2019). 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 Reference Emsellem, Cardoso and Tabeling1998), as a driving mechanism behind mitochondrial fission (Gonzalez-Rodriguez *et al.* Reference Gonzalez-Rodriguez, Sart, Babataheri, Tareste, Barakat, Clanet and Husson2015) as well as for pathological shapes of blood vessels during vasoconstriction (Alstrøm *et al.* Reference Alstrøm, Eguíluz, Colding-Jørgensen, Gustafsson and Holstein-Rathlou1999).

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.* Reference Reymann, Boujemaa-Paterski, Martiel, Guerin, Cao, Chin, De La Cruz, Thery and Blanchoin2012; Murrell *et al.* Reference Murrell, Oakes, Lenz and Gardel2015; Blackwell *et al.* Reference Blackwell, Sweezy-Schindler, Baldwin, Hough, Glaser and Betterton2016; Zhang *et al.* Reference Zhang, Kumar, Ross, Gardel and de Pablo2018), e.g. due to the formation of stress fibres (Tojkander, Gateva & Lappalainen Reference Tojkander, Gateva and Lappalainen2012). Accordingly, the tension in the cell cortex can be anisotropic (Rauzi *et al.* Reference Rauzi, Verant, Lecuit and Lenne2008; Mayer *et al.* Reference Mayer, Depken, Bois, Jülicher and Grill2010; Behrndt *et al.* Reference Behrndt, Salbreux, Campinho, Hauschild, Oswald, Roensch, Grill and Heisenberg2012; Callan-Jones *et al.* Reference Callan-Jones, Ruprecht, Wieser, Heisenberg and Voituriez2016), which is important for many biological phenomena such as cell-shape regulation (Callan-Jones *et al.* Reference Callan-Jones, Ruprecht, Wieser, Heisenberg and Voituriez2016), cell polarisation (Mayer *et al.* Reference Mayer, Depken, Bois, Jülicher and Grill2010), ingression formation (Reymann *et al.* Reference Reymann, Staniscia, Erzberger, Salbreux and Grill2016), the formation of a furrow during cell division (White & Borisy Reference White and Borisy1983; Salbreux, Prost & Joanny Reference Salbreux, Prost and Joanny2009) and the production of blood platelets (Bächer, Bender & Gekle Reference Bächer, Bender and Gekle2020). For a solid rod in the absence of any kind of fluid, Gurski & McFadden (Reference Gurski and McFadden2003) 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 Reference Ashgriz and Mashayek1995; Martínez-Calvo *et al.* Reference Martínez-Calvo, Rivero-Rodríguez, Scheid and Sevilla2020) under anisotropic interfacial tension. In Part 2 (Bächer, Graessel & Gekle Reference Bächer, Graessel and Gekle2021) 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 Reference Eggers and Villermaux2008; Boedec *et al.* Reference Boedec, Jaeger and Leonetti2014), 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)$.

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 Reference Drazin and Reid2004). 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

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 Reference Green and Zerna1954; Barthès-Biesel Reference Barthès-Biesel2016; Salbreux & Jülicher Reference Salbreux and Jülicher2017). 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 Reference Green and Zerna1954; Salbreux & Jülicher Reference Salbreux and Jülicher2017). 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

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 Reference Chandrasekharaiah and Debnath1994). 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

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

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

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

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 Reference Eggers and Villermaux2008), which relates characteristic scales of the interface and the surrounding fluid

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

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.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 Reference Bächer and Gekle2019). 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 Reference Pozrikidis2001; Zhao *et al.* Reference Zhao, Isfahani, Olson and Freund2010). 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 (Reference Bächer and Gekle2019). 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.* Reference Dmitrieff, Alsina, Mathur and Nédélec2017). For details on the implementation of the BIM we refer to Guckenberger & Gekle (Reference Guckenberger and Gekle2018).

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 Reference Aidun and Clausen2010; Krüger *et al.* Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2016) together with the IBM (Peskin Reference Peskin2002; Mittal & Iaccarino Reference Mittal and Iaccarino2005; Bächer, Schrack & Gekle Reference Bächer, Schrack and Gekle2017; Mountrakis, Lorenz & Hoekstra Reference Mountrakis, Lorenz and Hoekstra2017; Bächer *et al.* Reference Bächer, Kihm, Schrack, Kaestner, Laschke, Wagner and Gekle2018). Our LBM/IBM is implemented in the software package ESPResSo (Limbach *et al.* Reference Limbach, Arnold, Mann and Holm2006; Roehm & Arnold Reference Roehm and Arnold2012; Arnold *et al.* Reference Arnold, Lenz, Kesselheim, Weeber, Fahrenberger, Roehm, Košovan and Holm2013; Weik *et al.* Reference Weik, Weeber, Szuttor, Breitsprecher, de Graaf, Kuron, Landsgesell, Menke, Sean and Holm2019) and has been extensively validated (Gekle Reference Gekle2016; Guckenberger *et al.* Reference Guckenberger, Schraml, Chen, Leonetti and Gekle2016; Bächer *et al.* Reference Bächer, Kihm, Schrack, Kaestner, Laschke, Wagner and Gekle2018; Bächer & Gekle Reference Bächer and Gekle2019).

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

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 (Reference Plateau1873). 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 Reference Drazin and Reid2004). 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.

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 2*c*). 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 2*b*), 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

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 2*b* to figure 3*b*).

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

## 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 4*a*) and (3.2) (red line in figure 4*b*), 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 (Reference Tomotika1935) and Stone & Brenner (Reference Stone and Brenner1996). 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 4*a* and 4*b*). For decreasing ratio the opposite happens, the wavelength decreases (see inset from simulations at the top of figures 4*a* and 4*b*). 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.

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.* Reference Tinevez, Schulze, Salbreux, Roensch, Joanny and Paluch2009) and in the middle of the range of typical tensions reported by Winklbauer (Reference Winklbauer2015). 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.

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.

### 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 Reference Ashgriz and Mashayek1995; Martínez-Calvo *et al.* Reference Martínez-Calvo, Rivero-Rodríguez, Scheid and Sevilla2020), which we define based on (A 1) by

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.* (Reference Martínez-Calvo, Rivero-Rodríguez, Scheid and Sevilla2020) 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.

## 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 Reference Boussinesq1913; Scriven Reference Scriven1960; Whitaker Reference Whitaker1976; Hajiloo, Ramamohan & Slattery Reference Hajiloo, Ramamohan and Slattery1987; Powers Reference Powers2010; Yazdani & Bagchi Reference Yazdani and Bagchi2013; Narsimhan *et al.* Reference Narsimhan, Spann and Shaqfeh2015; Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018; Guglietta *et al.* Reference Guglietta, Behr, Biferale, Falcucci and Sbragaglia2020). The dispersion relation in presence of interface viscosity and tension anisotropy is derived in appendix B.3 and reads

with the abbreviation

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.* Reference Dimova, Aranda, Bezlyepkina, Nikolov, Riske and Lipowsky2006; den Otter & Shkulipa Reference den Otter and Shkulipa2007; Guglietta *et al.* Reference Guglietta, Behr, Biferale, Falcucci and Sbragaglia2020). 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.* (Reference Narsimhan, Spann and Shaqfeh2015) 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.

## 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 Reference Eggers and Dupont1994; Eggers Reference Eggers1997; Martínez-Calvo *et al.* Reference Martínez-Calvo, Rivero-Rodríguez, Scheid and Sevilla2020). 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 Reference Rutland and Jameson1970; Lafrance Reference Lafrance1975; Mansour & Lundgren Reference Mansour and Lundgren1990; Ashgriz & Mashayek Reference Ashgriz and Mashayek1995; Eggers Reference Eggers1997), 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.* Reference Martínez-Calvo, Rivero-Rodríguez, Scheid and Sevilla2020). For the ideal fluid jet, $\varXi$ is independent of the tension anisotropy.

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.* (Reference Martínez-Calvo, Rivero-Rodríguez, Scheid and Sevilla2020) 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 Reference Bächer, Graessel and GeklePart 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.* Reference Keber, Loiseau, Sanchez, DeCamp, Giomi, Bowick, Marchetti, Dogic and Bausch2014) or cell-laden hydrogels extruded from a nozzle during bioprinting applications (Snyder *et al.* Reference Snyder, Rin Son, Hamid, Wang, Lui and Sun2015; Mandrycky *et al.* Reference Mandrycky, Wang, Kim and Kim2016).

## Acknowledgements

C.B. and K.G. thank the Studienstiftung des deutschen Volkes for financial support. C.B. acknowledges support by the study programme ‘Biological Physics’ of the Elite Network of Bavaria. We acknowledge funding from Deutsche Forschungsgemeinschaft in the framework of FOR 2688 ‘Instabilities, Bifurcations and Migration in Pulsating Flow’, Project B3. We gratefully acknowledge computing time provided by the SuperMUC system of the Leibniz Rechenzentrum, Garching, as well as by the Bavarian Polymer Institute, and financial support from the Volkswagen Foundation.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Simulation analysis

In the following we explain in detail the analysis procedure of the BIM and LBM/IBM simulations. In order to analyse in simulations both the wavelength of the dominant mode and its growth rate we consider the interface shape over time for a given perturbation and track a local maximum corresponding to a later droplet/vesicle as shown in figure 10(*a*). In order to obtain the most unstable mode for given anisotropic interfacial tension in the simulations, we consider a cylindrical interface of fixed anisotropic interfacial tension immersed in a fluid of fixed properties and apply a series of perturbations with varying wavelength. One simulation corresponds to one perturbation with fixed wavelength. The initial amplitude of the perturbation is chosen as $\epsilon (t=0) = 0.02 R_0$ in all simulations. Note that a multiple of the wavelength has to fit in the simulation box and we thus vary the wavelength in discrete steps. Figure 10(*a*) shows an example of a series of simulations for given anisotropic interfacial tension ${\gamma ^z}/{\gamma ^\phi } = 0.4$ using the LBM/IBM. In this example we force a perturbation with 10–17 maxima onto a cylindrical interface in a box of length 650 grid cells. Each curve corresponds to a simulation with different wavelength of the perturbation. In each of the simulations the radius at the position of the different maxima is tracked over time. The curves show the local radius at the position of the maxima, averaged over all maxima, varying in time. From the initial amplitude $\epsilon (t=0) = 0.02 R_0$ all modes grow with a different speed. In order to determine the fastest growing wavelength, we define a threshold for the amplitude $\epsilon _{{crit}}= 0.1 R_0$, which can be considered as an upper limit of small deformations. The inset of figure 10 shows the growth of all modes up to the radius $R_0 + \epsilon _{{crit}}$. The mode which first reaches this threshold, is considered as the fastest growing mode and its corresponding wavelength as the most unstable wavelength. This procedure can be applied to different values of the anisotropy ratio ${\gamma ^z}/{\gamma ^\phi }$ and allows us to numerically determine the most unstable wavelength depending on tension anisotropy. With increasing wavelength we increase the box length in order to ensure a good resolution for this procedure. Nevertheless, the finite box size and the resulting discrete variation of the wavelength result in discretisation effects. In order to account for these discretisation effects we compute the range in wavelength which cannot be distinguished due to the discrete variation: adding/subtracting half the wavelength is not possible for a given box length, because we necessarily enforce an integral number of wavelengths in the simulation box. This range is given in the figures as error bars of the simulation results.

The growth rate is calculated as the quotient of the derivative of the maximum radius over time and the maximum radius itself, in a regime of linear growth. In BIM simulations we use the series of perturbations as detailed above, for the LBM/IBM simulations we consider one period of the dominant perturbation with increased resolution and initial perturbation amplitude $\epsilon _0 = 0.002$. Figure 10(*b*) shows the growth rate determined from simulations for varying wavenumber in comparison to the analytically obtained dispersion relation (3.2) for LBM/IBM simulations. An error is estimated from the average value of the quotient which determines the growth rate.

As a measure of the nonlinear behaviour covered by the simulations, we further extract the nonlinear correction to the linear break-up time ${\rm \Delta} t_{nl}$ defined in (4.1). To do so we determine the break-up time in the simulation $t_{b}$ by tracking a minimum of the interface shape over time and extrapolating the last time steps towards a radius of zero. From the analytical evolution of a perturbation based on the linear stability analysis as given in (B 1) the linear break-up time can be calculated at the position of a minimum by

The difference between $t_\mathrm {b}$ and the linear break-up time determines the nonlinear correction as given by (4.1).

## Appendix B. Analytical derivation of the dispersion relations

#### B.1. Dispersion relation for a Stokes fluid

In the following we consider the inner and outer fluid in the Stokes limit with same viscosity, i.e. $\eta ^o = \eta$. The following calculations are based on the method introduced by Stone & Brenner (Reference Stone and Brenner1996). We use an ansatz of a cylindrical interface which is slightly perturbed in a periodic fashion with time dependent amplitude $\epsilon (t)\ll 1$. For later convenience we write the equation describing the shape as

which is a slightly different notation but equivalent to the ansatz in (2.1). We calculate the traction jump at the interface ${\rm \Delta} f_n$ in (2.2), which is equal to the negative membrane force, using (B 1) and assuming that the magnitude of the interface perturbation is small, i.e. $\epsilon \ll 1$,

The traction jump can be decomposed into a constant pressure contribution $p_0={\gamma ^\phi }/{R_0}$ and a perturbation in the traction jump

which is evaluated at the position of the unperturbed interface due to linearisation (Tomotika Reference Tomotika1935; Stone & Brenner Reference Stone and Brenner1996). The goal now is to solve the Stokes equation (2.8) in the presence of the traction jump perturbation (B 3). This is done by considering a ring force representing the traction jump (Stone & Brenner Reference Stone and Brenner1996) such that the Stokes equation becomes

The ring force in the radial direction (third term entering in (B 4) with the radial unit vector $\hat {\boldsymbol {e}}_r$) accounts for the presence of the membrane, from which a force due to interfacial tension is acting on the fluid. The interfacial force is evaluated at the undeformed radial position of the infinitely thin interface. As a consequence the ring force enters with a delta distribution $\delta (r-R_0)$.

In line with the periodic perturbation of the radius of the interface in (B 1), a periodic ansatz for the velocity components and the pressure field is chosen (Stone & Brenner Reference Stone and Brenner1996)

where, due to the perturbation of the interface (B 1) and kinematic boundary conditions, the radial velocity is written with a cosine while the axial velocity due to continuity equation (2.6), which involves a derivative with respect to $z$, is written with a sine. The prefactor of the velocities ${\gamma ^\phi }/{\eta }$ is chosen such that $\bar {v}_r$ and $\bar {v}_z$ become dimensionless, where the same applies for the pressure prefactor. The velocity prefactor ${\gamma ^\phi }/{\eta }$ is identical to the viscocapillary velocity based on the azimuthal tension $\gamma ^\phi$. The amplitude of the perturbation varies in time, i.e. $\epsilon =\epsilon (t)$. For solving the Stokes equation we introduce the Hankel transforms (Poularikas Reference Poularikas2000) of the velocity amplitudes $V_r(s)$, $V_z(s)$ and the pressure $P(s)$

together with their inverse transforms

where $\textrm {J}_\nu$ is the Bessel function of the first kind and order $\nu$.

For the transformation of the Stokes equation (B 4) and continuity equation (2.6) into Hankel space we use the following identities for the Bessel differential operators (Poularikas Reference Poularikas2000):

and

Together with the Hankel transform of the ring force

we are able to write down the continuity equation together with the two components of the Stokes equation in Hankel space

We can solve this system of equations for the radial velocity in Hankel space and obtain

which leads to the radial velocity in real space using the inverse Hankel transform

As done in the case of the ideal fluid jet without ambient fluid, we consider a perturbation of the interface growing with rate $\omega$, i.e.

The kinematic boundary condition at the interface leads to

and using (B 21) we obtain the growth rate

The integral can be evaluated (Gradshteĭn, Ryzhik & Jeffrey Reference Gradshteĭn, Ryzhik and Jeffrey2007, 6.535), taking the derivative with respect to variable $k$ in the denominator) and we obtain the dispersion relation for a Stokes fluid in (3.1).

#### B.2. Dispersion relation for an ideal fluid

Below we perform a linear stability analysis in the ideal fluid limit. Since the perturbation of the interface between both fluids is small, also the perturbation of the velocity is small. Because we consider a co-moving frame, the velocity vector contains only the perturbation and thus is small, $v_i \ll 1$, as well. As a consequence, the nonlinear term in the Euler equation is of second order and can be neglected, leading to the linear Euler equation

We here solve the linearised Euler equation (B 25) which is valid for each position $(r,z)$ both in the inner and in the outer fluid of same density $\rho ^o = \rho$. Note that the continuity equation does not change compared to the previous section.

In the linearised Euler equation we add the traction jump that here equals a pressure disturbance occurring at the interface due to interfacial tension according to (2.4), which is expanded as done in (B 1) and enters in terms of a ring force. In turn we have to solve

In the linearised Euler equation (B 25) the density appears rather than the viscosity. In the perturbation ansatz for the velocity (B 5) and (B 6) the prefactor changes accordingly to $\sqrt{ {\gamma ^\phi }/({\rho R_0})}$. The prefactors are chosen such that $\bar {v}_r(r)$, $\bar {v}_z(r)$ are again dimensionless. In total, after transformation into the Hankel space, we obtain the analytical equations

which are solved for the radial component of the velocity, as done in the case of the Stokes equation in § 3.1. We obtain for the velocity

and using the kinematic boundary condition we obtain the squared growth rate

We can again evaluate the integral (Gradshteĭn *et al.* Reference Gradshteĭn, Ryzhik and Jeffrey2007, 6.535) and the resulting dispersion relation for an ideal fluid is given in (3.2).

#### B.3. Dispersion relation taking interface viscosity into account

We now derive a dispersion relation for a system where, in addition to the fluid viscosity $\eta$, a viscosity of the interface $\eta _{S}$ is considered. The solution method starts from the Stokes equation and is again based on the previous B.1. Here, the interface viscosity results in an additional force acting from the interface onto the fluid (Scriven Reference Scriven1960; Narsimhan *et al.* Reference Narsimhan, Spann and Shaqfeh2015; Sprenger *et al.* Reference Sprenger, Shaik, Ardekani, Lisicki, Mathijssen, Guzmán-Lastra, Löwen, Menzel and Daddi-Moussa-Ider2020). The component normal to the interface is given by

where surface incompressibility and for the second identity the velocity ansatz from (B 6) is used. Similar to Powers (Reference Powers2010) we do not consider any tangential component of the viscous force. This additional normal force contribution appears as an additional term in the ring force in (B 4) such that (B 18) in the Hankel space becomes

Analogously, the radial velocity in Hankel space from (B 20) becomes

which is then transformed back to obtain $v_r$, which still depends on $\bar {v}_z(R_0)$. In order to obtain $\bar {v}_z(R_0)$, we use the continuity equation (B 17) to relate $V_z$ to $V_r$ and insert (B 34), then transform $V_z$ back from Hankel space and thus identify

which is then solved for $\bar v_z (R_0)$. Following the same steps as in appendix B.1 we then calculate the growth rate $\omega$:

The integral in this equation has already been solved in appendix B.1. Combining (B 36) and the solution for $\bar {v}_z (R_0)$ and solving the integral in (B 35) (based on 49(13) and 49(14) of Erdelyi *et al.* Reference Erdelyi, Magnus, Oberhettinger and Tricomi1954) leads to the dispersion relation in (5.1).

## Appendix C. Dispersion relation for an ideal fluid jet without ambient fluid

In deriving the dispersion relation for an ideal fluid jet without ambient fluid we follow the ansatz by Eggers & Villermaux (Reference Eggers and Villermaux2008). Here, we consider the outer fluid to have no influence, i.e. $\eta ^o=0$ and $\rho ^o \ll \rho$. As we perform a linear stability analysis, we again use the linearised version of the Euler equation (B 25). By the linearised Euler equation (B 25) and the continuity equation (2.6) the pressure $p$ must fulfil a Laplace equation

Corresponding to the interface perturbation, a perturbation ansatz for the pressure distribution in the jet is chosen (Eggers & Villermaux Reference Eggers and Villermaux2008)

with constant pressure $p_0$ and a general perturbation $\delta p$. At the interface, the pressure is determined by the modified Young–Laplace equation which according to (2.4) is

The pressure perturbation can be separated into the periodic perturbation in the $z$-direction, its magnitude $F(r)$ depending on the radial position and the constant prefactor $\delta \bar {p}$

The pressure obeys the Laplace equation (C 1), which in cylindrical coordinates and using that $p$ does not depend on the angular coordinate $\phi$, reads

Now, using (C 2) for the pressure and inserting (C 4) leads to a partial differential equation for the function $F(r)$

This equation is solved by the modified Bessel function of the first kind and of order zero

Using the perturbation ansatz (2.1) for the interface and considering anisotropic interfacial tension, the perturbation in the pressure in (C 3) is in analogy to (B 3) given by

From (C 4) and (C 7) we have the relation $\delta p=\delta \bar {p} \textrm {I}_0(kr) \exp ({\omega t+{\textrm {i}} kz})$, which together with (C 8) leads to the magnitude of the pressure perturbation

Since $\epsilon _0\ll R_0$ we can evaluate the modified Bessel functions at the position of the unperturbed interface $R_0$.

In line with the ansatz for the pressure a perturbation ansatz for the velocity is chosen,

with $\boldsymbol {v}_0$ the velocity of the unperturbed jet, which vanishes in the co-moving frame. The motion of the interface is governed by the kinematic boundary condition

The velocity perturbation is linked to the pressure by the linearised Euler equation (B 25), its $r$-component in cylindrical coordinates reads

Combination of the derivative with respect to time $t$ of (C 11), where in linearised form the second term on the left-hand side vanishes, and (C 12) gives

Inserting the ansatz for the interface radius (2.1) on the left and the one for the pressure (C 2) on the right-hand side of (C 13) and using the fact that $\textrm {I}^\prime _0(x) = \textrm {I}_1(x)$ leads to

Evaluating the modified Bessel functions at $R_0$ and inserting (C 9) finally yields the dispersion relation for an ideal fluid jet without ambient fluid

This dispersion relation is shown in figures 11(*a*)–11(*c*) for different values of the anisotropy ratio. The general shape of the dispersion relation and the changes due to variation of the anisotropy ratio are similar to that discussed in § 3.2.

In figures 11(*d*) and 11(*e*) dominant wavelength and corresponding maximal growth rate for the ideal fluid jet without ambient fluid are drawn as blue lines. For isotropic tension ${\gamma ^{z}}/{\gamma ^{\phi }} = 1$ we recover the well-known result (Rayleigh Reference Rayleigh1878; Drazin & Reid Reference Drazin and Reid2004; Eggers & Villermaux Reference Eggers and Villermaux2008) for the dominant wavelength of $k_{m} R_0\approx 0.697$ (i.e. ${\lambda _{m}}/{R_0}\approx 9.015$). Simulations for the ideal fluid jet without ambient fluid based on the small-$k$ approximation are included as blue triangles. They agree well with the analytical results. For comparison the curves and simulations for the ideal fluid from 4(*b*) and 6(*b*) are shown as red lines/squares. In figure 11(*d*) we see that the curves are quite similar, i.e. the dominant wavelength changes only slightly when there is no ambient fluid. The corresponding growth rate is significantly influenced by an additional ambient ideal fluid at small and intermediate tension anisotropy, while the influence vanishes for large anisotropy as we can observe in figure 11(*e*).

## Appendix D. Long-wavelength description

#### D.1. Dynamic equations for jet simulations

In the following, we derive the fluid equations of motion in long-wavelength (or small-$k$) approximation (Weber Reference Weber1931; Eggers & Dupont Reference Eggers and Dupont1994; García & Castellanos Reference García and Castellanos1994; Eggers Reference Eggers1997) for a fluid jet without ambient fluid, i.e. $\eta ^o = \rho ^o = 0$. The final equations are then solved numerically as detailed below to obtain the dynamic nonlinear evolution of an ideal jet interface in absence of an ambient fluid. The perturbation of the interface is considered to have a wavelength considerably longer than the radius of the liquid jet. Therefore, a typical radial length is small compared to a typical axial length, i.e. terms with a dependency on the radial position enter with a factor $\hat \epsilon \ll 1$ (Eggers & Villermaux Reference Eggers and Villermaux2008). The axial velocity is expanded in radial direction using a Taylor series and considering axisymmetry

From the continuity equation (2.6) one can obtain the radial velocity

where a prime denotes the partial derivative with respect to $z$. Similarly, we consider an expansion of the pressure with respect to the radial position

We further calculate the force from the inner fluid of the jet onto the interface using the outward unit normal vector on the interface given by

with $\hat {\boldsymbol {e}}_z$ the unit vector along the jet axis. The normal component of the force is given by the projection of the three-dimensional fluid stress twice onto the normal vector (Chandrasekharaiah & Debnath Reference Chandrasekharaiah and Debnath2014)

and the tangential force evaluated at the interface $r=R$ is

with $\boldsymbol {t}$ pointing along the interface which is perpendicular to $\boldsymbol {n}$ in the plane. Combining the force from the internal fluid and the force due to anisotropic interfacial tension the pressure at the interface follows

In addition, we obtain from the tangential force across the interface

The Navier–Stokes equation determining the fluid dynamics finally becomes

The kinematic boundary condition is also considered up to terms in lowest order with respect to $\hat \epsilon$, which gives the dynamics of the interface

The closed system of coupled equations (D 9) and (D 10) can then be solved numerically. We consider a one-dimensional interface which is discretised by 220 points. We initially perturb the interface according to (D 11) with $\epsilon _0=0.001$ and initialise the velocity profile $v_0(z)$ with zero. Derivatives of the radius and velocity profile are calculated at each time step using quintic splines. This allows us to calculate (D 8) and the right-hand side of (D 9) and (D 10). For solving the equations a three step Runge–Kutta algorithm is used with a typical time step of ${\rm \Delta} t = 0.00025$.

#### D.2. Dispersion relation in the long-wavelength description

The aim of the following discussion is to derive a simple, analytical equation for the wavenumber and growth rate of the most unstable perturbation mode. This we use to prescribe the dominant mode in small-$k$ simulations of a fluid jet with varying Ohnesorge number. As ansatz for the periodic interface perturbation growing in time we use

and for the velocity

Inserting the ansatz into the kinematic boundary condition (D 10) leads to

in linear order in $\epsilon _0$. The Navier–Stokes equation leads to

Combining both equations results in the dispersion relation for a fluid jet without ambient fluid in the small-$k$ approximation

which now depends on tension anisotropy, fluid density and viscosity.

The dispersion relation in the small-$k$ approximation (D 15) can be solved analytically for the position of the maximum, i.e. the dominant mode

and the dominant growth rate is determined by

In the limit of an ideal fluid jet, i.e. for vanishing viscosity and vanishing Ohnesorge number we obtain for the maximum growth rate

which approximates the blue curve in figure 11(*e*). We observe that in the framework of the small-$k$ approximation the dominant growth rate scales with the tension anisotropy to the power of minus one half. Analogously, the dominant wavelength in the limit of an ideal fluid jet scales with the square root of the tension anisotropy, as can be seen from (D 16). We note that the full dispersion relation with contributions of the Bessel functions leads to deviations from this scaling behaviour.

In figure 12 we compare the full dispersion relation obtained analytically for an ideal fluid jet without ambient fluid (C 15) to the one obtained in the small-$k$ limit (D 15) for vanishing viscosity. The inset shows the two dispersion relations for an anisotropy ratio of ${\gamma ^z}/{\gamma ^\phi } = 0.5$. We analyse the deviation of the approximation quantitatively by calculating the squared difference averaged over all sample points relative to the maximum of the dispersion relation

as the error. We observe a strong increase of the error towards small anisotropy. Nevertheless, the deviation for ${\gamma ^z}/{\gamma ^\phi }=0.1$ is approximately 15 % and thus the approximation still reasonable accurate. Increasing the anisotropy results in the deviation approaching zero, i.e. the approximation becomes perfectly accurate towards large tension anisotropy.

## Appendix E. General dispersion relation

We use the insights from the detailed derivations of the dispersion relations in the appendices B and C to obtain the modifications due to anisotropic interfacial tension of more general dispersion relations. From the dispersion relation in presence of an ambient Stokes fluid (3.1) or an ambient ideal fluid (3.2) we follow Tomotika (Reference Tomotika1935) and introduce a general density $N_\rho = {\rho }/{\rho ^o}$ and viscosity contrast $N_\eta = {\eta }/{\eta ^o}$ between inner and outer fluid. We use the following abbreviations concerning the wavenumber and the modified Bessel functions

*a*–

*c*)\begin{equation} y^2 = (kR_0)^2 + \omega \frac{\rho}{\eta} R_0^2 \quad F(x)=x \frac{{\rm I}_0(x)}{{\rm I}_1(x)} \quad G(x)=x\frac{{\rm K}_0(x)}{{\rm K}_1(x)}, \end{equation}

where an additional superscript $o$ indicates definition with the corresponding parameters of the outer fluid. The dispersion relation in presence of an ambient medium with a density and viscosity contrast (Tomotika Reference Tomotika1935) including anisotropic interfacial tension then takes the form

This general dispersion relation includes the Stokes fluid in (3.1) as well as the dispersion relation for an ideal fluid (3.2) in the corresponding limits. These are obtained by first considering identical fluids inside and outside, i.e. $N_\rho = 1$ and $N_\eta =1$. Second, for the Stokes fluid an expansion of the Bessel functions in the limit of small density is necessary, while for the ideal fluid the viscosity is set to zero and identities for the Bessel functions are used to rewrite the remaining terms.

Starting from the dispersion relation for an ideal fluid jet without ambient fluid (C 15) we follow Chandrasekhar (Reference Chandrasekhar1961) to obtain the general dispersion relation for a jet in a passive ambient medium. The general dispersion relation including tension anisotropy in this case has the form

This contains the dispersion relation of an ideal fluid jet in (C 15) in the limit of negligible viscosity.