{'fraction_non_alphanumeric': 0.03337357373729507, 'fraction_numerical': 0.012558900268702052, 'mean_word_length': 4.53616470841867, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electrohydrodynamics of Three-Dimensional Vesicles: A Numerical Approach

A three-dimensional numerical model of vesicle electrohydrodynamics in the presence of DC electric fields is presented. The vesicle membrane is modeled as a thin capacitive interface through the use of a semi-implicit level set Jet scheme. The enclosed volume and surface area are conserved both locally and globally by a new Navier-Stokes projection method. The electric field calculations explicitly take into account the capacitive interface by an implicit Immersed Interface Method formulation, which calculates the electric potential field and the trans-membrane potential simultaneously. The results match well with previously published experimental, analytic and two-dimensional computational works.

1. Introduction. Giant liposome vesicles are enclosed bag-like membranes composed of lipid bilayers. These vesicles share many similarities in composition and size with other biological cells such as red blood cells. Experimental observations have demonstrated that lipid vesicles exhibit a striking resemblance to more complicated biological cells in terms of their equilibrium shapes and dynamic behavior in various fluid conditions [1,2]. These make liposome vesicles a robust model system to study the behavior of more complicated non-nucleated biological cells.
It is also quite easy to artificially create these soft particles in a laboratory setting; in an aqueous solution lipid molecules will self-assemble into two-dimensional sheets, which then fold and curve in three-dimensional space and form a fluid-filled vesicle [3]. The ease of liposome vesicle formation, in addition to their bio-compatibility, has led to vesicles being proposed as building blocks for various biotechnologies such as directed drug and gene delivery [4,5] or biological microreaction [6,7].
Vesicles interacting with external flows have been of major recent interest. Three well-known types of vesicle dynamics in shear flow are tank-treading, tumbling and a transient state called trembling or vacillating breathing. These behaviors have been extensively studied in theory [8,9,10,11,12] and experiments [13,14,15]. Several numerical simulations of vesicle dynamics have also appeared in the literature. Among the numerical studies in two dimensions are models using the boundary integral method [16,17], the phase field [18,19,20], a coupled level set and projection method [21,22], a coupled level set and finite-element method [23,24] and the lattice Boltzmann method [25]. A few three-dimensional studies have been also reported using the phase-field approach [26], boundary integral method [27,28,29] and fronttracking [30]. Studies have shown the dynamics of the vesicle depends on three major parameters: the viscosity ratio between the enclosed and surrounding fluids, the shear rate and the reduced-volume, which is defined as the ratio of the volume of the vesicle to the volume of a sphere with the same surface area as the vesicle.
The use of electric fields, in combination with fluid flow, has been proposed as a powerful method to direct the behavior of vesicles towards a wide range of biotechnological applications. Weak electric fields have found applications in cell manipulation techniques such as electrofusion [31], tissue ablation [32], wound healing [33], and in the treatment of tumors [34]. Strong electric fields induce electro-poration in vesicles through the formation of transient pores in the membrane and could play a role in novel biotechnological advances such as the delivery of drugs and DNA into living cells [35,36].
Recent experiments have reported on the topological behavior of vesicles subjected to either a DC or AC electric field with different intensity, frequency and duration of exposure to the field [37,38,39,40]. Experiments show that depending on the conductivity and permittivity differences between the enclosed and surrounding fluids the vesicle undergoes various shape transformations, such as prolate (major-axis aligned with the electric field) to oblate (major axes are perpendicular to the field). One very interesting phenomena in this context is the dynamics of an initial prolate vesicle in a strong DC field with a small enclosed fluid conductivity as compared to the conductivity of the surrounding fluid. It is theoretically expected that as time progresses the vesicle transitions to an oblate shape first, then evolves back into the prolate shape. During such transition nearly cylindrical shapes with high-curvature edges were observed for a vesicle subjected to strong pulses [38]. This behavior was attributed to the presence of salt in the solution and not due to the membrane characteristics. However in [41] the poration of the vesicle membrane was proposed as the possible explanation for both vesicle collapse and cylindrical deformations.
In addition to recent experimental works, theoretical models have also investigated the electrohydrodynamics of nearly-spherical vesicles. Leading-order perturbation analysis has been employed to obtain reduced models in the form of ordinary differential equations [42,43,44,45]. In another work a spheroidal shell model has been used to investigate the morphological change of the vesicle in AC electric field [46].
Despite numerous theoretical investigations, numerical studies of the vesicle electrohydrodynamics are rare. In a recent work, a boundary integral method was employed to study different equilibrium states of a two-dimensional vesicle in the presence of a uniform DC electric field [47]. A Prolate-Oblate-Prolate (POP) transition was captured for a vesicle with an inner-fluid conductivity smaller than the surrounding region.
Currently, there is still a gap between the morphological changes observed in experiments, what theoretical models predict, and methods to control this behavior. This gap needs to be addressed before vesicles can form a building block for electrohydrodynamic based microfluidic systems and technologies. Part of the difficulty arises from modeling the complex physics of the vesicle electrohydrodynamics and challenging phenomena such as membrane poration or fusion. Compared to hydrodynamics investigations, the electrohydrodynamics of the vesicle has substantially faster dynamics with much larger deformations which makes the numerical modeling nontrivial. Hence a lack of thorough numerical investigation with different material properties and electric field parameters still remains.
In a previous work by the authors the Immersed Interface Method (IIM) was utilized to solve for both the electric potential and trans-membrane potential around a vesicle in three dimensions [48]. The jump conditions for potential and its first and second derivatives on the interface were determined and utilized to obtain accurate electric potential solutions for a three-dimensional vesicle with an arbitrary shape. In this paper the electric field model presented in this recent work is combined with a projection-based hydrodynamics solver and a semi-implicit jet scheme for capturing the interface. This model is used to study the electrohydrodynamics of three dimensional vesicles in general flow. The dynamics of the vesicle is determined by the interplay between the hydrodynamics, bending, tension and electric field stresses on the membrane. This is one of the first attempts at investigating the electrohydrodynamics of vesicles in three-dimensions. As one may expect, for this kind of physics a three dimensional model will result in a much richer and wider range of topological changes. The differences compared to a two dimensional model will be subtle and important at the same time. In addition to this, a full Navier-Stokes system of equations is considered here unlike all the previous models which investigate the problem in the Stokes region. For vesicles in strong DC field, the approximate velocity may sometimes exceed 0.01 m/s [41], resulting in non-trivial Reynolds numbers, and therefore, the Stokes assumption for the fluid flow might not always hold. Moreover, the method presented here allows for studying deflected vesicles with smaller reduced volumes than what theoretical models normally address and is able to predict the type of deformations observed in experiments.
The remainder of this paper is organized as follows. In Section 2 the physical system and formulation behind the electrohydrodynamics of the vesicle are described. In Section 3 the numerical algorithm is presented. This will be followed by sample numerical results in Section 4. Final remarks and possible future work then follow.
2. Theory and formulation. Let a three-dimensional vesicle of encapsulated volume V have a surface area of A. Deviation from a perfect sphere is measured by a reduced volume parameter, v, defined as the ratio of the vesicle volume to the volume of a sphere with the same surface area: v = 3V /(4πa 3 ) where a = A/4π is the characteristic length scale. The typical size of the vesicle is a ≈ 10 − 20 µm while the thickness of the bilayer membrane is d ≈ 5 nm [49]. Due to the three-orders of magnitude difference between the vesicle size and the membrane thickness the membrane will be treated as an infinitesimally thin interface separating the inner and outer fluids. The vesicle membrane is assumed to be impermeable to fluid molecules and the number of lipids on the membrane does not change over time, while the surface density of lipids at room temperature is constant [49]. These two conditions result in an inextensible membrane with constant enclosed volume and global surface area, along with local surface incompressibility. Another important feature of the membrane is its electrical insulating property and impermeability against ionic transfer [40,42,47]. Therefore, in the presence of an external electric field the membrane acts as a capacitor. This capacitive property is an important factor in studying the dynamics of vesicles in the presence of electric field.

Electric Field Equations.
A schematic of a vesicle exposed to an external electric field is illustrated in Fig. 2.1. Different properties inside (−) and outside (+) of the membrane are shown in the figure. The embedded fluid (Ω − ) is separated from the surrounding fluid (Ω + ) by the vesicle membrane represented as Γ. This membrane is assumed to be made of a charge-free lipid bilayer with uniform and constant capacitance C m and conductivity G m . The vesicle is suspended in a fluid of conductivity s + and permittivity + . The enclosed region is assumed to have a different conductivity s − and permittivity − . It is assumed that properties are constant in each fluid, with a finite jump occurring across the membrane.
The leaky-dielectric model is assumed to hold in the whole domain [50]. This results in no local free charge in the bulk fluids, and thus the electric field (E) is irrotational and given as the negative gradient of the electric potential (Φ) However, the presence of the capacitive interface leads to a discontinuity in the potential across the membrane where V m is the transmembrane potential and can be obtained from the conservation of current density across the membrane [45,51], where n is the outward facing normal vector on the interface (pointing into the Ω + fluid), ∇ s is the surface gradient operator, u is the underlying fluid velocity and Q ± is the induced charge on the top and bottom side of the membrane. In the limit of fast bulk charge relaxation and negligible charge convection the time evolution of E ± and the surface convection of Q ± can be ignored [45]. Furthermore, continuity of the Ohmic current J = sE in the normal direction across the membrane can be written as Due to the difference in electrical conductivities of the inner and outer fluids charges accumulate on the interface at different rates. The effective charge density, Q, induced on the membrane as a result of an external electric field is then expressed as the jump in the normal component of displacement field D = E: (2.6) The electric force τ el acting on the membrane is obtained from the Maxwell Tensor T el : where [ ] in the above relation and through the rest of this paper shows the jump across the interface. For an arbitrary quantity f at location x Γ on the interface a jump is mathematically defined as Note that when discussing jumps in quantities at the interface a superscript + or − does not denote a bulk value, but the value as one approaches the interface.

Fluid Flow Equations.
Assume that both bulk fluids are Newtonian and incompressible with a matched density ρ. Consequently they both satisfy the Navier-Stokes equations where u is the velocity vector, D/Dt is the material (total) derivative, and T hd is the bulk hydrodynamic stress tensor defined as The fluid flow in each region is coupled via the conditions on the inextensible membrane. The velocity is assumed to be continuous on the surface, [u] = 0. However, due to the forces exerted by the membrane the hydrodynamic stress undergoes a jump across the interface of the two fluids. This condition is obtained by balancing the hydrodynamic and electric stresses with the bending and in-extension tractions of the membrane, τ hd + τ el = τ m + τ γ on Γ, (2.11) where τ hd = n · [T hd ] is the normal component of the hydrodynamic stress while τ m and τ γ are the bending and tension traction forces, respectively. These two forces together constitute the total membrane force per unit area and are calculated by taking the variational derivative of the total energy of the membrane [8]. For a three dimensional vesicle and neglecting spontaneous curvature the ultimate forms of the bending and tension stresses are found to be [52] In this relation, H = ∇ · n is twice the mean curvature (the total curvature), K = ∇ · [n∇ · n + n × (∇ × n)] is the Gaussian curvature, κ c is the bending rigidity, κ g is the Gaussian bending rigidity, γ is the tension and ∇ s = P ∇ where P = (I − n ⊗ n) is the surface projection operator. The last condition to consider is the local area incompressibility constraint on the membrane. This condition is enforced by ensuring that the velocity is surfacedivergence free on the membrane: (2.13) Note that the divergence conditions, ∇ · u = 0 and ∇ s · u = 0, should be sufficient to conserve both volume and area. However, in practice it is difficult to ensure that global surface area and enclosed volume are sufficiently conserved [20]. Therefore, when developing the numerical method in Sec. 3 global conservation will be explicitly enforced as well.
2.3. Description of the Interface. The motion of the interface is tracked using a gradient-augmented level set method. First developed by Osher and Sethian [53] the level set method is based upon an implicit representation of an interface as the zero level set of a higher dimensional function. This was later extended to explicitly include information about the gradients of the level set [54]. The use of gradient information allows for the accurate determination of the interface location and curvature information away from grid node locations [54].
Using the same notation for the interface as Fig. 2.1 the level set function φ is defined as (2.14) while the level set gradient field is defined as The level set value is chosen to be negative in the Ω − and positive in the Ω + domain. This representation has the advantages of treating any topological changes naturally without complex remehsing, and the ability to calculate geometric quantities from the level set function and its derivatives. For example, the unit normal to the surface, n, can be easily computed as The motion of the interface under the flow field u is modeled by a standard advection equation The evolution of the gradient field is obtained by taking the gradient of the level set advection equation, Using the level set function one is able to define the viscosity at any location in the computational domain, x, with a single relation, (2.19) where H is the Heaviside function. Different forms of the Heaviside approximations have been proposed [55,56,57]. Here an accurate method based on an integral calculation [58] is used. The Heaviside is defined as Similarly the Dirac delta function is expressed as These definitions will be used to localize the contributions of interface forces in Navier-Stokes equations.

Continuum Surface Force
Model. The definitions of the Dirac and Heaviside functions in Sec. 2.3 allow for the use of a single equation to describe the dynamics of the fluid over the entire domain [59]. This is accomplished by writing the singular contributions of the bending, tension, and electric field forces as localized body force terms, similar to what has been done for two dimensional vesicles [21,22]. This results in the following single-fluid formulation, Note that in this formulation all surface quantities, such as tension or the jump in the Maxwell stress tensor, are calculated on the interface and extended such that they are constant in the direction normal to the interface.

Dimensionless Parameters.
When an electric field is applied to the system charges from the bulk fluids migrate towards the interface. The time scale associated with this migration process is given by the bulk charge relaxation time [50,60], For an ion-impermeable membrane the characteristic time scale associated with the charging process is given by [36,61] where a is the characteristic length scale. The applied electric stresses act to deform the vesicle membrane on the electrohydrodynamic time scale given by where the applied electric field has a strength of E 0 and η = µ − /µ + is the viscosity ratio. The bending forces act to restore the membrane to an equilibrium configuration on the bending time scale, Finally, the response time of a vesicle to any externally applied shear flow is simply the inverse of the shear rate,γ 0 , It is beneficial to estimate the order of magnitude for the times at which different physical processes occur. Typical values of the physical properties have been reported [38,40,62,63].
Using these values the bulk charge relaxation time is t c ≈ 10 −6 s, the bending time scale is t κc ≈ 160 s, the membrane charging time scale is t m ≈ 2×10 −3 s, the time scale for shear is tγ = 1 s and for an electric field of E 0 = 10 5 V/m the electrohydrodynamics time scale will be t ehd ≈ 10 −4 s.
It is worth noting that the bulk charge accumulation on the interface happens in a much faster time than any other events. Therefore, it can be concluded that the electric field adjusts to a new configuration of the vesicle and fluid almost instantaneously and the quasi-static assumption for the electric field in Eq (2.1) is valid. It is also important to note that with this parameter set the electrohydrodynamics time scale, t ehd , is faster than the membrane charging time scale, t m . If this is not the case then the vesicle membrane will not be able to respond quickly enough to the applied forces, and only small deformations will be observed [40,43].
2.6. Nondimensional Model. Given a characteristic length a and time t 0 the characteristic velocity is given by u 0 = a/t 0 . Material quantities such as viscosity and permittivity are normalized by their counterparts in the (outer) bulk fluid (i.e. µ = µ/µ + andˆ = / + ). The dimensionless fluid equations are written as where dimensionless quantities are denoted by a hat. The velocity at the boundary of the domain is given as u ∞ = χŷ where χ =γt 0 is the normalized applied shear rate. The uniform DC electric field at the boundary is imposed as Φ ∞ =Êŷ whereÊ is the normalized strength of the applied electric field. The dimensionless parameters are defined as follows. The Reynolds number is taken to be Re = ρu 0 a/µ + , while the strength of the bending is given by a capillary-like parameter, Ca = t κ /t 0 and the strength of the electric field is given by the Mason number, M n = t 0 /t ehd .
The time-evolution equation of the transmembrane potential is nondimensionalized in a similar manner: where the dimensionless membrane capacitance is given asĈ m = (C m a) / (t 0 s + ), the dimensionless membrane conductivity isĜ m = (G m a) /s + and the conductivity ratio is expressed as λ = s − /s + . For simulations in the absence of shear flow the most appropriate time scale is the membrane-charging time, t m [40]. Hence for all the electrohydrodynamic computations with no imposed shear flow (χ = 0) in Sec. 4.2 the simulation time scale is set to t 0 = t m and the dimensionless parameters Re, Ca and M n are calculated accordingly. On the other hand, using the membrane charging time scale for the hydrodynamic simulations (M n = 0) will lead to abnormally huge number of iterations. Therefore for the hydrodynamic computations in Sec. 4.1 the simulation time scale is set to the time scale associated with shear flow (t 0 = tγ). This time scale has been also used for the investigation of the vesicle dynamics under combined shear flow and weak electric field in Sec. 4.3. This seems to be a rational choice as the time scale of the applied weak electric field is in the same range as the one for the shear flow.
Using the typical experimental values of the physical parameters discussed earlier in Sec. 2.5 the non-dimensional numbers for simulations in the absence of shear are Clearly, at this time scale the Reynolds number is not small enough to justify the Stokes approximation while the bending contribution is almost negligible. Despite the small contribution from the bending in this situation it is kept in the formulation as it may become important during other time-scales, such as when shear is applied. For simplicity the hat notation in the equations henceforth dropped. All the above mentioned parameters along with the viscosity ratio, η, and reduced volume, v, will be used to investigate the dynamics of a three-dimensional vesicle in different flow conditions and strength of electric field.
3. Numerical Methodology. In this section the specific numerical implementation of the non-dimensional evolution equations described in Sec. 2.6 will be given.
3.1. Electric field solution: Immersed Interface Method. Using a second order time discretization scheme for the time-evolution equation of the transmembrane potential in (2.31) and treating the term V m implicitly one obtainŝ The right hand side of Eq. (3.1) shows that the solution of the electric field at the current time step, t n+1 , is needed in order to compute V n+1 m . To solve the electric field equation with mismatched conductivity, a three-dimensional Immersed Interface Method (IIM) with time-varying solution jumps is employed. To achieve second-order accuracy in space, the jump conditions for the electric potential are derived up to the second normal derivative. First define r = −n · E α with α representing the domain with lower conductivity (i.e. if s − < s + then α = −). Using this definition and employing the relations given in Eqs. (2.3), (2.5) and (3.1) the complete set of jump conditions can be derived as (for details of the derivation see [48]) where s −α represents the higher of the two fluid conductivities. The field equation Eq. (2.2), along with the above set of jump conditions are used to solve for r and Φ simultaneously. The new trans-membrane potential is then updated using Eq. (3.1), noting that r provides the solution to the electric field in one of the domains. The technique to solve the resulting linear system plus further details of the implementation and sample results can be found in Ref. [48].
To calculate the electric field contribution to the Navier-Stokes equation, Eq. (2.30), the normal component of the Maxwell stress tensor is needed. This tensor is computed at the interface for each domain, using the corrections in jump conditions to obtain second-order accuracy. The difference between the normal component is then calculated, giving the electric field force on the interface.

Level Set Advection.
Many of the equations given in previous sections require high order geometric quantities, such as the curvature. To ensure that surface quantities are smooth a semi-implicit level set advection scheme is utilized. The method used here is an improvement of the semi-implicit, gradient augmented level set [64], and is briefly described here. Consider the advection of the level set field, φ, and the associated gradient field, ψ, due to an underlying flow field, u. Seibold et. al. generalized this type of advection by considering a sub-grid of points centered on a grid point, which provides a jet of level set information [65]. It was shown that using a sub-grid spacing of O(δ 1/4 ), where δ is the floating point operation accuracy, provided the optimal accuracy for linear advection equations. As the equations which model the electrohydrodynamics of vesicles require that a smooth level set field be maintained at all times, this high accuracy solution is relaxed. In this work the jet of level-set information surrounding a grid point is the cell centers, Fig. 3.1(a). The advection of the level set jet proceeds with the following three steps: 1) advection of cell-centers, 2) smoothing at cell-centers, and 3) projection onto grid points. Unlike the original semi-implicit, gradient-augmented level set method [64] this has the advantage of only requiring a single smoothing operation, as described below.
The advection of cell-centers proceeds by using a Lagrangian advection scheme: where x c is the cell-center location,φ c is the tentative level set value at the cell center, and I 3 (φ n , x d ) is a cubic interpolant of the previous time-steps level set values evaluated at the departure location, x d . The values at cell-centers are then smoothed using a semi-implicit technique first introduced for standard level set advection [66]: where β is a user-parameter, typically chosen to be 0.5 [64,66]. Note that the grid in this smoothing step is the cell-center (offset) grid. The boundary condition is the continuity of derivatives using the tentative level set values. For example, on the right edge the boundary condition is where N x and h x are the the number of grid points and the grid-spacing in the xdirection, respectively. Finally, the cell-center values are projected back onto the grid points by using second-order stencils. For example, the updated level set values at grid point (i, j, k) are given by while the gradient in the x-direction is given by The other gradient values can be similarly evaluated. The advantage of this method over that of Ref. [64] is that only a single smoothing operation is required, versus one for the level set and three for the gradient fields. Work is currently underway investigating the full properties of this semi-implicit jet scheme.

