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.


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 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 η or the membrane of a vesicle or cell immersed in a fluid in the limit of the Stokes equation, i.e. density ρ = 0. The fluid jet is immersed in an ambient fluid with η 0 , ρ 0 . The cylindrical interface of initial radius R 0 (dashed line) is subjected to a periodic perturbation with amplitude (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 γ z (orange) different from that in the azimuthal direction γ φ (green), both of which contribute to the membrane force acting onto the fluid with different curvature components (grey circles).
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.

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 δR(z, t), such that the local radius is given by R(z, t) = R 0 + δ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 2004). The perturbation of the interface is illustrated in figure 1: it modulates the radius in z-direction along the cylinder axis with amplitude (t) = 0 e ωt , a wavelength λ and a corresponding wavenumber k = 2π/λ of the wave vector pointing along the cylinder axis. The perturbation with initial amplitude 0 grows in time with growth rate ω. Accordingly, the interface of the jet can be described by its radius as (2.1) 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.

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 γ φ and the axial interfacial tension γ 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 φ is the inverse of the local radius of the interface R φ = 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 such that the curvature is negative at a neck and positive at a bulge (compare figure 1). The derivative R 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 Δ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 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 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 γ /R = p in the limit of isotropic surface tension γ = γ φ = γ z and vanishing curvature along z, i.e. R z → ∞. 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.

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 v, fluid density ρ, kinematic viscosity ν = η/ρ and shear viscosity η. Density and viscosity of the outer fluid are denoted by ρ o and η 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 ∇ · v = 0. (2.6) The Navier-Stokes and continuity equations together govern the motion of the fluid. Besides the well-known Reynolds number Re = ρR 0 V 0 /η 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 In our anisotropic scenario with γ z / = γ φ we define two distinct Ohnesorge numbers Oh z and Oh φ 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 0 = −∇p + ηΔv. (2.8) 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 (2.9) 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.

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 ρ 0 = ρ or viscosity η 0 = η. 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.

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 Re = 0 and an Ohnesorge number Oh → ∞. Membrane forces due to interfacial tension are calculated as detailed in Bächer & Gekle (2019). As the interfacial tension in the simulation we use γ φ ≈ 10 pN μ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 0 = 0.02. The fluid viscosity is chosen as η ≈ 1.2 × 10 −3 Pa s. We simulate for approximately 100 time steps with adaptive step size and a total simulation time of approximately 20 ms.
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 × 40 × 40 with some simulation lattices extending up to 800 × 40 × 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 -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 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 0 = 0.002 and increased resolution with a radius of 13 LBM grid cells. Here, a typical fluid lattice consists of 180 × 70 × 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.

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. Re 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. η o = η. 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 I ν (x), K ν (x) being the modified Bessel functions of the first and second kind, respectively, and of order ν. Positive values of ω correspond to growing perturbations (2.1), whereas perturbation modes with negative growth rate are dampened. Because of the positive prefactor ω S 0 = γ φ /(R 0 η), which is the inverse of the viscocapillary time based on the azimuthal tension γ φ , 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 − γ φ /γ z and γ φ /γ z . This range depends on the square root of the anisotropy of the interfacial tension. We recover for isotropic tension γ z = γ φ = γ the range of growing wavelengths between − γ φ /γ z = −1 and γ φ /γ 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 γ 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 ω is plotted only against positive kR 0 due to symmetry. We further distinguish the individual contributions from γ φ , the first term in (3.1), and from γ 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 γ z and γ φ that determines the dispersion relation. The azimuthal tension γ φ , 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 γ 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 γ φ and γ 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 γ φ rises compared to γ z as shown in figure 2(b), the range of growing wavelengths increases. On the other hand, if γ φ decreases relative to γ z , the range becomes smaller (see figure 2c). Because the destabilising γ φ contribution reaches a maximum at kR 0 = 1.59 and tends to zero for kR 0 → ∞, independent of the anisotropy ratio, also for vanishing γ z → 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 γ z = 0 with an extended kR 0 -range on the horizontal axis. In the limit γ φ = 0 all modes are stable. Furthermore, changes in the anisotropy ratio shift the position of the maximum of the dispersion relation. If γ φ > γ z (figure 2b), the position of the maximum of ω shifts to larger values of kR 0 , if γ φ < γ z the maximum is found at smaller kR 0 as shown in figure 2(c). This means for dominating axial tension γ z the instability wavelength increases.

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 FIGURE 2. Dispersion relation in the Stokes regime for η = η o . Curves are shown for (a) isotropic interfacial tension γ z /γ φ = 1.0 and for anisotropic interfacial tension with (b) γ z /γ φ = 0.5 and (c) γ z /γ φ = 2.0. We distinguish the contributions from γ φ (green) and γ 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. γ z = 0. The γ φ contribution (green) has its maximum at kR 0 = 1.59 in each of the panels, because γ φ 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.
Compared to the dispersion relation for a Stokes fluid in (3.1), we here obtain an equation for the squared growth rate ω 2 . According to the ansatz for the perturbed interface in (2.1), perturbations with real and positive ω will grow. Imaginary ω describe oscillatory perturbations of the surface which do not grow in time. Imaginary ω correspond to negative values of ω 2 . Each positive ω 2 has a positive and negative solution ω. The positive solution will grow while the negative is damped. Thus, we are interested in non-negative values of ω 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 γ z /γ φ enters the equation in the same way. However, the prefactor of the growth rate changes ω 2 0 = γ φ /(ρR 3 0 ), 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 γ φ . Furthermore, the geometrical factor containing the Bessel functions is remarkably different.
We distinguish the contributions from γ φ (green) and γ z (orange). While γ z is purely damping, γ φ 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).
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 γ φ contribution. Here for an ideal fluid the destabilising γ φ contribution no longer reaches a maximum at finite kR 0 but instead increases indefinitely. Thus, in the limit γ 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 γ φ /γ z and the maximum shifts towards larger wavenumbers (compare e.g. figure 2b to figure 3b).
In appendix E we further generalise our results to the dispersion relation including a general density and viscosity contrast as derived by Tomotika (1935).

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 λ 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 γ z /γ φ , 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 Oh ≈ 0.00025, i.e. very close to the ideal fluid limit. i.e. the dominant wavelength λ m , changes with the ratio γ z /γ φ . 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 ≈ 0.562 (i.e. λ m /R 0 ≈ 11.18) for isotropic tension γ z /γ φ = 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 γ z /γ φ 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 γ φ -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 λ m = 3.96R 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 γ φ or axial tension γ 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 λ m → ∞ for γ z → ∞ can be understood on the basis of the Young-Laplace equation for anisotropic interfacial tension in (2.4), as well. For infinite γ 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 γ z → 0. We start with considering the ideal fluid. Due to finite γ φ 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 γ z → 0 there is no counteracting contribution which balances this tendency. The monotonic increase of the growth rate for the γ φ -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 λ m → 0 for the ideal fluid. In the Stokes limit, however, the γ φ -contribution is finite for large kR 0 , possibly due to viscous stresses from the fluid which have a damping effect on the perturbation.

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 γ φ = 10 −4 N 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 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 Oh φ or Oh z varies, while the other can be kept constant. Here, we keep Oh φ ≈ 5.5 shown by the dark-green downwards-pointing triangles in figure 5(b). Consequently, Oh z changes from 30 to approximately 2.5, as shown by the light-green upwards-pointing triangles. The transition in the Ohnesorge number 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 Oh → ∞. 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.

