A theoretical investigation of the frisbee motion of red blood cells in shear flow

The dynamics of a single red blood cell in shear flow is a fluid–structure interaction problem that yields a tremendous richness of behaviors, as a function of the parameters of the problem. A low shear rates, the deformations of the red blood cell remain small and low-order models have been developed, predicting the orientation of the cell and the membrane circulation along time. They reproduce the dynamics observed in experiments and in simulations, but they do not simplify the problem enough to enable simple interpretations of the phenomena. In a process of exploring the red blood cell dynamics at low shear rates, an existing model constituted of 5 nonlinear ordinary differential equations is rewritten using quaternions to parametrize the rotations of the red blood cell. Techniques from algebraic geometry are then used to determine the steady-state solutions of the problems. These solutions are relevant to a particular regime where the red blood cell reaches a constant inclination angle, with its membrane rotating around it, and referred to as frisbee motion. Comparing the numerical solutions of the model to the steady-state solutions allows a better understanding of the transition between the most emblematic motions of red blood cells, flipping and tank-treading. Mathematics Subject Classification. 14-04, 74F10, 74H10, 74K15, 92C05, 92C17. Received July 25, 2019. Accepted February 22, 2021.

Investigating how blood flows is of prime interest for many medical and industrial applications, but the nature of blood makes this task particularly challenging. It has long been recognized that the behavior of an isolated RBC under flow is the foundation stone of the knowledge of hemodynamics. Physiological flows being shear-dominated, the case of an isolated RBC in shear flow has been particularly investigated over the last 50 years, in experiments [1, 4, 10, 14-16, 18, 21, 26, 33, 40] and more recently in numerical simulations of the RBC dynamics [8, 9, 12, 21-23, 28, 29, 36, 38, 45]. When subjected to shear flow, an RBC adapts to the flow conditions then experiences different long-term motions, which are the topic of this paper: the transient response is not investigated here. In the 1960s-1970s, it has been discovered that RBCs may behave similarly to rigid disks [18], flipping in the flow without deformation, or to droplets [33], the RBC reaching a fixed orientation with respect to the flow direction, with its membrane circulating around the cytoplasm, a movement referred to as tank-treading [16] by analogy with the movement of the tread of a tank. The dynamics depends on two parameters which are varied in experiments: the shear rate and the viscosity of the external fluid. More precisely, the non-dimensional capillary number, allows to quantify the ability of a shear flow to deform the RBC: Ca = µ extγ a 0 /G s , where µ ext is the dynamic viscosity of the external fluid,γ the shear rate, a 0 the characteristic size of the RBC and G s the surface shear modulus. In addition, the viscosity ratio λ between the internal viscosity µ int and the external viscosity λ = µ int /µ ext also controls the dynamics [21][22][23]. Recently, much progress has been made in the knowledge of the RBC dynamics under pure shear flow, for a wide range of Ca and λ [1,10,[21][22][23], see [24] for a recent review of these motions. Here, we focus on the cases where the shear rate is small, ie when Ca is small.
When Ca is typically lower than 0.05-0.1, experiments show that RBCs barely deform under shear flow [1,10]. However, they exhibit a rich dynamics that is still not well understood. At very low shear rates, an RBC flips like a solid disk [18]. Jeffery [19] showed that at zero Reynolds number, solid disks rotate in shear flow over orbits that are purely determined by their initial orientation. However, contrary to rigid disks, RBCs flip over orbits that shrink towards the vorticity axis, when increasing the shear rate. This is referred to as the orbital drift of RBCs [7,10,23]. At the end of this drift, the RBC has its small axis aligned with the vorticity axis and it spins at constant angular speed. This motion is called rolling. For higherγ, the dynamics depends on the viscosity ratio λ. At large viscosity ratio (for instance for an RBC in plasma), rolling is stable when the shear rate increases. This is not the case when the viscosity ratio is low: rolling is not stable and the RBC transitions to tank-treading, albeit with orientation oscillations (referred to as swinging), which decrease with the shear rate [1]. The transition between the motion of flipping over orbits and swinging is not understood. Some experiments showed hysteresis in this region, the behavior of the cell at a given shear rate being different if the experiment was performed at increasing or decreasing shear rate [10,26]. In that region, simulations predict a motion called frisbee, hovering or inclined rolling depending on the authors, where the RBC does not lie in the shear plane nor along the vorticity axis [8,23,36]. Examples of these dynamics computed by numerical simulations solving the Navier-Stokes equations coupled to a mechanical model of the RBC membrane [22,23,25,35] are provided in Figure 1a. Among those dynamics, the orbital drift of RBCs and the rolling motion [4,7,10,23,44], and the tank-treading/swinging [1,14,16,41,45] have been the object of dedicated studies. This is not case for the frisbee motion, which has only been mentioned in a few simulation papers [8,23,36].
The complexity of the RBC dynamics has motivated many experimental and numerical investigations, but theoretical works with simplified models have tremendously contributed in enlightening the dynamics of the RBC in shear flow [1,11,20,23,37]. Such models use the fact that at low shear rate, the RBC deformations remain small: they model the RBC as a fixed-shape ellipsoid with a membrane allowed to circulate around it. Note that other models managed to relax the fixed-shape assumption: Vlahovska et al. [42], for instance, predict the dynamics of quasi-spherical capsules with membrane circulation and shape deformations. In the fixed-shape models, the RBC is represented as an ellipsoid that encloses viscous fluid inside its membrane, and the membrane viscoelasticity is accounted for. The recent 3-D version of such a model has retrieved the dynamics of an RBC described in the former paragraph [23], as illustrated in Figure 1b. Figure 1b shows that in spite of the little interest that the frisbee motion has received, it plays a central role in the transition between other motions: it is the dynamics that allows transition between a dynamics where the symmetry axis of the cell is in the shear plane (swinging) to off-plane dynamics as flipping over orbits or rolling. One objective of this  [22,23,25,35] showing examples of Orbit, Frisbee, Rolling and Swinging dynamics (seen from the vorticity axis), with Lagrangian beads to track the membrane displacement: at rest, the white bead is on the dimple and the green and the blue beads are on the rim. The arrows on the left display the shear flow, in the shear plane formed by the flow direction and the velocity gradient direction. (b) Result of the model by Mendez and Abkarian [23] for an ellipsoid of semi-axes a 1 = a 2 = 4.2375 µm and a 3 = 1.2511 µm. Phase diagram as a function of the in-plane capillary number Ca * and the effective viscosity ratio λ ef f (the definitions of Ca * and λ ef f are provided in Sect. 2). The additional hatched region of the phase diagram corresponds to a range of (Ca * , λ ef f ) for which swinging may be obtained if the initial orientation is close enough to the shear plane. It is a zone where hysteresis occurs if the phase diagram is explored at fixed λ ef f and increasing or decreasing Ca * . study is actually to clarify these transitions by further analyzing the frisbee motions. In particular, we aim at determining under which conditions frisbee is possible and how the characteristics of the frisbee motion depend on the parameters of the problem.
The existing 3-D model [23] consists of a system of nonlinear coupled ordinary differential equations (ODEs) describing the rotation of the RBC and the circulation of its membrane along time. However, the model remains complex and further analysis is needed to better understand the transitions between the different motions. Following the work by Jeffery [19], Mendez and Abkarian [23] have parametrized the rotations of the particle with the Euler angles. However, Euler angles are not well suited for mathematical analysis, being notably limited by the gimbal lock problem, when the first and the third Euler angles become indistinguishable if the second Euler angle has some critical values. In the present paper, we show the limitations of the existing system based on the Euler angles and develop a system based on quaternions, to avoid any singularities. Moreover, writing the system in terms of quaternions yield polynomial equations. One may then use algebraic geometry tools to investigate the system. In this paper, in addition of providing the system itself, we use algebraic geometry techniques (elimination theory using the Gröbner bases) to determine which ranges of parameters the frisbee dynamics may be found, and compare them to the numerical results found by solving the system numerically. This notably evidences regions where frisbee solutions exist but are unstable.
In Section 2, the existing model of the RBC dynamics under pure shear flow [23] is first described and the corresponding notations are introduced. The problems relative to the use of the Euler angles will be highlighted. Section 3 develops the new version of the system, based on the parametrization of the RBC rotation by quaternions. Section 4 presents one example of the analysis allowed by the new system, with the explicit determination of the frisbee solutions in the phase diagram.
2. An analytical model for the dynamics of an RBC at low shear rates Mendez and Abkarian (referred to as MA) [23] recently proposed a model to describe the dynamics of an RBC at low shear rates (small values of Ca). This section describes the existing model, introducing the necessary notions and notations. As the RBC dynamics occurs without major deformation when Ca is small, the RBC is represented by an axisymmetric ellipsoid of prescribed shape, which never changes under flow. The choice of an ellipsoid to model the biconcave shape of an RBC is motivated by the result from Bretherton [5], who showed how rigid particles with symmetry of revolution behave as axisymmetric ellipsoids, whose dynamics was determined theoretically by Jeffery [19]. Indeed, Goldsmith and Marlow have shown that the flipping motion of RBCs in their experiment at low shear rates is well reproduced by the Jeffery model for an oblate ellipsoid of semi-axes 4, 4 and 1.5 microns [18]. While the shape of the object is fixed, it can freely rotate around its center of mass. In addition, the membrane of the RBC is represented as a layer that has exactly the shape of the ellipsoid, but that can circulate/slide around the fixed shape. The movement of the membrane also entrains the internal fluid in a consistent manner. The membrane being elastic, its circulation around the shape is associated with the in-plane deformation of the membrane elements and the existence of elastic stress in the membrane [1,2]. The degrees of freedom of the problem are of two types: those describing the orientation of the ellipsoid in space and those describing the circulation of the membrane relative to the ellipsoidal shape. The ellipsoid is subjected to pure shear flow in an infinite domain. When Ca is small, the Reynolds number for such a flow is also small, so that inertia may be neglected. The flow inside and outside the RBC is thus governed by the Stokes equations.

Fixed frame, main flow and geometry
The three-dimensional Euclidian space, V , is endowed with a fixed normal frame (O, e 1 , e 2 , e 3 ). This origin and basis correspond to the laboratory frame that never changes along time, referred to as the fixed frame. Vectors in V will be noted with bold characters, whereas coordinate vectors in the fixed basis ( e 1 , e 2 , e 3 ) will be noted with hats.
The external flow (unperturbed flow) is a pure shear flow defined as u 1 =γ x 2 , whereγ, the shear rate, has the dimension of a frequency. e 1 is thus the direction of the flow, e 2 the direction of the velocity gradient and e 3 the vorticity direction.
In the fixed frame the velocity gradient tensor is: The strain-rate tensor and the skew-symmetric part of the velocity gradient tensor read: The orientation for the fluid ellipsoid modeling the RBC is defined with respect to a reference orientation, called configuration at rest. In the position at rest, the surface of the ellipsoid is defined as 3)). The (x 1 ,x 2 ) plane is the shear plane. (b) The instantaneous membrane tank-treading rate is defined in the body frame by the tank-treading rate vectorΩ. A Lagrangian marker P located at x 1 = x 2 = 0 and x 3 = a 3 at rest allows to track how much the membrane has circulated.
The semi-axes lengths of this ellipsoid are (a 1 , a 2 , a 3 ), with a 1 = a 2 > a 3 , such that E rest is an oblate axisymmetric ellipsoid, centered at the origin, with a small canonical axis (symmetry axis) along e 3 and two long axes along e 1 and e 2 .