Level Set Reinitialization and Closest Point
Calculation. After the advection of the level set a reinitialization procedure is used. This is done for two reasons. First, it has been shown that reinitialization aids in the conservation of mass in level set-based simulations [67]. Second, the point on the interface with the shortest distance to a given grid point, called the closest-point, is required for the evaluation of the surface equations given below. The method employed here uses a variant of Newton's method and a tricubic interpolation to find the closest point on the interface from an arbitrary grid location in the neighborhood, see Ref. [68] for more details. Once the closest points are determined the level set function is replaced by the corresponding signed distance value while the gradient field is set to the unit outward normal vector, which is the normalized vector pointing from the closest point to the grid point, corrected so that it points outward. This procedure does not need to be done everywhere in the domain, but only near the interface. In the results below all nodes within five grid points of the interface have their closest point calculated and explicitly reinitialized. The remaining nodes are reinitialized using a first-order PDE based reinitialization scheme [69,70].

Hydrodynamics Solution: Projection Method.
There are four total conservation conditions that must be satisfied during the course of the simulation: 1) local surface area, 2) total surface area, 3) local fluid volume, 4) total fluid volume. Vesicle dynamics are extremely sensitive to any changes in these quantities [21,22,26,27] and thus high accuracy is required. Unfortunately, error in the solution of the fluid equations, in addition to the non-conservation properties of the level set method, can induce large errors over the course of the simulation. To account for this, the fluid equations will be modified to explicitly correct any error accumulation. The idea presented here is based on the work of Laadhari et. al. [24]. The difference is that in the previous work the corrections were not reflected in the fluid field, but instead modified the velocity used to advect the level set. Here, the corrections are included in the fluid field calculation and therefore a modified advection field is not required. Preliminary results for standard multiphase flow problems were presented in Ref. [71] and are extended here for vesicles.
A projection method is implemented to solve for the velocity, pressure and tension. First, a semi-implicit update is performed to obtain a tentative velocity field, u * − u n ∆t + u n · ∇u * = −∇p n + 1 Re ∇ · ∇u * + (∇u n ) T + f n H + f n γ + f n el , (3.11) where f H , f γ and f el are the bending, tension, and electric field forces localized around the interface, Eq. (2.22). The superscript n refers to the solution at the previous time step. Next, the tentative velocity field is projected onto the divergence and surfacedivergence free velocity space, where q and ξ are the corrections needed for the pressure and tension, respectively. Finally, the pressure and tension are updated by including the corrections,