Dominant growth rate
We now investigate the growth rate of the most unstable mode ω 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 γ 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 ω m ≈ 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 γ φ . 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.

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 from simulations with initial perturbation amplitude 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 ω 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 γ z /γ φ = 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. 2) (red line) is shown with corresponding BIM simulations (triangles) and LBM/IBM simulations (squares), respectively, depending on the tension anisotropy γ z /γ φ . 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.

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 Guglietta et al. 2020). In combination with typical fluid viscosities ( § 2.4) and sizes of the order of R 0 = 1 μm for vesicles and R 0 = 4 μm for red blood cells this leads to typical viscosity ratios 2η S /(ηR 0 ) = 0.1-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 η S have little effect on dominant wavelength and growth rate, but also around a fraction of 1 the anisotropy dominates. If the interface viscosity η 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.

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. η 0 = 0 and ρ 0 ρ. 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 Ξ = 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, Ξ 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 Ξ 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 ≈ 0.2, the volume assumes its maximum at γ z /γ φ = 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 γ z /γ φ = 0.5 to γ z /γ φ = 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.

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 One simulation corresponds to one wavelength, which is determined by the number of maxima (different curves) per box length. The radius averaged over all maxima divided by the unperturbed radius R 0 is shown over time. The inset shows the growth up to a radius of 110 % of R 0 . The first simulation to reach this threshold is considered as the fastest growing mode. (b) LBM/IBM simulations of a single period of the perturbation with increasing resolution and smaller initial perturbation amplitude 0 = 0.002 allow us to determine the growth rate (here shown for different perturbation wavelengths with γ z /γ φ = 2 in comparison to the analytical solution (3.2)) and the nonlinear correction of the linear break-up time.
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 (t = 0) = 0.02R 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 γ z /γ φ = 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 (t = 0) = 0.02R 0 all modes grow with a different speed. In order to determine the fastest growing wavelength, we define a threshold for the amplitude crit = 0.1R 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 + 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 γ z /γ φ 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 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 Δ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 b and the linear break-up time determines the nonlinear correction as given by (4.1).

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. η o = η. The following calculations are based on the method introduced by Stone & Brenner (1996). We use an ansatz of a cylindrical interface which is slightly perturbed in a periodic fashion with time dependent amplitude (t) 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 Δ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. 1, (B 2) The traction jump can be decomposed into a constant pressure contribution p 0 = γ φ /R 0 and a perturbation in the traction jump which is evaluated at the position of the unperturbed interface due to linearisation (Tomotika 1935;Stone & Brenner 1996). 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 1996) such that the Stokes equation becomes The ring force in the radial direction (third term entering in (B 4) with the radial unit vectorê 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 δ(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 1996) 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 γ φ /η is chosen such thatv r andv z become dimensionless, where the same applies for the pressure prefactor. The velocity prefactor γ φ /η is identical to the viscocapillary velocity based on the azimuthal tension γ φ . The amplitude of the perturbation varies in time, i.e. = (t). For solving the Stokes equation we introduce the Hankel transforms (Poularikas 2000) of the velocity amplitudes V r (s), V z (s) and the pressure P(s) where J ν is the Bessel function of the first kind and order ν.
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 2000): and Together with the Hankel transform of the ring force We can solve this system of equations for the radial velocity in Hankel space and obtain (B 20) 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 ω, i.e.
(t) = 0 e ωt . (B 22) 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 2007, 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 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 ρ o = ρ. 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 γ φ /(ρR 0 ). The prefactors are chosen such thatv r (r),v z (r) are again dimensionless.
In total, after transformation into the Hankel space, we obtain the analytical equations (B 29) 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. 2007, 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 η, a viscosity of the interface η 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 1960;Narsimhan et al. 2015;Sprenger et al. 2020). 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 (2010) 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 onv z (R 0 ). In order to obtainv 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 identifȳ ds, (B 35) which is then solved forv z (R 0 ). Following the same steps as in appendix B.1 we then calculate the growth rate ω: The integral in this equation has already been solved in appendix B.1. Combining (B 36) and the solution forv z (R 0 ) and solving the integral in (B 35) (based on 49(13) and 49 (14) of Erdelyi et al. 1954) leads to the dispersion relation in (5.1).
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 δp = δpI 0 (kr) exp(ωt + ikz), which together with (C 8) leads to the magnitude of the pressure perturbation .
Since 0 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 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 I 0 (x) = 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. . Results for an ideal fluid jet without ambient fluid. The dispersion relation for an ideal fluid jet without ambient fluid with ρ o = 0 is shown for (a) isotropic interfacial tension γ z /γ φ = 1.0 and for anisotropic interfacial tension with (b) γ z /γ φ = 0.5 and (c) γ z /γ φ = 2.0. We distinguish the contributions from γ φ (green) and γ z (orange). Depending on the tension anisotropy (d) wavelength and (e) growth rate (blue curves) are shown in comparison with the results presented above in figures 4(b) and 6(b) for an ideal fluid (red curves). While the wavelength is similar to the case of an ideal fluid, for small and intermediate anisotropy ratios the growth rate is visibly larger.
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 γ z /γ φ = 1 we recover the well-known result (Rayleigh 1878;Drazin & Reid 2004;Eggers & Villermaux 2008) for the dominant wavelength of k m R 0 ≈ 0.697 (i.e. λ m /R 0 ≈ 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).

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 1931;Eggers & Dupont 1994;García & Castellanos 1994;Eggers 1997) for a fluid jet without ambient fluid, i.e. η o = ρ 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ˆ 1 (Eggers & Villermaux 2008). The axial velocity is expanded in radial direction using a Taylor series and considering axisymmetry v z (z, r) = v 0 (z) + v 2 (z)ˆ 2 r 2 + O(r 4 ).
(D 1) From the continuity equation (2.6) one can obtain the radial velocity v r = − 1 2 v 0ˆ r − 1 4 v 2ˆ 3 r 3 + O(r 5 ), (D 2) 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 p(z, r) = p 0 (z) + p 2ˆ 2 r 2 + O(r 4 ).