Orientation of the ellipsoid during the dynamics
At any instant t, the position of the ellipsoid is defined by rotating the position at rest: where r t is a rotation centered at the origin. The vectors (e 1 , e 2 , e 3 ) = (r t ( e 1 ), r t ( e 2 ), r t ( e 3 )) form the basis of a time-dependent frame, called the body frame. Coordinate vectors in the body frame will be noted with normal characters, neither bolded nor hatted.
In MA [23], rotations of the ellipsoid are parametrized with Euler angles (shown in Fig. 2a), that is three angles (θ, ϕ, ψ) and their corresponding rotations of center O: -r e3,θ , the rotation of axis e 3 and angle θ, -r N,ϕ , the rotation of axis N = r e3,θ ( e 1 ) and angle ϕ, -r e3,ψ , the rotation of axis e 3 and angle ψ (e 3 = r N,ϕ ( e 3 ) = r N,ϕ • r e3,θ ( e 3 )), respectively represented by their matrices in the local frames: The rotation r t thus reads Table 1. Characteristic behavior of ϕ depending on the motion of the fluid ellipsoid (after a possible transient response). In the Orbit motion, ϕ oscillates periodically between ϕ min and ϕ max . In the other motions, ϕ is constant.

Dynamics
Orbit Frisbee Swinging Rolling Behavior of ϕ ϕ ∈ [ϕ min , ϕ max ] constant ϕ ∈ ]0, π/2[ ∪ ]π/2, π[ ϕ = π/2 ϕ = 0 or π and the associated matrix in the fixed frame (and in the body frame) can be expressed as: Any rotation can be written in this way, with ϕ ∈ [0, π] and θ, ψ ∈] − π, π]. The angle ϕ is the angle between the symmetry axis of the ellipsoid e 3 and the vorticity axis e 3 . The cases ϕ = 0 and ϕ = π have to be considered carefully since there is an infinite set of Euler angles representing the rotation: if ϕ = 0 or ϕ = π, the two Euler angles θ and ψ are indistinguishable and only their combination (sum or difference) matters.
The angle ϕ is particularly interesting to describe the dynamics of the fluid ellipsoid, as the characteristic movements at low Ca can be defined from ϕ, as shown in Table 1.