13)
γ n+1 = γ n + ξ. (3.14) The four conservation conditions can be written as [24] ∇ · u n+1 = 0 (local volume conservation), The use of only the pressure and tension is not sufficient to satisfy all four conservation conditions. There, the pressure and tension fields are split into a constant and spatially-varying component: wherep andγ are spatially varying while p 0 and γ 0 are constant. Note thatp,γ, p 0 , and γ 0 all vary in time. Conceptually, this splitting allows for the enforcement of local conservation throughp andγ while global conservation is enforced through p 0 and γ 0 . The corresponding corrections are now q =q + (1 − H(φ))q 0 , and ξ =ξ + ξ 0 , while the projection step, Eq. (3.12), is now written as Noting that the time derivatives of the volume and area are to correct any accumulated errors in the solution, and using Eq. (3.21), the four conservation equations can be written in terms of the four unknowns (q,ξ, q 0 , and ξ 0 ), the current area and volume, and the initial area and volume. Specifically, applying the local area conservation equation requires that −∇·u * = ∆t∇· −∇q + δ(φ)q 0 ∇φ + δ(φ) ∇φ ∇ sξ −ξH∇φ + ξ 0 H∇φ , (3.22) while the total volume conservation requires that where V n is the current volume and the time-derivative of the volume is chosen so that at the end of the time-step the volume equals the initial volume, V 0 . Conservation of local and global area results in the following two equations:

