Nearly spherical vesicles in an external flow

Tank-treading, tumbling and trembling are different types of the vesicle behavior in an external flow. We derive a dynamical equation for nearly spherical vesicles enabling to establish a phase diagram of the system predicting the regimes. The diagram is drawn in terms of two dimensionless parameters depending on the vesicle excess area, fluid viscosities, membrane viscosity and bending modulus, strength of the flow, and ratio of the elongational and rotational components of the flow. The tank-treading to tumbling transition occurs via a saddle-node bifurcation whereas the tank-treading to trembling transition occurs via a Hopf bifurcation. We establish a critical slowing near the merging point of the transition lines.


I. INTRODUCTION
Vesicles, which are closed membranes separating two regions filled up by generally different fluids, have very much in common with biological cells and have a number of applications in pharmaceutics. The truly nonequilibrium problem of describing the vesicle dynamics in external flows attracted lots of attention both from experimentalists [1,2,3,4,5,6,7] and theoreticians [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. In experiments different regimes of vesicle dynamics were observed, usually in share flows. In the tank-treading regime a vesicle shape is stationary. In the tumbling regime a vesicle experiences periodic flipping in the shear plane. Trembling, experimentally discovered in the work [7], is an intermediate regime between tanktreading and tumbling, in the regime a vesicle trembles around the flow direction. Theoretical studies [24] and [20] have predicted similar vesicle behavior and called it vacillating-breathing and swinging respectively. The different types of the motion are illustrated in Fig. 1.
In this paper we focus on constructing the phase diagram, which would predict what kind of regime corresponds to a given set of external parameters, such as vesicle non-sphericity, viscosity contrast, strength of external flow. Although this problem was addressed previously by different authors, it is still far from being closed. As long as no analytic solution of this problem exists, theoretical studies were based either on numerical simulations or on some approximations allowing analytical treatment. Numerical investigations of this problem involved several different computational schemes, including boundary element method [8,9], mesoscopic particlebased approximation [10,11,12,13,14], and an advected field approach [15,16,17,18]. These approaches have shown qualitative agreement with experiment. However, they did not solve the problem of constructing the vesicle phase diagram completely. Analytical studies of the problem can be divided in two major classes. In the first one, see Refs. [11,19,20], phenomenological models of vesicle dynamics based on the classical work of Keller and Skallak [21] were proposed and proved themselves to be rather efficient in explaining particular experiments. However, they cannot cover wide range of parameters and cannot, consequently, be used for constructing a phase diagram for the vesicle dynamics. In the second series of works, Refs. [22,23,24,25], the studies are focused on nearly spherical (quasi-spherical) vesicles whose shape is close to spherical one and can be parameterized by an expansion over spherical harmonics.
We propose a natural extension of the theory developed for quasi-spherical vesicles by accounting higherorder expansion terms. A perturbation scheme around the Lamb solution for spherical body in an external flow [39] allows one to derive the dynamic equations for the shape and the orientation of a vesicle and to investigate them analytically. We show that the high-order in the vesicle distortion terms produce a qualitative change in the phase diagram and make it significantly more complicated. The resulting diagram contains all three types of vesicle behavior which were observed experimentally. In our work, we analyze how the vesicle dynamics depends on different control parameters, such as viscosity contrast, vesicle excess area, internal membrane viscosity, strength of the flow, and ratio of the elongational and rotational components of the flow. We analyze also the vesicle orientation in the tank-treading regime. Some of the results, derived in this paper were already reported in [48]. Here we present significantly more detailed derivations and analyze several aspects of the problem which were not discussed in [48].
Speaking about membranes we have in mind lipid bilayers, the simplest type of the biological membranes. Physical properties of such objects have been extensively studied, both experimentally and theoretically (see e.g. the books [26,27,28] and the reviews [29,30,31]). There are several features of the membranes which are important for our analysis. First, we assume that the membrane is in a fluid (isotropic) state (is a 2d liquid), which is typical of lipid bilayers under normal conditions. Second, we assume that the vesicle has an excess area, that enables one to treat the membrane as incompressible. Third, we assume that the membrane is impermeable to the surrounding fluids, the condition is usually well satisfied in experiment. Finally, we take into account the membrane internal viscosity, which can play an essential role, say, in the vicinity of the lipid-bilayer melting point [32].
The structure of our paper is as follows. In Section II we expose basic theoretical facts concerning physics of membranes. A special attention is paid to dynamic properties of the membranes which can be treated as moving interfaces immersed in a 3d fluid. In Section III we formulate peculiarities of nearly spherical vesicles, analyze their equilibrium properties, and derive phenomenological equation which describe their dynamics in weak external flows. In Section IV we establish a dynamic equation for the vesicle evolution obtained in the framework of expansion over deviations from the ideal spherical shape. We introduce dimensionless parameters controlling the vesicle behavior. In Section V we consider the case of planar external velocity field. In this case the dynamic equations can be essentially simplified enabling one to establish the phase diagram of the system containing domains corresponding to tank-treading, trembling and tumbling. In Section VI we analyze limit cases of weak and strong external flows where it is possible to study the vesicle behavior in details. In Conclusion we discuss some outcomes of our work and its possible extensions. Some technical details of the perturbation expansion are presented in Appendix.

II. BASIC RELATIONS
We are interested in processes which take place at scales of the order of the vesicle size which are assumed to be much larger than the membrane thickness. This assumption is well justified for giant vesicles, usually examined in experiment. In the main approximation, the membrane can be treated as infinitesimally thin, that is as a 2d object (film) immersed into a 3d fluid. In this case the vesicle is characterized by its geometrical shape. In other words, in this limit the vesicles can be considered as active interfaces separating different pieces of the fluid.
The two membrane properties, its incompressibility and impermeability, imply that both the vesicle volume V and its surface area A are conserved, provided the vesicle has an excess area. The latter can be characterized by a dimensionless factor ∆, which is traditionally introduced as where r 0 is a vesicle "radius" determined by its volume. The excess area factor is non-negative, ∆ ≥ 0, and the minimal value ∆ = 0 corresponds to the ideal spherical geometry. The nearly spherical (quasi-spherical) vesicles are characterized by the condition ∆ ≪ 1.
The energy of the membrane is determined by its bending distortions and can be written as the following surface integral [34,35,36,37] over the membrane position. Here κ andκ are bending modules, H and K are the mean and Gaussian curvatures, respectively. They are related to the local curvature radii of the membrane R 1 and R 2 as In accordance with the Gauss-Bonnet theorem, the second term in the right-hand side of Eq. (2) (with the factorκ) is invariant under smooth deformations of the membrane shape. Therefore the term is irrelevant provided the vesicle topology is fixed. Besides the bending energy (2) the membrane is characterized by its surface tension σ. One should stress that for vesicles with excess area the surface tension is an auxiliary variable adjusting to other quantities characterizing the membrane and the surrounding flow to ensure the membrane incompressibility.

A. Flow near the vesicle
In this paper we consider the situation where both fluids, outside and inside the vesicle, are Newtonian. Furthermore, we assume that the Reynolds number associated with the vesicle dynamics is vanishingly small, which is the case in real experiments [1,2,3,4,5,6,7]. Under these assumptions the fluids can be described by the Stokes equation (4) has to be supplemented by the incompressibility condition ∇v = 0, which leads to the Laplace equation ∇ 2 P = 0 for the pressure.
We divide the flow near the vesicle into two parts: an external flow, which would be observed in the fluid in the absence of the vesicle, and an induced flow, which is excited as a result of the vesicle reaction to the external flow. The external flow is assumed to be stationary, its velocity is designated as V . One should remember that the vesicle is advected by the flow, and therefore the above assumption should be valid in the Lagrangian reference frame attached to the vesicle. Below, we neglect the term with the time derivative in Eq. (4) since the characteristic time scale associated with the vesicle dynamics is assumed to be large compared to the viscous relaxation time ̺r 2 0 /η. We assume that the characteristic spatial scale of the external flow is much larger than the vesicle size. In this case the external velocity V near the vesicle can be approximated by a linear profile, determined by a derivative matrix ∂ k V i . The incompressibility condition implies that the matrix ∂ k V i is traceless. Generally, the external flow has two contributions, elongational and rotational: whereŝ is the strain matrix (symmetric part of the matrix ∂ k V i ) and ω is the angular velocity vector. The strain can be characterized by its strength s, defined as s 2 = (1/2)trŝ 2 . Note that for a shear flow s = |ω| =γ/2, whereγ is the shear rate. The fluids inside and outside the vesicle are assumed to be different. We use the designations η for the external fluid viscosity, the viscosity of the internal fluid is designated asη. An important parameter which controls the tank-treading to tumbling transition is the viscosity contrastη/η. The limit where the viscosity contrast tends to infinity corresponds to a solid body behavior of the vesicle, that was concerned by Jeffery [46].
The membrane moves together with the fluid that is the velocity field v is continuous on the membrane and the field v determines the membrane velocity as well as the fluid velocity. For relatively slow processes, we are investigating, the membrane can be treated as locally incompressible, that leads to the condition to be satisfied on the membrane. Here δ ⊥ ik is the projector to the membrane, it can be written as δ ⊥ ik = δ ik − l i l k where l is the unit vector normal to the membrane. The 3d incompressibility condition ∇v = 0 together with Eq. (6) leads to the relation l i l k ∂ i v k = 0, to be satisfied at both sides of the membrane.

B. Membrane stress
The membrane reaction is characterized by its surface stress tensor T (s) ik . There are three contributions to the tensor related to the bending energy (2), to the surface tension of the membrane, and to the internal membrane viscosity: where σ is the surface tension coefficient and ζ is the membrane (2d) dynamic viscosity. Note that the surface tension σ plays an auxiliary role being adjusted to other stresses to ensure the local membrane incompressibility. An expression for the bending contribution to the surface stress tensor was found in the work [42] (see also the book [33]). It can be written as where H = ∇l is the membrane mean curvature and ∂ ⊥ i is defined by Eq. (6).
The surface force f (force per unit area) associated with the membrane stress tensor T (s) ik can be calculated ik . Then one obtains from the expressions (7,8) three contributions to the surface force where Here, again, H and K are the mean curvature and the Gaussian curvature of the membrane, and ∆ ⊥ is the Laplace-Beltrami operator, ∆ ⊥ = ∂ ⊥ i ∂ ⊥ i , associated with the mebrane. Note that the expression for the force (10) can also be derived by calculating the variation of the bending energy (2) due to infinitesimal membrane deformations [38].
The surface force f is compensated by the momentum flux from the surrounding medium to the membrane. This flux consists of two parts, related to the fluid pressure and the fluid viscosity. As a result of the balance, we find the following relations for the normal and tangential to the membrane components of the force. Here, we assumed that the unit vector l is directed outwards the vesicle and the subscripts in and out label regions inside and outside the vesicle, respectively. Thus, P in − P out is the pressure difference between the inner and outer regions that is the pressure jump on the membrane. Note that a fluid viscous contribution is absent in Eq. (13) due to the condition l i l j ∂ i v j = 0, following from the membrane incompressibility (see above). To find the velocity field at a given membrane shape one should solve the stationary Stokes equation η∇ 2 v = ∇P (inside and outside the vesicle) with the boundary conditions (6,13,14) on the membrane. An additional boundary condition reads that v → V far away from the membrane. Note that due to linearity of the equations and the boundary conditions for the velocity a solution of the equations can be written as a sum where v (s) is proportional to the gradient of the external flow (5), and v (κ) is proportional to the bending modulus κ. Of course, v (s) and v (κ) are complicated functions of the vesicle shape.

C. Membrane shape parametrization
Below, we use a particular parametrization of the vesicle shape which is where r 0 is determined by the relation (1). We use the spherical coordinate system with variables r, θ, φ and with the origin in the center of the vesicle. The dimensionless radial displacement u characterizes the deviations of the membrane shape from the spherical one.
The relations (17,18) are formally exact. However, they can be directly used only if u(θ, ϕ) is a single-valued function. Advection of the membrane by the surrounding fluid implies the following kinematic relation Here v r , v θ , v ϕ are spherical components of the velocity v taken at the membrane, that is at r determined by Eq. (16). Again, the relation (19) is formally exact, but can be directly used only if u(θ, ϕ) is a single-valued function. Closed equation which describes the dynamics of the membrane displacement u can be derived in two steps. First, one should find the fluid velocity profile for a given displacement u(θ, ϕ). Second, one should use the kinematic relation (19). This procedure results in a closed non-linear equation for u. Note that due to the property (15) the expression is a sum of two terms, which are proportional to the external flow gradient (5) and to the bending modulus κ. Of course, both terms are non-linear in u.

III. PERTURBATION EXPANSION: WEAK FLOWS
Below, we consider nearly spherical (quasispherical) vesicles, that is the excess area parameter ∆, introduced by Eq. (1), is considered to be small. In this case the dimensionless displacement u is small and it is possible to develop an expansion over u. This perturbation series is a basis for subsequent consideration.
It is natural to represent the function u(θ, ϕ) as a sum over spherical harmonics: The homogeneous contribution to u (its zero angular harmonic u 0,0 ) can be expressed via the inhomogeneous one (non-zero harmonics) from the relation (17) which reflects the volume conservation. Substituting the obtained expression for the zero angular harmonic into Eq. (18) we obtain an expression for ∆ whose expansion over u starts from the second order term. Therefore the displacement u can be estimated as √ ∆. That justifies the expansion over u.
Different angular harmonics in u play different roles. The zero harmonic u 0,0 can be excluded from the beginning, as we explained. First order harmonics u 1,m , describe a shift of the vesicle without changing its shape, and therefore do not play any important role in the vesicle dynamics. The most essential role is played by the second angular harmonic which determines mainly the vesicle shape (at small ∆). As to higher harmonics, they relax fast in comparison with the relatively slow dynamics of the second harmonic. Therefore, the higher harmonics also do not play an essential role in the vesicle dynamics. To avoid a misunderstanding, let us stress that the last assertion is valid for stationary external flows. As it was discovered experimentally [49] and explained theoretically [50], at some conditions (abrupt inversion of the external purely elongational flow) high angular harmonics are generated, the phenomenon is called wrinkling.
The vesicle shape depends on the strength of the external flow. In weak flows it is close to an equilibrium one whereas in strong flows it is determined by the velocity gradient matrix (5). In this section we consider the first case. We discuss the equilibrium vesicle shape and then develop phenomenology for the vesicle dynamics in weak flows.

A. Equilibrium
In the absence of the external flow, an equilibrium shape of the vesicle can be found by minimization of an effective free energy where the first term is determined by the expression (2) and r 2 0 ∆ is the membrane excess area expressed in terms of the displacement u. The Lagrange multiplierσ, related to a fixed value of the membrane area, coincides with the equilibrium value of the surface tension. The second Lagrange multiplier (related to the volume V) is absent in Eq. (21) since we imply that the zero angular harmonic in an expansion of the displacement u is expressed via other ones from the relation (17). Therefore, the volume conservation is automatically satisfied in our scheme.
If ∆ is small, the principal contributions to the energy (2) as well as to the excess area are of the second order in u. It is convenient to write the contributions in terms of the coefficients u l,m of the expansion (20) of u(θ, ϕ) over the angular harmonics: Note, that the first angular harmonic (with l = 1) is absent in the expansions. The reason is that it corresponds to a vesicle shift as a whole, which does not change the energy and the area of the vesicle. As it follows from Eq. (22), the free energy is minimal if only the second angular harmonic is excited. In this caseσ 0 = −6κ/r 2 0 , which is an equilibrium value of the surface tension.
Note that the second order term (22) is degenerate in m. Therefore, in order to determine the vesicle equilibrium shape, one should take into account terms of higher order in the expansion of the effective free energy (21), which violate the degeneracy. In the main approximation, it is possible to keep third order in u terms, and a contribution to u related to the second angular harmonics.
For a subsequent analysis, it is convenient for us to use the following real basis instead of the traditional angular functions Y 2,m . The functions ψ µ are normalized as In terms of the functions (23), the contribution to u related to the second angular harmonic, can be rewritten as follows where u ν are some real coefficients. Expanding the bending energy (2) and the excess area ∆ upto the third order in u and substituting there the expansion (25) one obtains where summation over repeated indices is implied and we designated Components of the object Ξ µνλ are of order unity, they can be found from the definition (28) after substituting the expressions (23).
Minimizing the free energy (26) over u µ and determining the Lagrangian multiplierσ from the condition ∆ = ∆ (3) , one obtains 1 +σ This is the correction related to the third order term in the expansion of the free energy. The minimum of the energy corresponds to a prolate uniaxial ellipsoid. If the principal axis of the ellipsoid is directed along the Z-axis its shape is determined by the expression Substituting the expression (29) into Eq. (26) we find that the coefficient before u µ u µ term in the effective free energy is estimated as κ √ ∆. It contains an extra small factor √ ∆ in comparison with the natural estimation κ. Thus, both, second order and third order, terms in the free energy (26) are of the same order. That gives a formal justification of taking the third order term into account despite the smallness of u.

B. Weak external flow, phenomenology
Here, we analyze the case of weak external flows that cannot significantly distort the vesicle equilibrium shape. As we established in the preceding subsection, the equilibrium shape of a nearly spherical vesicle is the prolate ellipsoid possessing the uniaxial symmetry. An orientation of such ellipsoid in space can be characterized by a unit vector n directed along the principal axis of the ellipsoid. If n is directed along the Z-axis then the vesicle shape is determined by the expression (30). Note that the vectors n and −n describe the same physical state since the ellipsoid is invariant under inversion.
One can formulate a phenomenological equation for the dynamics of n in a weak external flow: where ∂ k V j is the velocity gradient matrix of the external flow and D ijk is some tensor related to the vesicle orientation. Due to the symmetry n → −n the tensor D ijk contains only odd powers of n. Using the relation n 2 = 1, that is n i D ijk = 0, we arrive to the following equation containing a single dimensionless parameter D. Deriving this equation, we exploited the fact, that the vesicle dynamics should be purely rotational ∂ t n = ω × n in a case of an external flow ∂ j V i = −ǫ ijk ω k corresponding to a solid rotation. The factor D in Eq. (31) is dependent on relative significance of the viscous mechanisms and on the excess area parameter ∆, the dependence will be established further, see Section VI. For an external shear flow, it is convenient to utilize the following parametrization of the unit vector n n = (cos ϑ cos φ, cos ϑ sin φ, sin ϑ) .
The components here are written in the Cartesian reference frame attached to the flow: the X-axis is directed along the velocity and the Z-axis is directed opposite to the angular velocity vector ω (see Fig. 2 for clarification). Substituting the expression (32) into Eq. (31) and taking into account that the shear velocity gradient matrix has the only non-zero component ∂ y V x =γ one obtainṡ Note that the dynamics of the angle φ is separated. The equations (33,34) resemble equations for a single polymer dynamics examined in Ref. [44]. The equations (33,34) lead to either tank-treading or tumbling regimes of vesicle motion. For |D| > 1 the tank-treading regime is realized, with a steady tilt angle (between the vector n and the velocity direction) Otherwise, for |D| < 1, the tumbling regime takes place: the vector n experiences a time-periodic motion with an average rotation in the shear plane. Thus, D = 1 corresponds to the tank-treading to tumbling transition. As follows from Eq. (33), the transition is realized via the saddle-node bifurcation.

IV. PERTURBATION EXPANSION: GENERAL DYNAMICS
Here, we start to carry out the program (formulated in Subsection II C) leading to a dynamic equation for the dimensionless displacement u. The program can be realized for nearly spherical vesicles by using a generalization of the Lamb scheme. In accordance with Lamb [39] (see also Ref. [40]), a solution of the stationary Stokes equation can be explicitly expressed via the velocity field taken at a sphere both for the internal and external problems. The Lamb scheme can be directly applied to a spherical solid body immersed into a fluid or to a spherical cavity filled up by a fluid. Then the velocity field is expressed in terms of its surface value. For a nearly spherical vesicle the scheme has to be slightly modified. Namely, one should express the velocity field via its value on the sphere of the radius r 0 . The values can be obtained by an analytical continuation of the internal and of the external velocity fields and are slightly different for the internal and for the external problems. The next idea is to represent the boundary conditions as an expansion over the displacement u which is small parameter for nearly spherical vesicles.
In the zeroth approximation, one can ascribe the membrane velocity directly to the sphere r = r 0 ignoring deviations of the vesicle shape from the sphere. Keeping then the lowest in u terms in all expressions one passes to an equation for the displacement u equivalent to the one discussed in Refs. [24,25]. However, as we demonstrated in Ref. [48], such approximation is not self-consistent. The problem is that it leads to a dynamics sensitive to initial conditions. And one can overcome this sensitivity only by accounting for high order terms in u.

A. Closed equation
Here, we derive an equation for the displacement u in the approximation where the membrane velocity and the boundary conditions (13,14) are related to the sphere r = r 0 . Corrections to the equations associated with the deviations of the vesicle shape from the sphere are small in u. However, for the reasons formulated above, we keep the leading non-linear in u term in the expression for the boundary force (9).
Note, first of all, that the variational derivative of the effective free energy (21) can be represented as Therefore the boundary condition (13) can be rewritten as Here, we divided the surface tension, σ, into a homogeneous,σ, and an inhomogeneous,σ, parts. By definition, the zero angular harmonic is absent inσ. Next, for the sphere r = r 0 the average curvature is H = 2/r 0 and , that explains validity of the expression (36).
To find the inhomogeneous part of the surface tension, σ, one has to use the second boundary condition, (14). Taking the derivative ∂ ⊥ i of the relation (14) and relating a result to the sphere r = r 0 one obtains where σ l and v r,l are contributions to the surface tension and to the radial velocity, respectively, associated with the l-th order angular harmonic. As above, the subscripts in and out are related to the interior and exterior regions of the vesicle. Applying the Lamb scheme to the sphere r = r 0 , one finds for the internal problem Here, we used the condition ∂ r v r (r 0 ) = 0, following from the relation l i l k ∂ i v k = 0 where l i = (sin θ cos ϕ, sin θ sin ϕ, cos θ) is the unit vector perpendicular to the sphere r = r 0 . We substituted also v r (r 0 ) = r 0 ∂ t u, which is the kinematic relation (19) taken in the main approximation in u.
For the external problem, one should separate the external flow velocity V since it does not tend to zero as r → ∞. Writing v = V + w, one obtains from the Lamb scheme where, again, we used the incompressibility condition ∂ r w r (r 0 ) = 0 and the kinematic relation w r = r 0 ∂ t u. However, the expressions (40,41) should be modified for l = 2 because of the external flow. The modified incompressibility condition is ∂ r w r,2 (r 0 ) + s ik l i l k = 0 and the modified kinematic condition is ∂ t u 2 = w r,2 + r 0 s ik l i l k . The modified conditions lead to the expressions Collecting together the relations (36-43) one finds a closed equation for the displacement û whereâ is a dimensionless operator with angular components We included into Eq. (44) a dependence on the rotational part of the external flow, which can be established by an account of the non-linear term in the kinematic relation (19). The Z-axis of our reference frame is implied to be directed opposite to the angular velocity ω of the external flow.
Recall that the quantityσ entering the dynamic equation (44) through Eq. (21) is the surface tension σ averaged over angles. As previously,σ is an auxiliary quantity ensuring the surface conservation law. Let us stress thatσ is a function of time adjusting to the current vesicle shape. Note that the strain and the rotation parts of the external flow are separated: the angular velocity ω extends the time derivative (its effect is equivalent to passing to the rotating reference frame) whereas the strain matrix enters the term s ij l i l j playing a role similar to the free energy derivative. The reason is that the elongational part of the flow leads to some viscous dissipation whereas the solid rotation does not imply any dissipation.
One can further elaborate the equation (44). First of all, one can keep in the effective free energy (21) terms of the second and of the third order in u. Higher order terms in F are small since u ≪ 1. The reason why the third order term should be kept besides the second order one is explained above. Next, the expansion of the term s ij l i l j does not contain any angular harmonics with l > 2, so this term does not push u outside the l = 2 subspace. Therefore the higher order angular harmonics in u die out after a finite time due to the last term in (44). Therefore one can use a reduced equation where u contains only second order angular harmonics. The operatorâ in this case is reduced to a constant a = 16 depending on the viscous ratios. The constant a can be called a generalized viscosity contrast. Note that the limit a → ∞ (where the viscosity of the internal fluid or the membrane viscosity tend to infinity) should correspond to a solid body behavior of the vesicle. One can substitute the expression (26) into the equation (44) and then project the equation on the subspace l = 2. Then the equation (44) is reduced to where the subscripts in (∂ ϕ u) µ and (s ij l i l j ) µ designate projections to the basis (23) calculated in accordance with Eq. (24).

B. Rescaled equation
The factorσ in Eq. (46) should be extracted from the condition 2u µ u µ = ∆, which is the area conservation law written in the main approximation in ∆. Then one obtains where terms of the order su are neglected. In accordance with Eq. (15), the right-hand side of Eq. (47) is a sum of two terms, proportional to s ij and to κ. Note that the term proportional to κ, is of the second order in u (the first order term is absent), that justifies keeping this high-order term in the expansion of the free energy. Deriving the equation (47) we neglected terms of order su. They are much less compared to the term with κ provided s ≪ κ √ ∆/(ηr 3 0 ). This is the applicability condition of the equation. However, below we demonstrate that results obtained from the equation (47) can be extended to stronger flows.
The "vector" S µ in Eq. (48) has an absolute value S, given by Eq. (51), its "direction" is determined by the projections of the object s ij l i l j to the basis (23) that is the "direction" is determined by the structure of the strain matrix s ij . Therefore the two parameters, S and Λ, together with the "direction" of the "vector" S µ completely determine a character of the vesicle dynamics in the external flow.
The quantity τ * is the characteristic time scale of the vesicle relaxation. In comparison with the combination ηr 3 0 /κ, related to the external fluid, the time τ * contains an additional factor a, reflecting the contributions of the internal fluid viscosity and the membrane viscosity into the relaxation, and also the factor ∆ −1/2 . This extra factor reflects the slowness of the second order angular harmonic relaxation related to the degeneracy of the second order free energy in m. Due to the degeneracy the relaxation is determined by the third order term in u in the effective free energy, which contains a smallness ∆ 1/2 in comparison with the energy of higher angular harmonics. Therefore the adiabaticity condition, enabling one to use the stationary Stokes equation, can be written as τ * ≫ ̺r 2 0 /η that is aη 2 r 0 ≫ ρκ √ ∆. The inequality is valid because of the large value of the radius r 0 (in comparison with the molecular length) and the small value of ∆.
The parameter S characterizes the relative strength of the external flow. The expression (51) for S can be explained as follows. The external viscous surface force ηs should be balanced by the surface tensionσ times a variation of the vesicle curvature which is estimated as √ ∆/r 0 . Thereforeσ ∼ ηsr 0 / √ ∆. As we have established in Subsection III A, see Eq. (29), the characteristic equilibrium surface tension can be estimatedσ 0 ∼ κ √ ∆/r 2 0 . Ratio of the surface tension values is therefore estimated as S: S ∼σ/σ 0 . The applicability condition s ≪ κ √ ∆/(ηr 3 0 ) of the equations (47,48) can be rewritten as S ≪ 1/ √ ∆, in terms of S.
The parameter Λ determines the relative strength of the rotational part of the external flow. Note, that ωτ * ∼ SΛ. The condition Λ ∼ 1 corresponds to the angular velocity whose effect is comparable with the effect of the strain, the condition gives ω ∼ s/(a √ ∆). This characteristic angular velocity does not coincide with the characteristic value of s, that stresses again different roles of the rotational and of the elongational parts of the external flow.

C. Very strong external flows
Let us consider the case of very strong external flows, S ≫ 1/ √ ∆, that is s ≫ κ √ ∆/(ηr 3 0 ). As we already explained, the membrane surface tension, which is determined by the balance between the viscous surface force ηs and the product of surface tensionσ and the variation of the vesicle curvature √ ∆/r 0 , is estimated as σ ∼ ηsr 0 / √ ∆. Therefore for the external flows, we con-sider here,σ ≫ κ/r 2 0 . The inequality implies that the leading role in the vesicle reaction to the external flow is played by the surface tension.
Based on the inequalityσ ≫ κ/r 2 0 , one could try to neglect the terms with the module κ in Eq. (46) to obtain an equation independent of the curvature energy (2). That would correspond to omitting the contribution v (κ) in the decomposition (15). However, the resulting equation is incorrect. First, there are additional terms in the equation for u, originating from the regular expansion over u and related to an account of deviations of the vesicle shape from the sphere r = r 0 . The additional terms, which can be estimated as su, are of the same order as kept in Eq. (46). We ignored the terms in the equation (47) for u at s ≪ κ √ ∆/(ηr 3 0 ) since then the term ∼ su is small. Now the additional terms should be taken into account. Second, as we will see in the case of the planar flow, the truncated equation (without κ-terms) cannot be used for describing the vesicle dynamics because of its internal properties.
If the terms with the module κ are neglected in the boundary conditions (13,14) then they appear to be invariant under the transformation v → −v, P → −P, σ → −σ. The stationary Stokes equation η∇ 2 v = ∇P and the boundary condition (6) are also invariant under the transformation. The kinematic relation (19) becomes invariant under the transformation if one adds the rule t → −t. Therefore, in this approximation the backward in time evolution of the displacement u is equivalent to the direct evolution in the external flow with the velocity −V . Since the vesicle dynamics is determined by the matrix (5), the transformation V → −V is equivalent to space inversion. That produces some additional symmetry leading to essential consequences for planar velocity fields.

V. PLANAR VELOCITY FIELD
In this section we discuss in more details the case of a planar external velocity field, when the velocity vector V lies in a plane and is independent of a coordinate normal to the plane. Then the velocity gradient matrix ∂ j V i has nonzero elements only in the plane being a 2d matrix. We choose X and Y axes of our reference system parallel to the plane and assume without the loss of generality that the diagonal elements of the matrix ∂ j V i are zero. Two non-diagonal components of the matrix ∂ j V i completely determine the flow, they can be parameterized in terms of the strain s and the angular velocity ω as ∂ y V x = s+ ω and ∂ x V y = s − ω. In particular, for an external shear flow the only nonzero element of the matrix is ∂ y V x =γ and ω = s =γ/2.
For a planar velocity field, the only nonzero element of S µ in Eq. (48) is S 5 = S. Thus only two parameters in Eq. (48), S and Λ, completely determine the type of the vesicle dynamics. One of the main goals of this paper is to construct a "phase diagram" in the plane S − Λ, which indicates a type of the vesicle motion for a given pair of these parameters.
Very strong flows, characterized by S ≫ 1/ √ ∆ need a special consideration due to the reasons formulated above. It is presented in a separate subsection. Surprisingly, the case can be analyzed in the framework of the same equations (47,48).

A. General consideration
For the planar flow, it is convenient to pass from the variables U µ (49) to another set of variables, "angles" Φ, Θ, Ψ, J, defined as where Θ varies from −π/2 to π/2 and J varies from 0 to π/2. The representation (52) satisfies the normalization condition (49) and, correspondingly, consists of four variables instead of five components U µ . As follows from the equations (48), there is a solution with J = 0. Actually, it is a consequence of general symmetry θ → −θ, which is present in the vesicle equations in the planar flow. Due to this symmetry, there exists a symmetric in θ solution of u, that corresponds to J = 0 as it follows from Eqs. (23,25,52). A special analysis is required to analyze the stability of the solution. Our numerical simulations demonstrated the stability. That is why we mainly limit ourselves to the investigation of this case. At J = 0 the displacement u is Therefore the "angle" Φ characterizes the vesicle orientation in the X − Y plane whereas the "angle" Θ determines the vesicle shape. In the case J = 0 the system of equations (48) is reduced to the following two equations Note that the last, nonlinear in u, summand in (47) produces the term cos(3Θ) in Eq. (54) and no term in Eq. (55). Let us find a region of parameters S and Λ where the equations (54,55) admit stable stationary points, the case corresponds to the tank-treading vesicle motion. Equating to zero the right hand sides of the equations, one finds relations determining a stationary point (for given parameters S and Λ). In order to investigate its stability one should linearize the equations (54,55) near the stationary point to obtain τ * ∂ t (δΘ, δΦ) =B(δΘ, δΦ).
(56)    Fig. 3 the attractors of the system (54,55) are limit cycles. They correspond to either tumbling or trembling behavior. The difference is illustrated in Fig. 5 where the Θ − Φ atlas is plotted. The tumbling regime corresponds to a cycle separating the atlas into two regions each containing a pole Θ = ±π/2, and the trembling regime corresponds to a cycle separating the atlas into two regions one of which does not contain any of the poles. In the tumbling regime the "angle" Φ grows without a limit whereas in the trembling regime it varies in a restricted domain. A transition line from tumbling to trembling, obtained numerically, is depicted by a dashed line in Fig. 3.
The diagram has a complicated structure near the special point S = √ 3, Λ = 2/ √ 3. A vicinity of the point is depicted on Fig. 4, where the regions of coexistence of two different stable points and of a stable point and of a limit cycle are shown. More detailed description is given below.
To avoid a misunderstanding, note that there is an additional region in the S − Λ plane where stationary solutions exist, which are stable in terms of the variables Θ and Φ. However, a stability investigation in the framework of the complete equation (48) shows that these solutions are unstable in the extended space with the variables J and Ψ. Therefore these solutions cannot be realized as the tank-treading motion in real systems.

B. Tank-treading to trembling transition
The tank-treading to trembling transition is determined by the condition tr B = 0. The corresponding curve on the S − Λ plane starts from the above special point S = √ 3, Λ = 2/ √ 3 (the point e 1 on Fig. 4) and goes to the right. This curve, marked as red, is described by the equation where S varies from √ 3 to ∞. Above the red curve, at Λ > Λ * , the stationary point loses its stability via a Hopf bifurcation. Let us establish characteristics of the bifurcation.
Expanding the equations (54,55) near the point (57), one finds the equation for a complex variable z where The variable z is expressed via the deviations of the "angles" from their stationary values as Above the transition line, at Λ > Λ * , the vesicle motion is described by a limit cycle with the radius proportional to √ Λ − Λ * , near the transition curve. This motion corresponds to trembling since the radius is small, and the corresponding limit cycle cannot surround a pole, see the atlas in Fig. 5.
Note the critical dependence of all parameters in Eq. (58) on S − √ 3. Taking into account the critical dependence, one concludes that near the line the amplitudes of the Θ and Φ variations can be estimated as √ Λ − Λ * , without a critical dependence on S − √ 3. A vicinity of the special point e 1 needs an additional analysis since the frequency τ −1 * √ S 2 − 3 of the Hopf bifurcation tends to zero at the point and the approximation leading to the equation (58) is not valid there.

C. Tank-treading to tumbling transition
The transition from tank-treading to tumbling is determined by the condition det B = 0. The transition curve on the S − Λ plane, designated as orange, has a complicated shape, it can be described in a parametric form where the parameter ζ varies from 1/ To be more precise the condition det B = 0 determines a stability boundary of the tank-treading regime. The region of parameters 1/ √ 2 < ζ < ζ 0 , determining the part of the orange curve going from the point e 1 to the point e 0 , corresponds to a transition from one tank-treading regime to another one. That is why the dark green region on Fig. 4 corresponds to coexisting two tank-treading regimes. At passing from the region through the red curve the Hopf bifurcation occurs. However, the bifurcation takes place for one of two possible tank-treading regimes, the other one remains stable. Therefore there exists a region of coexistence of the tank-treading and trembling, the region is marked as fuchsia on Fig. 4. Its left boundary, corresponding to an instability of the trembling motion, is found numerically.
Let us consider now the region of parameters ζ 0 < ζ < √ 3/2, giving the upper part of the orange curve. We are interested in the dynamics of the deviations δΘ, δΦ from the stationary values of the "angles". An analysis shows that there exists a linear combination ξ of the deviations possessing a slow dynamics, that is the characteristic relaxation time of ξ scales as √ δΛ where δΛ = Λ − Λ(S) and Λ(S) is the value of the parameter Λ at the transition curve, determined by the expressions (59). Then, using adiabaticity, it is possible to formulate a closed equation for the parameter ξ which is valid at small δΛ and ξ. The expression (60) is characteristic of a saddle-node bifurcation. The parameters of the saddle-node bifurcation (60) have critical behavior near the boundary points ζ = ζ 0 and ζ = √ 3/2 which can be expressed as where f 1 and f 2 are functions of ζ varying less than by 25% as ζ runs from ζ 0 to √ 3/2. In the normalization where ξ 2 = (δΘ) 2 + (δΦ) 2 one finds f 1 ≈ 4.5 and f 2 ≈ 60 at ζ → ζ 0 . The function F 2 tends to zero as ζ → ζ 0 . Therefore higher order terms in the equation for ξ should be taken into account near the point.
The limit cycle which is a result of the saddle-node bifurcation, destroying the tank-treading regime, could be unstable in its turn. The situation is realized between the points e 0 and e 5 , see Fig. 4. Above the segment e 0 e 4 a final result of the instability is another tank-treading regime continuously continuing to larger S. Above the segment e 4 e 5 it is trembling, characterized by a limit cycle which does not pass through the stationary point. After the point e 5 (to the left from the point) the tanktreading regime is destroyed and a limit cycle passing through the stationary point is formed. Above the segment e 5 e 6 the cycle corresponds to trembling, otherwise it corresponds to tumbling.
Both functions, F 1 and F 2 , tend to zero as ζ → √ 3/2. This is the limit of weak external flows, where the vesicle relaxation rate is proportional to the strength of the flow. Then the right-hand side of the equation should be proportional to S which behaves like S ∝ 3 − 4ζ 2 as a consequence of Eq. (59). That explains the dependence F 1 , F 2 ∝ 3 − 4ζ 2 . However, the limit of weak flows needs an additional analysis since, in accordance with results of Subsection III B, there are two soft degrees of freedom corresponding to rotations of the equilibrium uniaxial ellipsoid. We postpone the analysis to the next section.

VI. LIMIT CASES
We established a general picture of the vesicle dynamics in an external flow which appears to be rich of different types of behavior. Particularly, the phase diagram depicted in Fig. 3 contains different domains and has a complicated structure. The situation is simplified for different limiting cases which can be analyzed in more detail.
For example, the equation (55) leads to the conclusion that in the case Λ → ∞ the vesicle rotates with the angular velocity ω. For external flows with comparable strain and rotational parts vesicle rotation is quite natural, since in accordance with the definition (51) the limit Λ → ∞ is achieved either at a → ∞ or at s → 0. The first case corresponds to a solid body behavior of the vesicle, so one should reproduce the classical Jeffery's result [46], which predicts, that for external flows with ω > √ ∆s the tumbling regime supersedes the tank-treading one. The second case corresponds to a purely rotational external flow, where the fluid rotates as a whole with all inclusions.
Below, we analyze more complicated limit cases.

A. Purely elongational flow
The purely elongational flow is realized provided the angular velocity ω is equal to zero, that is ∂ y V x = ∂ x V y = s. Therefore, in our designations, the elongation is directed along the main diagonal in the X − Y plane.
The condition ω = 0 leads to Λ = 0, in accordance with the definition (51). In this case the system of equations (54,55) has a stable stationary point Φ 0 , Θ 0 , determined by the relations The "angle" Θ 0 monotonically decreases from π/6 to zero as S increases from zero to infinity. The value Φ 0 = π/4 is quite natural since it corresponds to the vesicle orientation along the elongation direction, as is seen from Eq. (53). The stability of the point (62) can be easily established from the linearized equations In the limit S ≫ 1 both "angles" relax to their equilibrium values with the same rate 8 √ 10π s(3 √ 3∆ a) −1 . Recently a wrinkling effect was observed in purely elongational flows at a sudden invertion of the elongation direction, see Ref. [49]. The effect can be explained in the framework of our scheme, the corresponding analysis is presented in Ref. [50].

B. Weak external flows
Let us consider weak external flows characterized by the condition S ≪ 1. We have already discussed the case in Section III B from the phenomenological point of view. We can analyze the case in terms of the "angles" Θ and Φ and then establish a value of the phenomenological constant D, introduced in Eq. (31).
As follows from Eq. (54), for S ≪ 1 the "angle" Θ is close to π/6, which is a stable point of the equation. Substituting the value Θ = π/6 into Eq. (55) one obtains a closed equation for the "angle" Φ If Λ < 2/ √ 3 then the equation (65) has a stationary point which is stable. Otherwise Φ increases unlimited, that corresponds to the tumbling regime. Therefore Λ = 2/ √ 3 is the transition point from tank-treading to tumbling.
If Θ = π/6 then the expression (53) describes a prolate uniaxial ellipsoid with the principal axis directed along the vector (32) with φ = Φ and ϑ = 0. Comparing then the equation (65) with Eq. (33) (obtained for a shear flow with s = ω =γ/2) one finds As it should be, the transition point D = 1 from tanktreading to tumbling corresponds to Λ = 2/ √ 3. Let us stress that the value (67) is independent of the character of the external flow. Therefore the equation (31) with (67) is correct for any external flow V .
Note that in the "solid body" limit a → ∞ (where the viscosity of the internal fluid or the membrane viscosity tend to infinity) the quantity D tends to zero. Account of higher order terms in ∆ gives, that in the solid limit D stops decrease at value order of √ ∆. Diminishing of D leads to a solid rotation of the vesicle in particular case of external shear flow as follows from Eqs. (33,34). The behavior corresponds to the classical result of Jeffery [46], who demonstrated that a solid ellipsoid rotates in an external planar flow, provided ω > √ ∆s.

C. Strong external flows, truncated equations
In the case of strong external flows, where S is large, the last (non-linear in U ) term in the right-hand side of Eq. (48) is small in comparison with the first one. If to omit this last term we pass to a truncated equation. In terms of the variables introduced by Eq. (52) the truncated equation is written as a system of equations homogeneous in S. The system (68) corresponds to the limit case considered by Misbah [24] and Vlahovska and Gracia [25]. In the subsection we examine solutions of the system (68). A relation to observable behavior of vesicles, which is not straightforward, is discussed in the next subsection. The system (68) leads to conservation of two quantities, J and an additional integral Υ, which can be introduced via the relation sin Υ Λ − cos Υ = sin Θ Λ − cos Θ cos(2Φ) .
For definiteness, we choose a root of the equation (69) lying in the domain |Υ| < arccos(1/Λ). Existence of two integrals of motion, J and Υ, implies that a character of an evolution, described by the system (68), depends on initial conditions (determining the values of the integrals). Thanks to existence of two integrals of motion, the system of equations (68) can be completely integrated. For the purpose we introduce a variable It turns out, that for arbitrary initial conditions a solution passes through a point, where U 5 = 0. It is convenient to choose initial time, t = 0, as a moment, when U 5 = 0 and U 4 = cos Υ > 0. Then the initial conditions are ρ = 1 and ∂ t ρ = 0 and one derives from the system (68) the following equation which can be obviously solved explicitly. The parameters Φ and Θ are expressed via the variables ρ and Υ as U 4 = cos Θ cos(2Φ) = Λρ − Λ + cos Υ (72) U 5 = cos Θ sin(2Φ) = (τ * /S)∂ t ρ/ρ.
A solution of the equation (71) behaves differently at Λ < 1 and at Λ > 1. At Λ < 1 the variable ρ trends to infinity as time grows, thus U 5 has its utmost value be equal to √ 1 − Λ 2 , that corresponds to a stable point (tank-treading behavior). At Λ > 1 the variable ρ experiences oscillations. Then one conclude from Eqs. (71-73), that the vesicle evolution is described by a limit cycle, its characteristics are determined by values of the integrals J and Υ. If the angle Φ grows (decreases) unlimited, the vesicle is in the tumbling regime. On the contrary, if Φ is bounded, one deals with the trembling regime. It is convenient to represent vesicle dynamics as a geographic atlas, where Θ and Φ play roles of latitude and longitude, correspondingly, see Fig. 5. The trembling regime corresponds to a closed curve on the Θ-Φ atlas, which does not surround a pole. The tumbling regime corresponds to a curve separating the poles. In Fig. 5 such curves terminate at the boundaries of the atlas having the same Θ on the ends, since points on the right and on the left boundaries of the atlas with the same latitude are physically identical. For the truncated system (68) the tumbling and trembling regimes coexist at any Λ > 1. A choice between the regimes is determined by a value of Υ. If cos Υ > 2Λ/(Λ 2 + 1), then the limit cycle corresponds to tumbling, otherwise the cycle corresponds to trembling. If Υ takes one of its boundary values, that is if cos Υ = 1/Λ, then the limit cycle degenerates into a point Φ = 0, cos Θ = 1/Λ.
Thus, the truncated system of equations (68) has two stationary points corresponding to the tank-treading regime.

D. Strong external flows, slow dynamics
We demonstrated that the truncated system of equations (68) can be completely integrated. However, the system cannot be directly used for an analysis of the vesicle dynamics in the limit of strong shears S ≫ 1. Indeed, the system leads to a dependence of its solution on initial conditions and admits different limit cycles for any Λ > 1. Both these properties contradict obviously to the results obtained in Subsec. V A. The contradiction is resolved if one restores the terms omitted in the truncated system (68) and originating from the last (non-linear in U ) term in the right-hand side of Eq. (48). The restored terms provide a relatively slow evolution of both integrals of motion, J and Υ, which leads to a well defined behavior, independent of the initial conditions.
Our nearest goal is to deduce the equations of motion for J and Υ. We consider the case Λ > 1 where the truncated system of equations (68) leads to limit cycles. Then, analyzing the complete system of equations, one can separate fast motion along the limit cycles and relatively slow evolution of the integrals of motion on scales larger than the cycle period. Note, that a typical time of the fast dynamics is τ * /S whereas a typical time of the slow dynamics is τ * , the large ratio of the times justifies the separation. The equation controlling the slow evolution can be found by averaging over the cycle period of the expressions for the time derivatives of J and Υ obtained from the complete system (48). At this averaging, one can use the fast dynamics described by the system (68). The result can be schematically represented in the form ∂ t Υ =Υ(Υ, J, Λ), ∂ t J =J(Υ, J, Λ). Explicit expressions forΥ andJ are quite cumbersome, so we do not present here its final form.
One can check, that the ultimate value of J is equal to zero at any Λ. That is why below we consider the case J = 0, and find the slow dynamics for the quantity Υ at the condition. ThenΥ is a function of Λ and Υ, the function can be found from Eqs. (54,55,69) where angular brackets mean averaging over the cycle period. The evolution of Υ described by the equation ∂ t Υ =Υ leads to a stationary value which can be found from the equationΥ = 0. The value corresponds to a limit cycle, which is stable provided ∂Υ/∂Υ < 0.
In accordance with the analysis, made in the previous subsection, a character of the limit cycle depends on the value of Υ. If cos Υ > 2Λ/(Λ 2 + 1), then the limit cycle corresponds to tumbling, otherwise the cycle corresponds to trembling. A numerical investigation based on Eq. (75) shows that the boundary value of Υ, Υ = arccos[2Λ/(Λ 2 + 1)], is achieved at Λ ≈ 1.52, larger Λ correspond to tumbling, smaller Λ correspond to trembling. At Λ < √ 2 all the limit cycles appear to be unstable. That means that the condition the stationary point (74) appears to be an attractor. Thus, the system falls to the stationary point that corresponds to the tank-treading regime. Note that the position of the stationary point is slightly shifted due to the presence of the additional term in the equation for Θ (54). The new position is in the main approximation in 1/S. A correction to the value Θ = arccos(1/Λ) is of the order 1/S 2 .

E. Extremely strong flows
The above analysis is, strictly speaking, correct for flows characterized by S < 1/ √ ∆. For stronger flows, 1/ √ ∆ S, our consideration should be extended. Some additional terms of the higher order in u should be taken into account in the equation for u, which for 1/ √ ∆ S are larger than those kept in Eq. (47). Leading terms of such kind can be estimated as su, they originate, say, from account of deviations of the vesicle shape from a spherical one. The terms are associated with the contribution v (s) to the membrane velocity, whereas the second term in the right hand side of Eq. (47) is associated with the contribution v (κ) to the membrane velocity, see Eq. (15).
However, there is an essential difference between the terms v (s) and v (κ) . If one omits the term v (κ) then there appears a symmetry of the vesicle dynamics expressed in terms of the equation for u, which is invariant under simultaneous time and space invertions, see Subsection IV C. Since the axes of our reference system are attached to the eigen vectors of the strain matrixŝ, the space inversion is equivalent to the transformation ϕ → −ϕ, θ → θ, that is it can be written as Φ → −Φ and Ψ → −Ψ in terms of the "angles" (52). Therefore the equations for Θ and Φ should be invariant under the transformation One can easily check that the truncated equations (68) are invariant under the transformation (77). However, our analysis demonstrated that the symmetry survives even though higher in u terms will be taken into account provided κ → 0. If small in u corrections to the truncated equation (68) will be taken into account then the system of limit cycles, characteristic of this equation at Λ > 1, will survive with small perturbations. Next, if the term v (κ) in the expression (15) is neglected all the limit cycles will be neutral (they are no stable, no unstable). Indeed, the cycle passing through a point Φ = 0, Θ cannot be stable, since then due to the symmetry (77) it should remain stable after the time inversion. Thus, the corrections to the truncated equation (68) keep its main property which is an existence of an additional integral of motion (leading to the neutrality of the limit cycles).
Thus only the term v (κ) proportional to κ can destroy the integrability corresponding to the neutral limit cycles. It produces a selection leading to stable limit cycles or to stable stationary points. Since the selection is produced among limit cycles slightly disturbed in comparison with the ones corresponding to the truncated equation (68), the results will be the same in the main approximation in ∆. Thus the above results obtained for strong flows can be immediately extended to the extremely strong flows, where S > 1/ √ ∆.

VII. CONCLUSION
We have investigated dynamics of nearly spherical vesicles in an external stationary flow, being mainly focused on the shear flow. The general calculational scheme developed for nearly spherical vesicles enabled us to analyze the dynamics in detail. The scheme is based on solving the 3d hydrodynamic (Stokes) equations with boundary conditions posed on the membrane. Besides the membrane bending elasticity we have taken into account the internal membrane viscosity, leading to an additional dissipation mechanism.
There are essentially different regimes of the vesicle dynamics dependent on a relation between the strain of the external flow and the vesicle relaxation rate. In the weak flows the vesicle shape is close to an equilibrium one which is a prolate ellipsoid. Then a role of the external flow is reduced mainly to an orientation of the ellipsoid. In the strong flows the vesicle shape and orientation are determined by the flow. We established the vesicle behavior for different strengthes of the external flow including its orientation relative to the external velocity.
The most interesting phenomenon in the vesicle dynamics is transition to tumbling regime which occurs at increasing the generalized viscous contrast (45) or at increasing the rotational component of the external flow. We demonstrated that in weak flows a direct transition from tank-treading to tumbling occurs whereas in strong flows an intermediate regime, trembling, is realized. The behavior is in accordance with experiment of Kantsler and Steinberg [7] and corresponds qualitatively to numerics of Noguchi and Gompper [20]. The phase diagram we obtained is plotted in Fig. 3 where the variables S and Λ are some combinations of observable quantities (51).
The possibility of an intermediate regime between tank-treading and tumbling was discussed theoretically, first by Misbah [24] and then (qualitatively) by Noguchi and Gompper [20]. Note, however, that the calculational scheme used by Misbah [24] and then by Vlahovska and Gracia [25] is not self-consistent though formally the authors have taken into account principal terms of the equation for the vesicle distortions. The situation cannot be improved by introducing higher order terms related to the external flow. And only after introducing higher order terms related to the membrane bending elasticity into the equation it enables one to find, say, boundaries of the region where the trembling regime is realized.
The transitions tank-treading to tumbling and tanktreading to trembling have essentially different scenarios. The first transition is described as a saddle-node bifurcation whereas the second one is described as a Hopf bifurcation. Our theory predicts existence of a special point on the diagram where the two above transition lines merge. One expects an essential "critical" slowness of the vesicle dynamics near the point. Therefore it plays a role analogous to some extent to the critical point on fluid phase diagram.
Near the transitions thermal fluctuations are relevant which smear the transitions. One expects that the effect is especially strong near the special point on the phase diagram. Role of the thermal fluctuations, which can be examined in the spirit of the works [43,44,45], constitutes a subject of special investigation, to be done separately. Here we note only, that due to fluctuations one should be careful in comparison an experiment with the phase diagram obtained in our work since the diagram is deduced ignoring the fluctuations.
It is worth to mention the effect, related to the thermal fluctuations, recently discovered by Kantsler et. al. [49]. It was shown that the relaxational dynamics of a vesicle in external elongational flow is accompanied by the formation of wrinkles on a membrane. Theoretical investigation of this effect presented in [50] was based on the theory, developed in this paper. It was shown that the formation of wrinkles is related to the dynamical instability induced by negative surface tension of the membrane.
We investigated nearly spherical vesicles assuming that the excess area factor ∆ is small. It enables one to formulate a powerful calculational scheme enabling to find details of the vesicle dynamics analytically. We have doubts that the scheme can be generalized for a general case ∆ ∼ 1. Probably, the case can be investigated only numerically. However, the general approach should be the same: one has to solve the stationary Stokes equation inside and outside the vesicle with the boundary conditions (6,13,14) at a given vesicle shape and then to use the equation (19) to formulate an equation for membrane distortions.
1). Also we denotê We need in pressure, velocity and some components of symmetric part of velocity gradient S ij on sphere with radius r 0 . Here we list all quantities on inner side of the sphere, omitting index in. Quantities on outer side of the sphere can be obtaining by changing l → −l − 1 and plugging functions (A1) with index out. Pressure inside is where The velocity is given by Tangential-radial components of velocity strainŜ is In (A16) necessary for us coefficientsK are K sx rτ (l) = l − 1 l(l + 1) , K sy rτ (l) = 2l + 1 l(l + 1) .
Tangential-tangential components ofŜ are In the interesting case of Y = Z = 0 one obtains from (A17) and h αβ S αβ s = 0.
The radial derivative is given by ∂ r S rr = 2∂ 2 r v r =K sx+ rrr X + . . . , where we kept the only necessary contribution from the unction X with K sx rrr (l) = 2(1 − l 2 ).
In this Subsection we return to the derivation of Eq. (44). We have seven unknown scalar functions, six in (A1) and surface tension σ, which depend on spherical angles θ, ϕ. To find these functions, one should satisfy all the boundary conditions on the membrane. Three boundary conditions come from velocity continuity in the whole space. For us it is convenient to use three equivalent conditions: continuity of normal component of velocity l i v i , continuity of normal derivative of normal component of velocity l i l k ∂ k v i and normal component of vorticity l k ω k . Another three conditions come from the continuity of the momentum flux, (13,14). Seventh condition corresponds to the surface flow incompressibility, see (6). After the seven scalar filed are found, one should use the relation between the velocity and the temporal derivative of function u (16) presented at (19) and obtain required equation of motion. In this section we put κ = 0, and do not take into account the terms arising from the bending force of the membrane.
To obtain the dynamical equation with the accuracy up to n-th order of √ ∆, one should satisfy all boundary condition with the same accuracy. It is convenient to find the next in √ ∆ correction to the equation using the recursive procedure. In accordance with the scheme we represent any quantity, for example X in , as a series where X in n ∝ ∆ n/2 . For vector quantities which have low indices we put the index n on the left from the main letter to avoid mixing of different indices.
We represent the local value of the surface tension as σ =σ +σ, whereσ does not depend on angles θ, ϕ and the integral ofσ over the angles is zero. The reason of the division is different scaling laws of partsσ andσ with ∆: expansion of the quantities in series over ∆ arē where lower index n near a term corresponds to the term scaling as ∆ n/2 . Let us find equation of motion in zero order in √ ∆. Surface incompressibility, i.e. continuity of l i v i reads Continuity of l i l k ∂ k v i leads to In this subsection we find the first corrections in √ ∆ to the equation of motion (A35). One should repeat the steps (A24,A26,A28, A29,A32,A33), accounting next order in √ ∆ and find correctionu 1 order of √ ∆ tou (A35). In the main approximation we required (A24), that l i v i out 0 s = l i v i in 0 s on the sphere with radius r 0 . Here low index "0" stands for contribution from X 0 , . . ., see definition (A3), whereas low index "1" in (A38) corresponds to contribution from X 1 , . . .. Quantity l i v i calculated on the membrane surface through found coefficients X 0 , . . ., differs from that on the sphere: where indices "s" and "m" stands for sphere surface and membrane surface correspondingly. Condition can be rewritten as Analog of (A29) now iṡ In (A41) velocity whereu 0 is from (A35). Corrections to boundary conditions (A32,A33) arē ∇ ⊥ ασ 1 = η l i S iα 1 s + η δ 1 l i S iα 0 m .
In (A43) In (A44) It follows from evolution equation, written in the main approximation over √ ∆, that at long times all excess area becomes confined in harmonics l = 2. Coupling of higher order harmonics with sector l = 2 occurs only due to nonlinear terms, concerned in Subsection A 4. Relative part of excess area, confined in higher order harmonics is order of ∆. Hence, back influence of higher order harmonics on sector l = 2 appears, if one accounts terms order of s∆, that is the influence is negligible in our approximation. Thus, to obtain corrections to dynamics in sector l = 2 in the approximation, it is sufficient to account only coupling of mode l = 2 with itself.