Circulation of the membrane
As already explained, the membrane is able to move along the surface E t . One may define the position of the membrane elements at rest as a reference position for membrane deformation. The membrane elements are then tracked in a Lagrangian way, by defining x(x rest , t), the position of the membrane elements x at time t which is initially at x rest at rest. When the membrane is displaced along the surface, the elastic energy stored in the membrane depends on how much the membrane has circulated. Consider for example the membrane elements initially at the small axis of the ellipsoid x rest = t 0 0 a 3 . If this element remains at the same position, in the body frame, x t 0 0 a 3 , t = t 0 0 a 3 : the membrane is in its equilibrium position and the elastic energy is minimum. On the contrary, if this element is displaced to the rim of the ellipsoid, the elastic energy in the membrane is maximum. This will modify the movement of the ellipsoid.
Following Keller and Skalak [20], MA impose the form of the circulation of the membrane using the following expression for the membrane velocity in the body frame as [23]: where the Einstein summation convention is used for repeated indices. ijk is the Levi-Civita symbol in three dimensions: ijk = 1 (resp. −1) if (i,j,k) is an even (resp. odd) permutation of (1,2,3), and 0 if any index is repeated.Ω j is the tank-treading rate around axis e j . x j are the coordinates in the basis (e 1 , e 2 , e 3 ) of the body frame (see Fig. 2b). Any point on the ellipsoid surface thus remains on this surface, as the movement is only tangential. It is useful to map this movement on a sphere by manipulating normalized coordinates. We define X i = x i /a i . The non-dimensional velocity of any point of the membrane is written as the tank-treading rate being now interpreted as the angular speed of the membrane non-dimensional coordinates. Note that this membrane motion is rather simplistic and does not respect the area-conserving character that the membrane velocity field should have, as the membrane of an RBC is almost inextensible. Secomb and Skalak [34] discussed this issue and proposed more sophisticated membrane velocity fields to respect this constraint, but the complexity of the flow field has not allowed to obtain quantitative results as for the motion of equation (2.5) [40].
The particular case of the axisymmetric ellipsoid with an axisymmetric membrane can be further simplified. Any circulation around the symmetry axis cannot be distinguished from a pure rotation around the symmetry axis:ω 3 andΩ 3 have exactly the same role for ellipsoids of revolution axis e 3 . As a consequence, MA have assumed thatΩ 3 = 0: the motion around the symmetry axis is viewed as a solid-body rotation, not as a membrane circulation.
In addition, the effect of the membrane on the movement is actually only determined by the position of one membrane point initially located on the small axis of the membrane: x 1 = x 2 = 0 and x 3 = a 3 . This point, denoted by P , referred to as the marker, follows the membrane movement, and may be described through its non-dimensional coordinates along time, ξ t . At rest, ξ 0 = t (0, 0, 1).
For an axisymmetric ellipsoid, the MA model for the dynamics of an RBC in shear flow is entirely described with two variable states: where SO 3 (R) is the rotation matrices group of R 3 , and S 2 in the unit sphere of R 3 . MA's model of the RBC dynamics predicting R t and ξ t along time, as a function of the parameters of the problem, is described in the remainder of the section.

Parameters of the problem
Two kinds of parameters are involved in the model: geometrical and fluid/flow parameters. The first parameters we introduce are geometrical and depend on the semi-axes lengths of the ellipsoid (recall that a 1 = a 2 ): .
(2.6) k 1 and k 2 are two geometric ratios characterizing the ellipsoid. The ratio of viscosity between the internal and the external fluid λ = µ int /µ ext is another parameter controlling the dynamics [21][22][23]. In the model, it contributes to the dynamics through a non-dimensional number also function of the geometry, referred to as Λ: The definition of the geometric factors f 1 , f 2 and f 3 is provided in Appendix A. Note that the definition of Λ is based on an effective viscosity ratio λ ef f which accounts for the membrane viscosity (see Appendix A). In the absence of membrane viscosity, λ ef f = λ. Finally, the dynamics depends on the shear rate and more precisely on the competition between the external shear stress imposed, µ extγ , and the in-plane elasticity of the membrane. As shown by Mendez and Abkarian [23], this competition is controlled by a specific capillary number denoted by Ca * and different from Ca: while Ca determines the ability of the flow to change the shape of the RBC, Ca * determines the ability of the flow to make the membrane circulate [23]. The detailed definition of Ca * with respect to the geometry of the problem is provided in Appendix A.
In the remainder of the paper, we fix a 1 = a 2 and a 3 , so that k 1 , k 2 , f 1 , f 2 and f 3 are fixed. The dynamics as a function of two parameters, λ ef f and Ca * is discussed.

Model for the RBC dynamics
The model itself is now described. The aim is not to derive the full model, but to provide the reader with the equations and some indications of how they were obtained. The reader is referred to [2,20,23] for details about the calculations. In the Euler angle formulation, the orientation of the ellipsoid through the rotation r t is determined by the three Euler angles (θ, ϕ, ψ). First, it is easy to show that the derivatives of the Euler angles are related to the angular speed of the ellipsoid around its axes expressed in the body frame and denoted by ω i : ω 1 =θ sin ϕ sin ψ +φ cos ψ,ω 2 =θ sin ϕ cos ψ −φ sin ψ,ω 3 =θ cos ϕ +ψ. (2.8) Note also that equation (2.5) may be applied to the marker P , so that: Then, as the total moment acting on a freely suspended particle is zero [20], it may be shown that the spins of the ellipsoid around its axes are [23]: where e 0 ij denote the components of the strain-rate tensor of the external velocity field expressed in the body cos ψ sin 2θ sin 2ϕ + cos 2θ sin ψ sin ϕ e 0 13 =γ 2 1 2 sin ψ sin 2θ sin 2ϕ − cos 2θ cos ψ sin ϕ .
(2.12) Equation (2.10) expresses that the particle rotates by the action of the skew-symmetric part of the velocity gradient, but also the strain rates, which tend to reorient the particle towards the axes of strain. In addition, the tank-treading of the membrane modifies the particle spins. An equation for theΩ i is needed to close the system. The tank-treading rates are obtained by writing the kinetic energy balance for the particle, assuming that the membrane is a prestressed Kelvin-Voigt viscoelastic material [1,11,23]. The reason for using a Kelvin-Voigt material is mainly practical: it is the simplest model to account for viscoelastic effects in the membrane while keeping the analysis tractable. The calculation gives an equation that can be projected in the different directions to yield an expression of the tank-treading rates in each direction which is similar to the one found in one dimension by [1,11,37]: (2.13) These equations for the tank-treading rates show that tank-treading is favored if λ ef f is small (Λ is large) and depends on the strain rates e 0 ij to which the particle is subjected. A second term appears in the expressions The time between two images is 2.5γ −1 . The blue bead marks point P, the so-called marker, located on the symmetry axis at rest. Point P is fixed in the frisbee motion, with respect to the fixed frame. The symmetry axis (A) of the ellipsoid is marked by the red bead (it is fixed with respect to the ellipsoid, whatever the motion). The black bead is a Lagrangian point L of the membrane, whose movement follows the circulation of the membrane. It is initially along the x 1 axis. The frisbee motion is thus characterized by a fixed orientation with respect to the fixed frame, a fixed position of the marker P in the fixed frame, and a constant angular speed around the axis of symmetry.
of theΩ i and it involves the coordinates of the marker ξ i and the capillary number Ca * . It is an elastic term that notably prevents the membrane from tank-treading if membrane elasticity is strong enough. If Ca * is large (weak elastic effects), this elastic term is negligible. The system constituted by equations (2.8)-(2.13) is the MA model for the dynamics of an RBC in shear flow. MA [23] have solved this nonlinear ODE system numerically. The orientation of the ellipsoid and the circulation of the membrane are predicted along time, as a function k 1 , k 2 , Λ and Ca * and of the initial values for θ, ϕ, ψ, ξ i . After a possible transient response, orbit, frisbee, swinging or rolling dynamics are obtained, depending on the parameters and the initial conditions. Examples are provided in [23]. One example of phase diagram, ie the diagram of dynamics in the plane of the two varying parameters for a given particle, the (Ca * ,λ ef f ) plane, is displayed in Figure 1b. Figure 3 illustrates the dynamics of the ellipsoid in the frisbee regime, for the particular case a 1 = a 2 = 4.2375 µm, a 3 = 1.2511 µm, λ ef f = 1.0 (so that k 1 = 0.543, k 2 = 0.840 and Λ = 0.647) and Ca * = 1.20 [23]. The frisbee regime is particularly interesting, as it is obtained in the intermediate range of the parameters, so it is the only motion that can be obtained by transition from any other dynamics (see Fig. 1).