24) and
where A 0 is the initial surface area and A n is the current surface area. Letq represent the vector holding the values ofq in the entire, discretized domain andξ is the vector holding the values ofξ in the entire, discretized domain. It is now possible to represent the conservation relationships, Eqs. (3.22)-(3.25), as linear operators acting onq, q 0 ,ξ, and ξ 0 . Local volume conservation, Eq. (3.22), can be represented as where D is the discrete divergence operator and where E V ≈ (V 0 − V n )/∆t − Γ (n · u * ) dA is the discrete form of the global volume correction needed. The surface-conservation equations are written in a similar manner: where D s is the discrete surface-divergence operator and for the local surface-area conservation, Eq. (3.24) and where E A ≈ (A 0 − A n )/∆t − Γ (Hn · u * ) dA is the discrete form of the global area correction needed and Using the notation above the set of four linear equations can be written in matrixvector form: Note that the specific spatial discretization have not yet been specified, as the general concepts are discretization-independent. While it is possible to form a globally-assembled matrix for Eq. (3.38), it would be computationally expensive. Not only would the size of the system be large, but it would also have to be re-formed every time step as many of the components, such as L s ξ , l s ξ and s T H depend on the location of the interface, which changes over time. Instead, write the complete system in the following simplified notation, where S is the Schur Complement of the original partitioned matrix and is give by The application of the Schur Complement is accomplished by a matrix-free iterative solver, such as GMRES.
Let the focus now turn to the calculation of D −1 . This inverse is given implicitly through the solution of the following, generalized linear system: Turning to the Schur Decomposition the solution can be written as where the Schur Complement is a rank-1 update on the matrix L s ξ : Using the Sherman-Morrison formula the inverse of this Schur Complement is given by So long as b ξ = s T H L s ξ −1 l s ξ this inverse is defined. While this inequality will not be proven here, a check during the course of the simulations presented below is performed every time step and this condition has never been violated.