(D 3)
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 − n · σ · n +p 0 + ηv 0 (D 5) and the tangential force evaluated at the interface r = R is − t · σ · n −η 2Rv 2 − 1 2 Rv 0 − 3v 0 R , (D 6) with t pointing along the interface which is perpendicular to 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ˆ , which gives the dynamics of the interface (D 10) 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 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 Δ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 Comparison of the small-k approximation with the analytical results. The inset shows the dispersion relation obtained by the small-k approximation (orange line) compared to the analytically obtained accurate dispersion relation (C 15) for an ideal fluid jet without ambient fluid (blue line) for a tension anisotropy of γ z /γ φ = 0.5. Systematically varying the tension anisotropy shows that the error (D 19) between exact theory and approximation strongly decreases towards large tension anisotropy. Therefore, the small-k approximation is less accurate for small anisotropy, but becomes a very accurate approximation for large anisotropy.
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 γ z /γ φ = 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 disp = ω 2 − ω 2 ska 2 ω 4 max (D 19) as the error. We observe a strong increase of the error towards small anisotropy. Nevertheless, the deviation for γ z /γ φ = 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 (1935) and introduce a general density N ρ = ρ/ρ o and viscosity contrast N η = η/η o between inner and outer fluid. We use the following abbreviations concerning the wavenumber and the modified Bessel functions , (E 1a-c) 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 1935) 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 ρ = 1 and N η = 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 (1961) 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 2(kR 0 ) 2 ω ρ η R 2 0 (2F(kR 0 ) − 1) + 4(kR 0 ) 4 (F(kR 0 ) − F(y)) + ω 2 ρ 2 η 2 R 4 0 F(kR 0 ) This contains the dispersion relation of an ideal fluid jet in (C 15) in the limit of negligible viscosity.