Limitations of the existing model
As already stated, the model can be solved numerically, but the parametrization of the orientation of the fluid ellipsoid with Euler angles leads to singularities in the advancement of the equations. For instance: where the (sin ϕ) at the denominator cannot be eliminated. As a consequence, the model is not well suited for mathematical analysis, due to the presence of singularities, in particular for ϕ = 0 or π, as already noted (Sect. 2.2). This has first motivated the development of a version without singularities, through the use of quaternions to parametrize the orientation of the ellipsoid. More significantly, the use of quaternion converts the model into a polynomial differential system. A whole set of new tools related to algebraic geometry is thus at our disposal and the model seems to be more tractable under this polynomial form. For example, in this article, we specifically make use of Gröbner basis to find the roots of polynomial equations to determine the fixed points of the system. In addition, we saw thatω 3 andΩ 3 are indistinguishable for ellipsoids of revolution axis e 3 . MA have arbitrarily chosen that the motion around the symmetry axis was a solid-body rotation, not a membrane circulation (Ω 3 = 0). However, this assumption prevents the existence of fixed points in the system: rolling, for instance, where ϕ = π/2 and the marker is fixed, is not a fixed point as the other Euler angles vary due to the non-zero value ofω 3 . As a consequence, to simplify the analysis, the opposite choice aboutω 3 andΩ 3 is made here: Assumption 2.1. Letω be the vectorial angular velocity of the rotation r t in the body frame, we assume thaṫ ω 3 = 0. The motion around the symmetry axis is then viewed as a pure membrane circulation.
Our starting point to introduce the quaternions is thus the system written in Section 2.5, albeit with: (2.14) instead ofω 3 = ζ 0 21 andΩ 3 = 0 in Section 2.5 [23]. This slight change has strictly no impact on the results, but allows to directly write frisbee solutions as fixed points of the system, as will be shown later.
The new version of the model is the object of the next section.

The model for the dynamics of an RBC without singularities
In this part, the model introduced in the former section is reformulated in the framework of quaternions for the parametrization of the rotations. Basics about quaternions are provided in Appendix B.

From Euler angles to quaternions, and vice versa
As explained in Appendix B, a unitary quaternion q = q 0 + q 1 i + q 2 j + q 3 k defines a rotation r q such that: r q (u) = quq. Conversely, a rotation r u,θ of angle θ along an axis u (where u is a unitary vector) is defined by the unitary quaternion cos(θ/2) + sin(θ/2)u.

Prerequisite: Lie bracket and hodge operator
Let us first introduce some specific notations, leading to more compact and manageable equations: Let A, B be two square matrices. The Lie bracket of (A, B) is the matrix: The Lie bracket is a bilinear, skew-symmetric operator. If A, B are symmetric matrices, [A, B] is skew-symmetric. There is a canonical isomorphism between the three-dimensional space A 3 (R) of (3, 3) skew-symmetric matrices and R 3 . This isomorphism will be called Hodge operator, and noted by a , in both directions: : The Hodge operator satisfies the following properties: Proposition 3.2 (Change of basis by rotation). Let R ∈ SO 3 (R) be a rotation matrix, and U ∈ R 3 . Then:

The differential system in the quaternions formulation
The ellipsoid being axisymmetric, the direction x 3 in the body frame will play a different role than x 1 and x 2 . The following diagonal matrix D is introduced, to handle this specificity: We may alternatively consider the difference between D and the identity matrix I: be a function such that r q is the ellipsoid's rotation and ξ the unitary coordinates of the marker in the body frame along time. Then (q, ξ) is solution of the differential system: where: -Ω ∈ R 3 is the coordinate vector, in the body frame, of the angular velocity of the tank-treading:

9)
-ω ∈ V is the angular velocity of the rotation r q . Its coordinate vector in the body frame is:

Frisbee solutions
The so-called frisbee motion is one of the possible dynamics for an RBC. In the model, it refers to a particular motion where the small axis of the ellipsoid e 3 reaches a fixed inclination, intermediate between the alignment with the vorticity axis e 3 and the shear plane ( e 1 , e 2 ). In terms of Euler angles, it means that ϕ is constant and different from 0, π/2 and π. In addition, in the frisbee motion, the marker P is fixed in the laboratory frame.

Definition
Since e 3 is fixed, and since, by assumption 2.1, the angular velocities around e 3 is zero, the whole solid body is motionless and the rotation r t is constant. The membrane can still circulate around this ellipsoid, rotating around the axis defined by the fixed marker P . However, considering only the marker and the solid body rotation one gets: Definition 4.1. A frisbee solution is a fixed point (q, ξ) : R → Q × S of the differential system (3.8).
In this section, the frisbee motion is investigated. We first derive the system of polynomial equations satisfied by a frisbee solution. Algebraic geometry tools are then used to compute frisbee solutions and compare them with results obtained by solving the dynamical system numerically.

Polynomial system for frisbee solutions
A frisbee solution (q, ξ) ∈ Q × S is defined by 7 parameters (q 0 , q 1 , q 2 , q 3 , ξ 1 , ξ 2 , ξ 3 ). We first reduce the number of parameters before presenting the equations they satisfy.

4.2.1.
The assumption q 3 = 0 Given one frisbee solution (q, ξ), it is possible to define a whole family of frisbee solutions in the following way: Let α be an arbitrary angle and consider the unitary quaternion b = cos(α/2) + sin(α/2)k, defining a rotation of angle α along the axis e 3 . The rotation matrix R b and the diagonal matrix D defined in equation (3.7) commute; it follows that (qb, t R b ξ) is also a solution of the differential system Consequently, we can choose any α to represent this whole family of solutions. A small computation shows that it is always possible to choose α such that q 3 = 0. We thus choose to limit the research of frisbee solutions to the cases where q 3 = 0.
In view of equations (3.1)-(3.3)), quaternions q = q 0 + q 1 i + q 2 j with q 3 = 0 correspond to the following Euler angles: which can be summed up in: Such a quaternion corresponds to a rotation of angle ϕ around the axis (cos(θ) e 1 + sin(θ) e 2 ) (see Appendix B).