Evaluation and Advection of Surface Quantities.
A word must be said about the solution of surface equations. In the projection method the pressure is defined everywhere in the domain and thus standard techniques can be used. The tension, on the other hand, is only defined on the interface and must be handled using special methods. Unlike the immersed boundary method or other front-tracking techniques Lagrangian points are not tracked over the course of the simulation. Instead, all surface equations and quantities, such as the tension and trans-membrane potential, are evaluated using the Closest Point Method [72,73]. The general idea is to replace surface equations by localized equations in the embedding space by extending quantities from the interface in a systematic way. This is done by replacing the value at a grid point in a finite difference approximation by it's interpolated value at the closest point on the interface. Additional information can be found in References [72,73]. All surface quantities, such as the tension and trans-membrane potential, are advected along side the level set. After the level set is advected surface quantities are advected by using the total-derivative form of the advection equation for a material quantity f : where I 3 (f n , x d ) is the evaluation of the cubic interpolant using the previous timesteps value evaluated at the departure location x d . The values are then replaced by the values at the closest point to a grid point, where x cp is the closest point on the interface to grid point x.
3.6. Overall algorithm. The overall algorithm to solve the electrohydrodynamics of the vesicle is performed using a time-staggered approach: Step I: Solve for the updated electric field and trans-membrane potential using Sec. 3.1.
Step II: Solve the hydrodynamics equations and obtain the updated velocity field, pressure and the tension using Sec. 3.4.
Step III: Using the updated velocity field from Step II advance the level set (interface) using Sec. 3.2.
Step IV: Reinitialize the level set values around the interface and recalculate the closest point data, Sec. 3.3.
Step V: Advect surface quantities such as tension and trans-membrane potential using Sec. 3.5.
Step VI: One time step has been completed. Return to Step I and repeat until final time is reached.