Equations for a frisbee solution
As a fixed point, a frisbee solution (q, ξ) satisfies the following equation system: Since q is a unitary quaternion, it is invertible and 0 = 1 2ω q givesω = 0. In the body frame, we also getω = 0. The vectorsω andΩ are expressed in equations (3.10) and (3.9), respectively. Using equations (3.11)-(3.13) together with the assumption q 3 = 0 and q 2 0 + q 2 1 + q 2 2 = 1 yields for the angular velocity vectorω: and for the tank-treading rate vectorΩ: Remark 4.2. The termΩ 3 =γ 1 2 − q 2 0 is the tank-treading rate of the circulating membrane along the small axis e 3 . In terms of Euler anglesΩ 3 = −γ 2 cos ϕ (Eq. (4.1)). Proposition 4.3. Let q = q 0 + q 1 i + q 2 j be a fixed quaternion with coefficient of k equal to zero, and ξ = t ξ 1 ξ 2 ξ 3 be a fixed unitary vector, then (q, ξ) is a frisbee solution if and only if: Proof. The first two equations come from equation (4.3). The three following equations come from the computation of 1 γΩ × ξ = 0 using equation (4.4), the last two being the unitary conditions for q and ξ. Remark 4.4. Note that the frisbee condition does not depend directly onγ, but it depends on Ca * .

Symmetric solutions
As already mentioned, the symmetry of the flow yields symmetric solutions with respect to the shear plane. The aim of this section is to express the relations between two solutions, symmetric to the shear plane.
First, recall that due to the assumption q 3 = 0, the axis of the rotation r q is in the shear plane. In terms of Euler angles, that means: ψ = −θ. Consider then a frisbee solution with initial Euler angles (θ, ϕ, −θ) corresponding to a rotation matrix R, and a unitary marker ξ. As a consequence of the symmetry with respect to the shear plane, the frisbee solution with Euler angles (θ, π − ϕ, −θ) does exist (note that we cannot set −ϕ since we imposed 0 ≤ ϕ ≤ π). Let ξ be the unitary marker in the body frame and R the rotation matrix for this frisbee motion with angles (θ, ϕ − π, −θ). Then we have:  is the matrix of symmetry about the shear plane.
Proof. Relations on Euler angles are expressed in quaternions, yielding polynomial equations satisfied by the q i . We can then write down matrix R q , and compute ξ = SR q R q ξ.
When q 0 tends to 1, which corresponds to ϕ tending to 0, the ellipsoid tends to its position at rest, with the equatorial plane in the shear plane. The symmetric body with respect to the shear plane does the same. Both ellipsoid eventually merge with the ellipsoid at rest.
However, these formulas are not valid for q 0 = 1. In that case, ϕ = 0, and the rotation is the identity, which may be defined by an infinite set of Euler angles (θ, 0, −θ), θ ∈ R.

Critical cases
Before tackling the cases usually referred to as frisbee, ie where ϕ ∈]0, π/2[, we remark that the definition of frisbee solutions as a fixed point does not exclude critical cases where ϕ = 0, π/2. These cases are studied in this section and their characteristics are gathered in Table 2. Angle between e 3 and e 3 ϕ = π/2 ϕ = 0 Real part of the quaternion Height of the unitary marker P ξ 3 = 0 ξ 3 = 1 in the body frame Position of the unitary marker P ± e 3 ± e 3 in the laboratory frame

Description
Equator orthogonal to shear plane, Rolling in the shear plane, marker on the equator, on the e 3 axis marker at its position at rest

Conditions
No solution if k 1 Λ +
Proof. We are left to prove the two last assertions. As for the unitary marker, write ξ = ξ 1 i + ξ 2 j + ξ 3 k. We can check that qkq = ±ξ which means, as wishes, that r q ( e 3 ) = ±ξ (see Appendix B.2).
As for the tank treading rate, equation (3.9) givesΩ = Remark 4.9. In these equilibrium states, the marker point P should stay motionless on the equator. These states are not stable in the numerical simulations of the dynamical system.
These equilibrium states correspond to the well-known tank-treading as described by Keller and Skalak [20], for instance. The ellipsoid has its symmetry axis lying in the shear plane and its angle with respect to the flow direction is constant. Actually, the membrane circulates around the vorticity axis e 3 . As the marker is also aligned with e 3 , the membrane circulation takes place without change in the membrane stress, thus without angle variations [30]. This is thus not a swinging dynamics, but a pure tank-treading, that is a circulation with constant tank-treading rate and without angle oscillations, as would be obtained for an inelastic membrane. Note also that we retrieve the fact that this motion is only possible when the viscosity ratio is low enough (condition K ≥ 1), as already shown by Keller and Skalak [20] for ellipsoids with an inelastic membrane.

General solutions
We now investigate the cases ϕ ∈]0, π/2[. It is a priori not feasible to find explicit solutions for the polynomial system (4.5). Our approach is to find numeric solutions for fixed values of (Ca * , λ ef f ). Such a calculation can be realised with Elimination Theory and Gröbner basis [13] for polynomial systems (see Appendix C, Section C). All computations were made with the open-source formal calculus software SageMath (http://www.sagemath.org). The algorithm is provided in a SageMath procedure file provided as Supplementary Material.

Formal definition of the algebraic set of solutions
Let us first define a polynomial ring in seven variables A = R[q 0 , q 1 , q 2 , ξ 1 , ξ 2 , ξ 3 , u]. The additional variable u is used to exclude the already known critical cases q 0 = 0, q 0 = 1, ξ 3 = 0. For given values of Ca * and λ ef f , we define the polynomials: The first six polynomials come from System (4.5). Note that, since ξ 3 will not be zero, it is possible to remove the fifth equation in (4.5), which is a combination of the first two ones. As for the polynomial P 7 , if q 0 = 0, q 0 = ±1 or ξ 3 = 0, it is not possible to find u ∈ R such that P 7 = 0. In other cases, setting u = 1 q0(q 2 0 −1)ξ3−1 gives P 7 = 0. Let Sol(Ca * , λ ef f ) ⊂ R 7 be the set of common zeros of polynomials P 1 , . . . , P 7 . A point (q 0 , q 1 , q 2 , ξ 1 , ξ 2 , ξ 3 , β, u) in Sol(Ca * , λ ef f ) corresponds to a frisbee solution such that q 0 = 0, q 0 = 1, ξ 3 = 0.

Finite set of solutions
The system of interest consists in 7 polynomial equations in 7 variables. One may hope to obtain a finite set of solutions. Note that excluding the critical cases is essential: rolling solutions at ϕ = 0 may for instance be obtained for an infinite number of solutions with q 0 = 1 and q 2 1 + q 2 2 = 1. It turns out that fixing q 0 = 0, q 0 = ±1, ξ 3 = 0 is sufficient to get a finite set of solutions Sol(Ca * , λ ef f ) for any couple of parameters (Ca * , λ ef f ).

Possible values for q 0
We first start by determining possible values for q 0 for a frisbee solution at given (Ca * , λ ef f ) using the elimination algorithm described in Appendix C.3, for the ideal J = P 1 , . . . , P 7 of the ring A. We first eliminate every variables but q 0 , by providing ring A with the lexicographic order ≺ 0 such that: In the ordering above, we set q 0 as the smallest variable; the ordering of other variables does not matter. Using Buchberger's algorithm [43], implemented in the method J.groebner basis() of SageMath, we obtain a reduced Gröbner basis G 0 generating the ideal J.
As explained in Appendix C.3, we now set: Then, H 0 consists in polynomials of G 0 with no other variables than q 0 . Since G 0 is a reduced Gröbner basis and R[q 0 ] is a principal ideal, only three cases may occur: In our study, we are always in the last situation, and get a non-constant polynomial Q 0 ∈ R[q 0 ]. In other words, the calculation of the reduced Gröbner basis for the chosen order, combined with the elimination algorithm, yields a polynomial Q 0 in q 0 whose set of roots contains the values of q 0 for the possible frisbee solutions. We then define: C 0 = {real roots of Q 0 between 0 and 1}.
Thus C 0 is a finite set of 'possible values' for q 0 in a frisbee solution. Roots of Q 0 are computed with arbitrary precision using the method .roots(RR) in SageMath.

Possible values for q 1
In a second step, we provide ring A with another lexicographic order ≺ 1 such that: Here again, q 1 is carefully set as the smallest variable. Proceeding in the same way as for q 0 , we get a polynomial Q 1 ∈ R[q 1 ] and define: C 1 = {real roots of Q 1 between −1 and 1}.
Thus, C 1 is the finite set of 'possible values' for q 1 in frisbee solutions.

4.5.5.
Recovering solutions from (q 0 , q 1 ) One could perform the elimination algorithm described above for all other variables. However, it turns out that given a couple (q 0 , q 1 ) ∈ C 0 × C 1 , there exists at most one frisbee solution (q 0 , q 1 , q 2 , ξ 1 , ξ 2 , ξ 3 , β) such that ξ 3 ≥ 0. This can be proven by hand using the polynomials of System (4.7). Considering polynomials P 5 , k 1 P 4 − ξ 3 P 1 , and P 2 we necessarily have: If ξ 2 3 < 0, there is no solution with these values of (q 0 , q 1 ). If not, we set ξ 3 = 1 − 4Ca * k1 q 2 0 q 1 q 2 and obtain from further calculations: We are left to systematically check that these values of q i , ξ i , u = 1 q0(q 2 0 −1)ξ3−1 are roots of the seven polynomials of (4.7). Proceeding this way, we know which couples (q 0 , q 1 ) correspond to a solution. The whole process leads to a complete numerical computation of the set Sol(Ca * , λ ef f ) of frisbee solutions.

Algorithm
To summarize, the algorithm to calculate the frisbee solutions reads: -Define the polynomial System (4.7) and fix the parameters of the problem: the a i (hence k 1 and k 2 ), Ca * and λ ef f (hence Λ). -Order the variables such that q 0 ≺ 0 q 1 ≺ 0 ... and calculate the associated Gröbner basis G 0 using Buchberger's algorithm. Then use elimination theory for reduced Gröbner bases to obtain polynomial Q 0 (see Appendix C.3). -Calculate C 0 the set of the roots of Q 0 in the imposed range.
-Order the variables such that q 1 ≺ 1 q 0 ≺ 1 ... and calculate the associated Gröbner basis G 1 using Buchberger's algorithm. Then use elimination theory for reduced Gröbner bases to obtain polynomial Q 1 (see Appendix C.3). -Calculate C 1 the set of the roots of Q 1 in the imposed range.

Results
The geometry of the ellipsoid is fixed: a 1 = a 2 = 4.2375 µm and a 3 = 1.2511 µm (k 1 = 0.543, k 2 = 0.840) in order to match the cases already studied by MA. For this particular ellipsoid, we use the algorithm described in Section 4.5.6 to calculate the frisbee solutions. The computation of the sets Sol(Ca * , λ ef f ) is performed for every Ca * from 0.2 to 4, with a step of 2 × 10 −2 , and every λ ef f from 0 to 4, with a step of 2 × 10 −2 . It is found that we always get zero, one or two solutions with ϕ ∈]0, π/2[).
We first present the diagram showing the number of frisbee solutions, i.e. the cardinality of Sol(C * a , λ ef f ), as a function of Ca * and λ ef f , in Figure 4: -Light gray points corresponds to parameters with one frisbee solution with ϕ ∈]0, π/2[ (two solutions for the complete range of ϕ). -Dark gray points corresponds to parameters with two frisbee solutions with ϕ ∈]0, π/2[. -Uncolored regions correspond to parameters without frisbee solution.  Figure 1. For each value of (Ca * , λ ef f ), the number of frisbee solutions found by the algorithm based on the Gröbner bases is shown: in the light grey region, one frisbee solution is found in the range ϕ ∈]0; π/2[; in the dark grey region, two frisbee solutions are found. No solutions are found outside these regions. Figure 4 correspond to previous numerical predictions made from the MA model advanced numerically [23] and already reported in Figure 1b. Several conclusions may be drawn from the comparison between the solutions of the polynomial system and the long-term behavior of the system solved numerically. First, the agreement between the two methods in reproducing the frontiers between the frisbee dynamics and the others is a verification of the calculations from the original model to its new version using quaternions. The frisbee-rolling and frisbee-swinging frontiers are predicted by both methods. However, the frontier between the orbit region and the frisbee region in the simulations does not match the direct calculation of frisbee solutions: there is a region of the (Ca * , λ ef f ) plane where frisbee solutions exist, but are not found numerically (1.10 ≤ Ca * ≤ 1. 40 and λ ef f ≤ 1.0). In addition, the calculation of the frisbee solution shows a region where two different frisbee solutions exist for given (Ca * , λ ef f ). In the simulations, when a frisbee dynamics is obtained, it has been found to be unique. In order to clarify this qualitative comparison, the frisbee solutions themselves are examined. Figure 5 shows how the system goes from flipping over orbits (small Ca * ) to swinging (large Ca * ), by displaying for each solution the corresponding value of ϕ, the angle between the symmetry axis of the ellipsoid e 3 and the vorticity axis e 3 . For a given value of λ ef f , during the orbital drift, it is known that ϕ tends to 0 when Ca * increases. On the other hand, for large values of Ca * , swinging is found, for which ϕ = π/2. It is thus expected that in the region where frisbee solutions are found, the associated angle ϕ goes from ϕ = 0 to ϕ = π/2 when increasing Ca * . Three series of data are plotted for two values of λ ef f : 0.1 (Fig. 5a) and 1.0 (Fig. 5b). First, the value of ϕ found by solving the polynomial equations system (4.5) for the frisbee solutions are shown (solid lines). Then, two series of results from the ODE system are obtained, for two different initial orientation ϕ 0 : it is shown that the long-term dynamics may depend on the initial orientation by comparing the solutions obtained with ϕ 0 = 0.1π (+) and ϕ 0 = 0.4π (•).

Curves in
First, the behavior in terms of frisbee solution is different depending on the value of λ ef f . If λ ef f is small enough (here λ ef f = 0.1, Fig. 5a), there is a range of Ca * where two frisbee solutions exist. This is not the case for λ ef f = 1.0 (Fig. 5b), for which no more than one solution exists at a given value of Ca * . When solving the ODE system for low values of λ ef f , rolling is never found stable: only flipping over orbit, frisbee or swinging are observed. Finally, note that at fixed λ ef f and Ca * , the long-term dynamics may depend on the initial conditions. For instance, it is shown in Figure 5 that when solving the ODE system, frisbee and swinging solutions coexist and may be reached or not depending on the initial orientation.
We show that when a frisbee solution is reached, the agreement between the numerical solution and the theoretical frisbee solution in terms of the value of ϕ is perfect. It is also shown that only the branches where ϕ increases with Ca * seem stable. Figure 5a shows that for λ ef f = 0.1, no stable frisbee solutions with ϕ ∈ [0.325π; 0.5π] have been found. When the initial orientation of the ellipsoid is in this range (• series in Fig. 5a), it does not reach a frisbee dynamics but goes to the shear plane and swinging is obtained. This part of the solution branch seems unstable. Note also that the part of the solution branch for which Ca * < 1.369 has not been obtained. These frisbee solutions also seem unstable and flipping over orbits are obtained in this range of Ca * in numerical simulations. For the unstable parts of the solution branch, a series numerical simulations have been performed by imposing the frisbee solution found with the Gröbner bases as initial conditions. It has been found that the system does not remain at this position and reaches after some transient period a stable dynamics different from what was imposed. However, a systematic study of the stability of the solutions of Sol(Ca * , λ ef f ) would be needed to confirm these tests and better characterize the bifurcations of the system. Note that the frisbee solutions allow to clarify the singular behavior of the simulations at λ ef f = 0.1, where the solution suddenly transitions from orbits to frisbee and from frisbee to swinging, exploring a limited range of ϕ in the frisbee motion. Our analysis shows that frisbee solutions exist for all values of ϕ, but may be unstable.

Conclusion
This work has been performed in the context of the modeling of the movements of a red blood cell in shear flow. An existing family of models [1,11,20,23,37], representing the RBC as a fixed-shape fluid ellipsoid with a membrane allowed to slide around the shape, had been shown to provide results in good agreement with experiments and simulations and enable explanations of many features of the RBC motion. Their most complete version, recently proposed by Mendez and Abkarian (MA) [23] is a 3-D model able to predict the off-shear plane dynamics of RBCs, which was not the case for the former models [1,11,20,37]. Note also that the 3-D MA model is an extension of the former 2-D models [1,11,20,37] and retrieves all the dynamics predicted by those 2-D models, if the symmetry axis is forced to lie in the shear plane. However, this 2-D dynamics is not always stable: when the intermittent regime (alternation of swinging and tumbling movements) is found in 2-D, the 3-D model predicts the existence and stability of the frisbee motion. Note that the intermittent regime has also been found to vanish for quasi-spherical particles if the fixed shape assumption is relaxed [42]. However, whatever its predictive capabilities, the complexity of the MA model, which is a system of 5 coupled nonlinear differential equations, makes its physical interpretation quite difficult. One strategy when analyzing dynamical systems is to study the fixed points of the systems, its limit cycles and their stability. However, the MA model is not well suited for such an analysis, notably because the equations are singular.
The first result of this paper is to propose a new version of the MA model, where the orientation of the ellipsoid, originally parametrized using Euler angles [23], is parametrized using quaternions. A new framework for predicting the dynamics of a RBC is thus proposed, based on quaternions, for the first time. In addition, some assumptions are slightly changed to feature steady-state solutions, which did not exist naturally in the previous system. Such a model being well suited for mathematical analysis, this work can be viewed as a first step towards a complete characterization of the system. One advantage of using quaternions is that the resulting equations are polynomials. This opens the way to the use of powerful techniques from Algebraic Geometry.
An example of this use is given in the paper. With the idea of further characterizing the MA model [23], the fixed points of the system are sought for. Such dynamics actually correspond to a specific motion of RBCs in shear flow, referred to as 'frisbee'. Numerical advancement of the dynamical system [23] had shown that frisbee is obtained for certain values of the parameters of the problem, the viscosity ratio and the capillary number. However, from numerical simulations, it is impossible to know if the frisbee dynamics found were the only ones or if other frisbee solutions exist. We first show that frisbee solutions are the roots of a system of 6 polynomials. After isolating an infinite set of trivial solutions, frisbee solutions are sought for by calculating the Gröbner basis of the ideal generated by the 6 polynomials of the system. Then, by elimination, the frisbee solutions are determined. To the best of our knowledge, this is the first time that non-trivial off-shear plane solutions relevant to RBC dynamics are determined analytically. Such solutions notably provide reference results for the validations of numerical simulations. Note that the method presented here can be used for any situation which can be formulated as a fixed point, for instance for other external flows.
Results show that the predictions of the numerical advancement of the MA model and the result of the calculation of the frisbee solutions by elimination are consistent, but interesting differences are obtained. All the solutions found numerically are indeed frisbee solutions of the system. However, there exist other solutions than those found numerically, showing the existence of unstable frisbee solutions. When studying the transition from the orbit dynamics to the swinging dynamics at fixed values of the viscosity ratio, numerical simulations show different series of behaviors when increasing the capillary number, depending on the viscosity ratio. When the viscosity ratio is high, frisbee is not obtained. But for intermediate values of the viscosity ratio, for instance, an RBC first experiences the flipping dynamics and the associated orbital drift towards rolling, then frisbee, then swinging. Transitions are smooth between the different dynamics, the frisbee solutions continuously reorienting from the rolling position (the result of the orbital drift) to the swinging position, with increasing capillary number. On the contrary, for small values of the viscosity ratio, the transitions from flipping to frisbee and from frisbee to swinging is abrupt. The calculation of the frisbee solutions explains this discontinuous transition by the existence of a branch of solutions passing continuously from the rolling position to the swinging position, but with parts of this branch being unstable.
We view the present work as a first and rather technical step towards a thorough stability analysis of the model, which is now all the more interesting that the present study has evidenced the existence of unstable solutions. With respect to that objective, the change of formulation to write a system based on quaternions is an important and interesting milestone. We have notably shown how the resulting dynamical model is well suited for the application of Algebraic Geometry techniques to determine some solutions of the problem. Our intention is to undertake a stability analysis of these solutions and fully characterize the bifurcations of the system to complete the investigation of this fluid-structure problem modeling the dynamics of a red blood cell under flow. Finally, the method used in this work is not limited to the case of the shear flow and the dynamics of RBCs or axisymmetric capsules in any general linear flow could be treated in the same way.
Appendix A. Parameters of the model for the dynamics of a fluid ellipsoid with circulating membrane Many geometrical factors are involved in the model. For the sake of completeness, they are defined in appendix, to lighten the main text. The ellipsoid of interest is an axisymmetric oblate ellipsoid of semi-axes a 1 = a 2 > a 3 . The quantity a 0 = (a 1 a 2 a 3 ) 1 3 is introduced, allowing the definition of α i as α i = a i /a 0 . We define f 1 , f 2 and f 3 as: and An explicit formula for g 1 can be given: setting a = a1 a3 , we get α 1 = α 2 = a 1/3 , α 3 = a −2/3 , and g 1 (a) = ∞ 0 ds (a 2/3 +s) 2 (a −4/3 +s) 3/2 . Substituting s by t 2 − a −4/3 gives g 1 (a) = ∞ a −2/3 2dt (t 2 +a 2/3 −a −4/3 ) 2 t 2 , which can be computed by partial fraction decomposition. A formal calculus software, such as SageMath, can be used to get the primitive function and calculate the integral of interest. We eventually obtain: In addition, the in-plane capillary number Ca * in the model is not the classical capillary number Ca = µ extγ a 0 /G s . Ca * results from an explicit calculus of the stress in the membrane [1,2,23] and has a nontrivial expression. In the model [23], the volumes of the ellipsoid and of the membrane are denoted by V and V m , respectively. G denote the shear modulus of the membrane, described as a Kelvin-Voigt viscoelastic 3-D material, of small yet non-zero thickness. C is a constant that varies between 0 and 1 and allows to account for the membrane prestress. Details about the C constant are provided by Dupire et al. [11]. The in-plane capillary number Ca * is expressed as (the minus sign is to obtain a positive number, as f 1 and f 3 are of opposite sign): It characterizes a competition between the external viscous stress µ extγ and the elastic restoring forces that prevent a prestressed membrane from circulating. Finally, the other main parameter of the problem, λ ef f , is involved in the parameter Λ, which reads A simple computation gives the matrix R q of r q in the basis (i, j, k) = ( e 1 , e 2 , e 3 ): Note that, since e i = r q ( e i ), the matrix of r q is also R q in the basis of the body frame.

B.3 Vectorial angular velocity
Let q : t ∈ R → q(t) ∈ Q be a C 1 family of unitary quaternions parametrized by time t. At any time t, the angular velocity vectorω of the rotation r q(t) is uniquely defined by the property: We can show that:ω = 2 dq dt q = 2qq. (B.3) Note that the product has to be taken in H, butω is a pure quaternion, an element of the Euclidian space V .
Appendix C. Elimination theory and Gröbner basis for polynomial systems C.1 The aim of elimination Consider a system of r polynomial equations: P 1 (x 1 , . . . , x n ) = 0, · · · , P r (x 1 , . . . , x n ) = 0 (C.1) where the P i are elements of the polynomial ring R[x 1 , . . . , x n ]. If each P i is of degree one, this is a system of affine equations. It is then possible to express one variable (say x n ) in terms of the others. Then, substituting x n by this expression gives a new system in the variables (x 1 , . . . , x n−1 ) only. The unknown x n has been eliminated, and can be recovered from the solution of a smaller system with less variables; no information is lost during this process.
The aim of elimination theory for general polynomials is, given the system (C.1), to write down a new system: . . . , x n−1 ) = 0, · · · , Q s (x 1 , . . . , x n−1 ) = 0 (C. 2) in which the variable x n no longer appears, such that if (x 1 , . . . , x n ) is solution of (C.1) then (x 1 , . . . , x n−1 ) is solution of (C.2), and such that knowing the solution of (C.2) gives as much information as possible on the solution of the initial system.
Example C.1. As a small example, consider the intersection of a circle and two lines in R 2 : P 1 (x, y) = x 2 + y 2 − 4 = 0, P 2 (x, y) = xy − y 2 = 0 (C.3) Elimination by hand shows that, if (x, y) is solution of (C.3), then y is a root of the polynomial Q 1 (y) = y 3 − 2y (elimination of x), whereas x is a root of Q 1 (x) = x 4 − 6x 2 + 8 (elimination of y). Factorizing these polynomials, we get at most three possible values for y (0, and ± √ 2/2) and four possible values for x (± √ 2, and ± 4). We then can test all the possible couples (x, y) to find the correct solutions.
One may also consider the elimination of several variables, which can be done successively, or simultaneously depending on the algorithm. Eliminating all variables but one, say x 1 , gives a single polynomial in x 1 , and thus a finite set of possible values for x 1 , as in Example C.1. Repeating the process for each variable yields a finite set of possible values for x 1 , . . . , x n . We then can check numerically which of these potential solutions (x 1 , . . . , x n ) really satisfy the original system.