Results.
In this section the dynamics of the vesicle is studied in pure shear flow, in the presence of an electric field, and in combination of the two effects. The vesicle surface area is fixed to 4π in all the situations while the enclosed volume is varied. For all results presented here a collocated, Cartesian mesh with uniform grid spacing in each direction is used. Periodicity is assumed in the x-and z-directions while wall boundary conditions are given in the y-direction. Unless otherwise stated the domain is a box covering the domain [−4.5, 4.5] 3 . Any externally applied shear flow is given by imposing a velocity of u bc = (χy, 0, 0) on the wall boundaries, where χ is the normalized shear rate. Electric fields are applied by fixing the electric potential in the y-direction at E bc = E 0 y, where E 0 is the normalized electric field strength. It is assumed that the trans-membrane conductance is zero and thusĜ m = 0. Finally, compact finite difference approximations for all spatial derivatives have been implemented (e.g. the pressure-Poisson term is approximated by ∇ · ∇q ≈ ∇ 2 q at the discrete level). For all results surface area and volume are all conserved to within 0.01% error. Here two major modes of motions for the vesicle in shear flow are reported: tanktreading and tumbling. The tank-treading motion happens if the viscosity ratio, η, is less than a critical value η c , which depends on the reduced volume of the vesicle [18]. In the tank-treading regime the vesicle reaches an equilibrium angle and stays at that position. When the viscosity ratio is above the critical viscosity the behavior of the vesicle changes, and it starts to tumble end-over-end. Studies have shown that the dynamics depend primarily on the viscosity ratio, η and the reduced volume, v, with no remarkable dependence on the shear rate strength, χ [17].
Sample tank-treading results are illustrated in Fig. 4.1(a) for two initial conditions: a disk-like vesicle and an ellipsoidal vesicle. Both initial conditions have the same reduced volume of v = 0.85 and the dimensionless shear rate at the boundary is set to χ = 1. The inclination angle with respect to the shear-flow axis, θ, the semimajor axis L, semi-minor axis B, and the half axis length in the vorticity direction, Z, are reported. These values are computed through the use of the inertia matrix of the vesicle [74]. The eigenvalues of the interia matrix correspond to the axis lengths, while the angle between the eigenvector associated with the largest eigenvalue and the y−axis is the inclination angle.
It is clear from Fig. 4.1(a) that both initial conditions result in the same equilibrium shape and inclination angle. Figure 4.1(b) shows the vortices being formed in the interior of the vesicle at the equilibrium state. The streamlines demonstrate that the fluid velocity is tangent to the membrane which indicates the tank-treading behavior.
The equilibrium inclination angle as a function of reduced volume and viscosity ratio is shown in Fig 4.2(a). A dimensionless shear rate of χ = 1 is applied at the boundary and results are compared to experimental data from Ref. [75]. The angles from the numerical model are in good agreement with experimental measurements. Small discrepancies in the case of η = 4.9 could be related to the role of thermal fluctuations in the dynamics of vesicles with larger viscosity ratios [75]. This effect is absent in the proposed model and further investigation is needed to fully understand this phenomena.
The transition between tank-treading and tumbling happens at the critical viscosity ratio η c . As v increases a larger η c is required to transition from tank-treading to tumbling. This behavior is shown in Fig 4.2(b). The present method compares well to experimental data from [14], particularly for larger reduced volumes.