C.2 Monomial order and Gröbner basis
Several Elimination algorithms, using resultants or the Gröbner bases for example, are at our disposal. In this section, we shortly present the notion of Gröbner bases [13], that will be used for our elimination algorithm.
A monomial order for the polynomial ring R[x 1 , . . . , x n ] is a total ordering ≺ on the monomials x k1 1 · · · x kn n , compatible to multiplication and such that 1 = x 0 1 · · · x 0 n is the smallest element. If n = 1, there is only one monomial order, given by the degree. If n > 1 one may choose a lexicographic order, and the monomial order is obtained by simply ordering the n variables.
Suppose now that R[x 1 , . . . , x n ] is endowed with a monomial order ≺. For any polynomial P , we can consider its leading term lt(P ). If n = 1, this is the term with the highest degree.
Let I =< P 1 , . . . , P r > be an ideal generated by r polynomials, we define: lt(I) = {ideal generated by leading terms of the polynomials in I}.
The ideal lt(I) is a monomial ideal, that is, an ideal generated by monomials. Such ideals are very tractable in formal computing, and lots of algorithms are constructed from them. However, it is not true in general that lt(< P 1 , . . . , P r >) =< lt(P 1 ), . . . , lt(P r ) > .
We only have: Definition C.2. Let I be an ideal of A. There exists a finite set of polynomials G = {Q 1 , . . . , Q s } ⊂ I such that G generates I and lt(I) =< lt(Q 1 ), . . . , lt(Q s ) > . Such a set of polynomials is called a Gröbner basis of I.
One may also define a reduced Gröbner basis. This is a Gröbner basis G such that, if Q i , Q j ∈ G, then lt(Q i ) does not divide any term of Q j . Reduced Gröbner basis always exist. Given a finite set of generator (P 1 , . . . , P r ) of I, the Buchberger Algorithm [43] provides a reduced Gröbner basis for I.

C.3 Elimination with Gröbner basis
When used with a lexicographic order, the Gröbner basis provides a simple elimination algorithm for polynomial systems: Consider System (C.1), and define the ideal I =< P 1 , . . . , P r > of R[x 1 , . . . , x n ] generated by I. Assume that R[x 1 , . . . , x n ] is endowed with the lexicographic order ≺ such that: x 1 ≺ · · · ≺ x n , and let G = (Q 1 , . . . , Q s ) be a Gröbner basis of I.
By definition of the lexicographic order, any polynomial Q of R[x 1 , . . . , x n ] such that lt(Q) ≺ x k is a polynomial in the variables x 1 , . . . , x k−1 only: the variables x k , . . . , x n do not appear. This simple observation, and the use of Gröbner basis leads to the Definition: Definition C.3. For 1 ≤ k ≤ n − 1, consider the set G k = {Q ∈ G, lt(Q) ≺ x k+1 }.
Then the ideal I k =< G k > generated by G k is called the k-th elimination ideal of I; it is an ideal of R[x 1 , . . . , x k ].
Thus, in order to eliminate a single variable, one may consider the ideal I n−1 and the system: {Q(x 1 , . . . , x n−1 ) = 0, Q ∈ G n−1 . This is the system (C.2). In our cases, we prefer to eliminate all variables but the smaller one, x 1 . Then, we consider the ideal I 1 and the set of polynomials G 1 of R[x 1 ]. Since G is a reduced Gröbner basis, and R[x 1 ] is a principal ring, it follows that G 1 is either empty (in which case the System (C.1) has no solution) or contains only one polynomial. System (C.2) thus reads: The case Q = 0 is the worst possible, since no information is obtained on x 1 . In order to avoid this, it is necessary to remove infinite set of solutions. This is what is done in Sect. (4.5).
If Q is not zero, it is possible to numerically find roots of Q and a finite set of possible values for x 1 . Repeating the process with different lexicographic orders gives a finite set of possible values for each variable, and define a finite set of possible solutions (x 1 , . . . , x n ) that have to be tested numerically.