4.2.
Vesicle dynamics in strong DC electric field. The response of a vesicle to an external electric field will depend on a number of factors, including the membrane properties, the viscosity ratio, the electrical properties of the fluids, and the applied electric field. The shape of the vesicles is described by the deformation parameter: where a y is the dimension of the vesicle in the direction parallel to the applied electric field (the y-direction) while a x is the dimension of the vesicle in the direction per- pendicular to the electric field, the x-direction. Note that these parameters are being calculated directly using the interface location and differ from the semi-axis values defined in Sec. 4.1.
If the conductivity of the inner fluid is lower than the outer fluid, λ < 1, and the electric field is strong enough such that t ehd < t m , an initially prolate vesicle will undergo a prolate-oblate-prolate (POP) transition [40,43]. An oblate shape is given by D < 1 while a prolate shape has D > 1.
As verification of the method simulation results are compared to the analytic work of Schwalbe et. al. [43]. The reduced volume of the vesicle is v = 0.98. The initial shape is given by a second-order Spherical Harmonics parametric surface, see Ref. [43] for details. In Figs. 4.3-4.5 the influence of the time step, grid size, and domain size is explored and results are compared to Ref. [43]. Note that the membrane charging time is chosen as the characteristic time in these simulations. Recall that using the typical experimental values of the physical parameters discussed earlier in Sec.2.5 this results in t 0 = t m = 2.1 × 10 −3 s,Ĉ m = 0.095,Ĝ m = 0, Re = 0.19, Ca = 3.8 × 10 4 , M n = 18, E 0 = 1, and χ = 0. Also a matched dielectric ratio is used in all the electrohydrodynamic simulations in this paper.
Several points need to be made about the results. First, the use of a time step of size ∆t = 5 × 10 −4 , grid size of N = 120 3 with a domain of [−4.5, 4.5] 3 results in a converged solution and is therefore used in the following simulations. Second, there is a discrepancy between the simulation results and the analytic results of Schwalbe et. al. While the general results are similar, for example the approximate time at which the prolate-oblate transition occurs and the time in the oblate shape, there are important differences. The transition from the initial prolate to oblate shape is delayed and is sharper for the current simulation versus the analytic result. Also, the deformation parameter increases slightly during the oblate shape (at a time of 0.4t m ) before decreasing again. Finally, the final equilibrium shape is slightly different between the two methods, as demonstrated by the difference in the final deformation  [43]. Convergence test is done for three different grid spacing to check the grid independence. The domain size is [−4.5, 4.5] and the time step is set to ∆t = 5 × 10 −4 . The investigation reveals that that h = 0.075 is the maximum acceptable grid spacing.
parameter. This difference is most likely due to the use of a mode-2 spherical harmonic description of the interface in the analytic work. The current simulation has access to a wider range of vesicle shapes than the analytic work and thus additional deformation modes will be accessible. This can be seen by comparing the cross-section of the vesicle in the x − y plane for the two results, Fig. 4.6. The analytic works always remains ellipsoidal in shape while the numerical simulation developed here demonstrates the availability of additional shapes. The other possible explanation for the discrepancy between the two results would be related to the use of Stokes assumption in the analytic solution. The present method uses a full Navier-Stokes solver and therefore possible inertia effects are taken into account.
To demonstrate the axisymmetric POP behavior an initially prolate vesicle with reduced volume of v = 0.93 and matched viscosity (η = 1) is placed in a strong electric    trans-membrane potential on the interface is also shown, which demonstrates the charging process.
The results clearly demonstrate the transition from a prolate shape to an oblate shape, and then back up to prolate. The initial prolate-oblate transition occurs at about the membrane charging time. This type of behavior has been previously predicted numerically [43,47] and verified experimentally [40]. The cylindrical deformations observed in two-dimensional vesicle electrohydrodynamic simulations is also observed [47]. It is worth noting that the dimple at the center of the vesicle seen at times 0.35 t m and 1.40 t m are more pronounced than in the two-dimensional simulations. This could be due to the different parameters used here and the fact that three-dimensional vesicles are typically more flexible due to the additional curvatures and deformation modes.
The area and volume were also tracked during this symmetric POP simulation and the errors, as a percentage of the initial area and volume, is given in Fig. 4.8. As can be seen both the area and volume are conserved extremely well.
The methods developed here allow for differing viscosity between the inner and outer fluids. As has been demonstrated previously an increase in the viscosity ratio can dramatically change vesicle dynamics in the absence of electric fields. Here the result seen in Fig. 4.7 is compared to a vesicle with the same parameters except that the viscosity ratio is increased to η = 2. The resulting deformation parameter evolution is presented in Fig. 4.9. The increase in the viscosity ratio results in a vesicle which no longer undergoes a full prolate-oblate-prolate transition. Instead, the higher viscosity vesicle has a slight flattening, but never reaches an oblate shape. In both cases the final equilibrium shape is the same.  Fig. 4.7. This is compared against a vesicle with the same conditions but a viscosity ratio of η = 2.0. While a full POP transition is observed for the case with matched viscosity ratio, the vesicle with larger η never evolves into the oblate shape.

4.3.
Vesicle dynamics in combination of shear flow and DC field. When both a shear flow and DC electric field are applied the forces will compete to determine the dynamics of the vesicle. Previous works have demonstrated that electric fields damp the tank-treading and tumbling motion of vesicles in shear flow [43,47]. If the electric field and shear flow directions are perpendicular the vesicle will have forces which act in perpendicular directions, leading to a competition between the two forces.
To demonstrate this consider a vesicle with matched viscosity, η = 1 and a reduced volume of v = 0.93 in a normalized shear flow of strength χ = 1. The characteristic time is taken to be the shear flow timescale,t 0 = tγ = 1 s, and thus using standard membrane parameters the dimensionless coefficients becomeĈm = 2×10 −4 ,Ĝ m = 0, Re = 0.001, and Ca = 10. The behavior of this vesicle under the influence of three electric field strengths, M n = 0, M n = 3 and M n = 15 are shown in Fig. 4.10. Note that due to the membrane charging time, t m ≈ 2.1 × 10 −3 s, being much faster than the characteristic time it is assumed that the vesicle membrane begins fully charged. During each time step the trans-membrane potential is updated using a sub-step iteration until a pseudo-steady state is reached. As expected the application of an electric field results in the equilibrium inclination angle increasing. A comparison in the x − y plane for the vesicles with M n = 0 and M n = 15 are seen in Fig. 4.11, which clearly demonstrates the higher inclination angle.  Next consider the application of an electric field to a vesicle in the tumbling regime, η = 10 > η c . The other parameters are the same as the tank-treading case: v = 0.93, χ = 1,Ĉm = 2 × 10 −4 ,Ĝ m = 0, Re = 0.001, and Ca = 10. As in the tanktreading case the trans-membrane potential is iterated until a pseudo-steady state is reached every time step.
In the absence of an electric field the vesicle will tumble end-over-end, with the inclination angle undergoing periodic repetition, Fig. 4.12. The application of a weak electric field, M n = 3 results in periodic behavior more akin to trembling. The vesicle does not undergo a rigid-body-like rotation, but instead the vesicle poles retract and the vesicle reaches a nearly-spherical shape, see Fig. 4.13 for the threedimensional representation of the vesicle over time. A comparison of the vesicles in the x − y plane at the elevated viscosity ratio for the cases of M n = 0 and M n = 3 is shown in Fig. 4.14. This figure clearly shows the the tumbling behavior is due to membrane deformation and not and end-over-end rotation. As the electric field strength is increases to M n = 15 the vesicle no longer undergoes a tumbling behavior. Instead the vesicle reaches an equilibrium, tank-treading inclination angle. Compared to the standard tumbling case, the simulation with M n = 3 shows a lagged tumbling behavior and experiences a lot more topological changes during the transition as it will be illustrated below. An interesting observation has been made for the vesicle with M n = 15. In this situation, the vesicle undergoes a tank-treading motion and stays at that position permanently.

Conclusion.
In this work a numerical model of the electrohydrodynamics of three-dimensional vesicles has been presented. The work uses a semi-implicit version of the Jet scheme to implicility track the interface. This is coupled to a full Navier-Stokes projection method with explicit corrections to ensure volume and surface area conservation. The electric field contributions are computed through an Immersed Interface Method solver which takes into account the jump in the membrane voltage across the interface.
This work was verified by comparing the hydrodynamic results (in the absence of electric fields) to previously published experimental works. The influence of the time-  blue in the online version). The dynamics is compared against a tumbling vesicle with the same viscosity ratio and no external field, M n = 0 (dashed line; red in the online version). The rest of the parameters are Re = 0.001, Ca = 10 and χ = 1. The presence of the electric field affects the inclination angle and causes the vesicle to reach the equilibrium condition at a smaller θ. This behavior is due to the vertical alignment of the electric field in the y direction. The electric field force induces a resistance against any angular movement about the initial vertical position of the vesicle.
step, grid-spacing, and domain size on the electrohydrodynamics of vesicles was investigated. The simulation was capable of capturing the prolate-oblate-prolate dynamics observed experimentally and predicted by previous analytic and two-dimensional numeric works. The combined effects of fluid flow and electric fields were also presented, with weak electric fields dramatically changing the behavior of vesicles in shear flow.
This work is a step towards understanding how to control vesicle dynamics through the use of electric fields and fluid flow. Future work will use the methods developed here to perform a detailed investigation on how electric fields, fluid flow, and material parameters can be used to influence the dynamics of vesicles.