Finite element methods for director fields on flexible surfaces

DOI 10.4171/IFB/281

Finite element methods for director fields on flexible surfaces

S ¨ORENBARTELS

Universit¨at Bonn, Institut f¨ur Numerische Simulation, Wegelerstr. 6, 53115 Bonn, Germany

E-mail:bartels@ins.uni-bonn.de

GEORGDOLZMANN

Fakult¨at f¨ur Mathematik, Universit¨at Regensburg, 93040 Regensburg, Germany RICARDOH. NOCHETTO

Mathematics Department and Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, U.S.A.

ALEXANDERRAISCH

Universit¨at Bonn, Institut f¨ur Numerische Simulation, Wegelerstr. 6, 53115 Bonn, Germany [Received 14 June 2011 and in revised form 16 March 2012]

We introduce a nonlinear model for the evolution of biomembranes driven by the L²-gradient flow of a novel elasticity functional describing the interaction of a director field on a membrane with its curvature. In the linearized setting of a graph we present a practical finite element method (FEM), and prove a priori estimates. We derive the relaxation dynamics for the nonlinear model on closed surfaces and introduce a parametric FEM. We present numerical experiments for both linear and nonlinear models, which agree well with the expected behavior in simple situations and allow predictions beyond theory.

2010 Mathematics Subject Classification:Primary 35K55, 74K15.

Keywords:Finite elements, surfaces, biomembranes, surfactants, director fields

1. Introduction and discussion of the proposed model

The question of how to predict the shape of a cell bounded by a lipid bilayer membrane has inspired a significant body of research in the past twenty years ranging from purely mechanical descriptions to advanced mathematical analysis. We refer, e.g., to the papers [10,19,21,22] for the discussion of the shape of a red blood cell and the basic models developed for this purpose. Excellent reviews of the topic can be found in [27,30]. Almost all of these models share the basic structure given by an energy functionalE. /depending on the shape of the cell which is identified with a surface ,

E. /DZ 0

2 .H H0/²C G

2 K

d : (1.1)

HereH is the mean curvature andK the Gauss curvature of and0 andGare the associated moduli of elasticity. By Gauss-Bonnet, the integral ofKis a topological constant on a closed surface and can be neglected for evolutions in one topological class of surfaces. The quantityH0is usually referred to asspontaneous curvatureand describes the preferred value of curvature induced by the ambient space on a membrane in equilibrium. SuchH0 might be constant, the usual choice, but

c European Mathematical Society 2012

it might also depend on another variable such as the bilipid concentration [7,13,17,18,28,33].

Alternatively,H0might be induced by an underlying director field as in [26] and our models below.

The first variation of (1.1) with respect to is given in [4] forH0 constant and in [12] for H0

depending on position.

In a broader context, the question to find the shape of a cell is surprisingly similar to the related problem of determining the shape of an interface between two immiscible liquids with or without surfactants. The prediction of the structure and the elastic properties of such interfaces is still a challenging problem in applied mathematics and physics and has been investigated by a wide range of techniques reaching from molecular dynamics simulations to continuum descriptions in coarse grained models; see [26] that inspired this work. These similarities motivate us to explore model energies for membranes which combine classical elasticity terms like those present in (1.1) with terms which couple the local orientation of the surfactants or the lipid molecules with the curvature of the interface or the membrane, respectively; we started this investigation in [6]. These energy contributions are relevant in the gel phase of the membrane.

In this paper we investigate a novel model for the shape of a lipid bilayer membrane which takes into account a coupling between the curvatureH D div of the membrane and the local orientation of the lipid molecules, described by a director fieldn. The nonlinear model is governed by the energy

E. ; n/DZ 1 2

ˇˇdiv ıdiv nˇˇ²C 2

ˇˇr nˇˇ²

d; (1.2)

where div andr are the tangential divergence and gradient operators and > 0. Comparing with (1.1), we can interpretH0 D ıdiv n withı 2 Ras an induced spontaneous curvature on due to the coupling withn. In order to develop an effective approximation scheme, we first linearize this model locally in a flat region of and represent as a graph with heightu(Monge gauge). The resulting model is a special case of that introduced by Laradji and Mouritsen [26]. We study this model, propose a practical FEM for anL²-gradient flow of the linearized functional, and show a priori estimates, which lead to existence of a limiting solution pair.u; n/. We also explore the dynamics of defects using our FEM. We next return to the nonlinear functional (1.2), derive an L²-gradient flow, propose a practical FEM for its solution, and show simulations of defects. The insight gathered from the linearized graph case turns out to be useful in understanding the nonlinear regime.

1.1 A model for surfactants

The starting point of our analysis is the Ginzburg–Landau model in Laradji and Mouritsen [26]

which was originally developed for surfactant monolayers at liquid-liquid interfaces with a locally varying density of surfactants. The formulation assumes that this interface is given by a two- dimensional surface in the three-dimensional ambient space described by a height functionu W

˝ ! R. The model in [26], which is discussed below, is an attempt to match deviations from the bending energy model (1.1) with H0 D 0 for low wave numbers, which were detected via molecular dynamics computations. The total energy of a configuration is assumed to be given by

(see Appendix A in [26]) F.u; ; n/D

Z

˝

p

1C jruj²C

2jdivj²C a 2 ²C c

2jrj² s C g./

2 jnj² h./nC k./

2 jdivnj² `./

2 divdivn dx with suitable constants,, a,c,s and nonnegative functions,h, g, k, and `. Here r and div denote the planar differential operators gradient, i.e., rz D .@1z; @2z/ for a scalar function z, and divergence, i.e., divF D @1F1C@2F2 for a vectorfieldF D .F1; F2; F3/, whereas D . ru; 1/=p

1C jruj²is the normal to the graph ofu. In [26] it is shown that the surface tension is vanishingly small for densitiesclose to one. Therefore we may assume that0and thatis nearly equal to1and discard all terms depending onandinF. As a further simplification and in order to focus on the interaction of orientation and curvature, we assume thatnis a unit vector and we omit for the moment the term proportional tonwhich favors alignment ofnalong. This leads to the following model which contains the essential features

F.u; n/D Z

˝

2jdivj²C k

2jdivnj² `

2 divdivn dx;

with constant parameters; k; `. Upon completing the square, one obtains F.u; n/D

Z

˝

2

div `

2 divn2

Ck 2

`² 8

jdivnj² dx:

Comparing with (1.1) and (1.2), we interpret this formula as saying that the local arrangement of the surfactants leads to a (position dependent) spontaneous curvature

H0D ` 2 divn ;

which becomes less important for large values of the bending rigidity. We finally observe that in order to bound the energy it is sufficient to assume that

k 2

`² 8 >0 :

The corresponding positive convex termjdivnj²in the energy gives us coercivity of the functional F.u; n/but it is insufficient for devising a practical numerical scheme, deriving a priori bounds for discrete solutions which allow passing to the limit, and showing existence of a minimizing pair .u; n/. We thus modify the model upon replacingjdivnj²with the usual Frank energyjrnj²of the director fieldn, which is ubiquitous in the theory of liquid crystals. As the maximal mesh size tends to zero we may pass to the limit forn, in view of the enhancedH¹regularity, as well as to enforce the unit length constraint onnvia a projection method due to Alouges [1], and extended in [5] to FEM. Such a projection does not increase the energy of the Dirichlet integral, but the analogous assertion is not true for the energy of the divergence.

1.2 A model for biomembranes

Our interest in augmented Canham–Helfrich models originates in the search for models that allow one to predict the experimentally observed coarsening mechanisms in membranes in the gel phase based on recombination of topological defects [23]. Related models, based on the assumption that this recombination is driven by an interaction between the director field and the curvature, have been proposed in [32] and analyzed in [6]. See also [20] for a closely related approach.

In the model in [32] the lipid monolayer is considered in the gel phase and it is assumed that the director field is oriented in a fixed angle relative to the surface normal [29]. Therefore it suffices to study the tangential partmof the director field which is itself a vector field of fixed length. The related energy functional in a linearized setting is

E.u; m/D 2

Z

˝juj²dxCCq

2 Z

˝jrmj²dx ı Z

˝

D²uW.m˝m 1 2I /dx;

subject to a length constraint for m. Our numerical experiments for a rigidly imposed length constraint show that the coupling betweenuandmis too weak in the regime of parameters which define a well-posed minimization problem in order to simulate the observed recombination of defects [6]. The coupling proposed in the present model is stronger in the sense that it involves one more derivative. It also allows a direct extension to the nonlinear model (1.2) on closed surfaces;

cf. Section1.6.

1.3 Linear model on graphs

We give a linearized version of (1.2) for surfactants and augment the obtained energy by a term which penalizes deviations of the out of plane component from a given value to model biomembranes. We focus on the local situation in which the surface is described by the graph of a functionuW˝ !R³with˝ R²convex. Moreover, we assume that the displacements are small,

jruj 1:

This yieldsp

1C jruj²1as well as . ru; 1/, whence div u. Moreover, we have r n rn; div ndivnp;

wherenp stands for the tangential part of the director fieldnD.n1; n2; n3/, that is,np D.n1; n2/, andr, div are the planar differential operators. We are now ready to write the linearized energy:

findu 2 H².˝/withu D uD on@˝,uD 2 H².˝/,n 2 H¹.˝IR³/withn D nD on@˝, nD2H¹.˝IS²/and2L¹.˝/as stationary points of the integral

E.u; n; / D 1 2

Z

˝juCıdivnpj²dx C 1

2 Z

˝jrnj²dxC 1 2

Z

˝

jnj² 1 dx

Z

@˝

g@udS : (1.3) Note that 2 L¹.˝/is the Lagrange multiplier for the nonlinear constraintn 2 S² and thatg is related to the boundary values for the mixed method we discuss below. This model captures the essential features of the simplified linear model of Section1.1with energyF.u; n/.

To model surfactants we do not impose an angle betweenandn, which typically tend to align in the gel phase of the membrane. To model biomembranes, instead, we penalize the deviation of nn3from a prescribed value0via

1 2"²

Z

˝ jn3j² ₀²2

;

with small parameter" > 0. This term being lower order does not cause difficulties in the numerical method or the passage to the limit and will thus be ignored for the subsequent discussion until Section4.

1.4 Relaxation dynamics for surfactants

To detect critical points we suggest a relaxation dynamics given by anL²-gradient flow, i.e., we assume that there exist constants uand n> 0such that

h@tu; vi D ^uhıE

ıu; vi for allv2H².˝/\H₀¹.˝/ ; h@tn; mi D ⁿhıE

ın; mi for allm2H₀¹.˝;R³/ :

For simplicity we assume in the following that the units are chosen in such a way that uD ⁿD1.

If we include the equilibrium condition for the Lagrange multiplier in our equations, then we obtain the following coupled system of partial differential equations: for allv2H².˝/\H₀¹.˝/, for all m2H₀¹.˝IR³/, and for all2L¹.˝/,

h@tu; vi D uCıdivnp; v C

Z

@˝

g@vdS ; h@tn; mi D uCıdivnp; ıdivmp

.rn;rm/ .n; m/ ; 0 D 1

2 ;jnj² 1 :

(1.4)

Hereafter we write for simplicity .;/ for the inner product in L². We impose the following boundary conditions provided by the setting of the model,

uDuD; nDnD on@˝ ;

and we need to choose a second boundary condition for the fourth order equation involvingu. Such a condition is implicit in the equation forutabove because integration by parts gives formallyzQDg with

Q

zDuCıdivnp:

Despite the fact that this quantity is a priori only inL², we prove that there exists a solution with Q

z 2 H¹so that the boundary conditionzQ D g is well-posed. Note that this is a natural condition at first sight in the energy minimization but it becomes essential for the operator splitting: we use a mixed method for the variablesuandzD Qz gwith homogeneous Dirichlet boundary conditions.

Finally we collect the equations in their strong form:

Q

z DuCıdivnp; zQˇˇ_@˝ Dg ; (1.5)

@tu D z ;Q uˇˇ

@˝ DuD; (1.6)

@tnp Dır QzCnp np; np

ˇˇ

@˝ DnD;p; (1.7)

@tn3 Dn3 n3; n3

ˇˇ

@˝ DnD;3; (1.8)

jnj² 1 D0 ; a.e. in˝ : (1.9)

The essential difference with respect to the model proposed by Uchida [32], and analyzed in [6] for a rigid constraintjnj D 1, is the additional derivative of divnp in the coupling termz. This leadsQ to additional difficulties in the stability analysis of the numerical scheme as compared to [6]. We propose in Section3.3a semi-implicit algorithm for the computation of approximate solutions in finite element spaces and prove uniform bounds for a suitable energy of the system. We then present in Section4several numerical experiments displaying quite interesting dynamics of defects. Note, that theL²-gradient flow for biomembranes is obtained by adding the term_"¹2.jn3j² ₀²/n3to the equation (1.8).

1.5 Qualitative analysis of defect-shape interaction

In order to understand the interaction of defects and shape in the biomembrane case, i.e., when the angle between the director and surface normal is fixed, we consider in Sections2 and 4 a decomposition of the director field n into a tangential and normal part. The normal part is a fixed multiple of the surface normal and the tangential, planar part np has a fixed length. This decomposition allows us to construct in Section2formal stationary solutions with uDdivnp. The proposed director fields are (infinite energy) limits of energy-minimizing configurations for a Ginzburg–Landau regularization of the Frank energyR

˝jrnpj²dxsubject to their own boundary data, cf. [8]. This approach allows a precise characterization of the shape corresponding to different defects and provides insight on the long time asymptotics of.u; n/. In the numerical experiments for the linear model on graphs reported in Section4we allow the tangential part of the director field to develop an out-of-plane component, so that the full director field violates the angle condition and finite energy minimizers are possible. We observe that for defects of degree˙1the asymptotic behavior is dictated by the solutions found in Section2. It is important to realize that, in contrast to [6], our new model with rigid constraintjnj D1admits defects in the limit becausenpis allowed to go out of plane near point singularities, a feature fully documented in Section4. The numerical results for the full model on closed surfaces reported in Section 6show that the theoretical and practical predictions of the interaction of defects with the membrane shape in the simplified case explain the interesting dynamics occurring in the full biomembrane model for which the presence of defects is unavoidable if the angle betweenandnis fixed.

1.6 Nonlinear model on closed surfaces

For a smooth embedded surface R³, a director fieldnW !S²and constantsı; "and, we consider the energy (1.2) augmented as follows

E. ; n/D 1 2

Z

jdiv ıdiv nj²dC 2

Z

jr nj²d C1

2 Z

.jnj² 1/dC 1 2"²

Z

f .n/d;

whereis the Lagrange multiplier for the rigid constraint jnj D 1andf is given byf .x/ D .x² ₀²/², for0 2Œ0; 1. The last term penalizes the deviation of the three-dimensional director fieldnfrom the cone of all vectors that have a given angle with respect to the unit normal to the surface, as discussed already in Section1.3. Thus, for"D 1, which corresponds to neglecting the last term, we obtain the surfactant case, while"1results in the modelling of biomembranes. In Section5, we derive a variation of the energy with respect to andn, which is the first step on the way to discover critical points ofE. ; n/. We also introduce a semi-implicit algorithm based on parametric finite elements of Barrett, Garcke and N¨unberg [3] to model theL²-gradient flow ofE. ; n/, see also [2]. As we are interested in the simulation of cells and biomembranes, side conditions like conservation of the enclosed volume and/or the surface area are important. For this purpose we use a Newton-iteration method, as proposed in [9]. In Section6we explore the behavior of the nonlinear model via simulations. We first show that forı D1and without angle penalization, the vectorsandntend to align since this minimizes.div . n//². We also display the dynamics of defects of degree˙1and observe that locally the membrane shape is similar to that discovered earlier in the graph case. We conclude that defects of the director fieldnhave a dramatic effect on the shape of , as observed in experiments, e.g., reported in [24].

2. Qualitative Behavior of Graphs

To build intuition about the mechanisms introduced by the coupling term in the model, we fix stationary tangential director fieldsnpof unit length and compute a functionu2H₀¹.˝/for which the first term in (1.2) vanishes. For ease of readability we omit the subscriptp throughout this section. We thus impose that the auxiliary variable

Q

zDuCdivn vanishes, thereby giving the relation

uDdivn:

Motivated by experimental observations we are particularly concerned with the surface structure when the director field represents a defect of positive or negative degree-one, i.e.,

nDexp.˙i /Dcos˙isin

in polar coordinates .r; / and complex notation. Notice that for such a field n we have R

˝jrnj²dx D 1, soncannot be a minimizer of the energy which involves the Dirichlet integral ofn, but it arises as the limit of minimizers of a corresponding Ginzburg-Landau regularization that penalizes the unit-length constraint, cf. [8]. Therefore our calculations are only meant to explain the

structures observed in our experiments of Section4which necessarily involve regularizations of the corresponding fields. We first compute the divergence

divnD@xcos˙@ysin D sin @x˙cos @y; and recall that Darctany=x, whence

@x D y

x²Cy² ; @y D x x²Cy² : We insert this result into the expression for divnand obtain

divnD y²˙x²

r³ D sin²˙cos²

r :

2.1 Positive degree-one defects

We now takenDexp.i /DcosCisin. We thus seekusuch that the inhomogeneous equation (in polar coordinates)

uD 1

r@r r@ru C 1

r²@²uD 1 r

holds. It is natural to look for a radial solutionu.r/D r^˛and the expression foruimplies the necessary condition˛D1and the cone-like surface (see Figure1, top):

u.r/D r:

Consider now the director fieldnDe^i.C=2/rotated by an angle=2. Such annsatisfies divnD 0, whenceuD0; this is depicted in the lower plot of Figure1. Any other rotationnDe^i.C0/by an angle0can be expressed asnDcos0n1Csin0n2withn1; n2the director fields in Figure1.

The corresponding solution is thus

uD rcos0: 2.2 Negative degree-one defects

We now takenDexp. i /Dcos isin. We thus seekuas a solution of the inhomogeneous equation

uD 1 r@r

r@ru

C 1

r²@²uD sin² cos²

r D cos.2 /

r :

We try a solution of the formu.r; /D C r^˛cos.2 /for suitable constantsC; ˛ and evaluate the partial differential equation to obtain the necessary condition

uDC ˛² 4

r^{˛ 2}cos.2 /D cos.2 / r ; whence˛D1; C D1=3and

u.r; /D 1

3rcos.2 /:

−1 0 1

−1 0 1

−1 0

1

−1 0 1 0 1

−1 0 1

−1 0 1

−1 0

1

−1 0 1 0 1

FIG. 1. Positive degree-one defect: director fieldnDeⁱand cone-like surfaceuD1 r(top), and rotated director field nDe^i.C=2/and functionuD0(bottom) related by (finite element discretizations of) uDdivnwithuj@˝D0.

This solution is a saddle and is depicted in Figure 2(top). Consider now the director field n D e ^{i. 0/}which can be written asn D e ^{i. 0=2/}eⁱ⁰⁼². We thus realize that the value ofnat results from reading the value at 0=2and rotating clockwise by0=2, an effective rotation of e ⁱ by the angle0=2. The corresponding solution thus reads

u.r; /D 1

3rcos.2 0/:

Figure2(bottom) displays such a pair.u; n/for0D=2.

3. A semi-implicit scheme for graphs

For simplicity we suppress in this section the indexhin connection with all finite element spaces and functions, that is, we write, e.g.,T,Vand.u; n/instead ofT_h,V_hand.u_h; n_h/, respectively.

We use upper indices for the functions at discrete time steps. In particularn⁰ 2 Vis a suitable approximation of the initial datanD. We point out that the definition and the analysis of theL²- gradient flow presented here is for surfactants. The linearized model for biomembranes, which will be experimentally investigated in Section4, is obtained by adding a lower order term that can be neglected in the analysis.

3.1 Finite element spaces

We letT be a regular triangulation [11] of˝ into triangles of maximal diameterh > 0. We denote byV D V.T/the space of all continuous functions on˝ that are affine on the elements in the triangulationT and we setV₀ D V\H₀¹.˝/. We say thatT isweakly acuteif the sum of every

−1 0 1

−1 0 1

−1 0

1

−1 0 1

−0.1 0 0.1

−1 0 1

−1 0 1

−1 0

1

−1 0 1

−0.1 0 0.1

FIG. 2. Negative degree-one defect: Director fieldnDe ⁱ and saddle-like surfaceu.r; / ¹₃rcos.2 /(top), and rotated director fieldn D e ^{i. =2/}and corresponding rotated saddle-like surfaceu.r; / ¹₃rcos.2 =2/

(bottom) related by (finite element discretizations of) u Ddivnanduj@˝ D 0. Due to the boundary conditionu mimics the exact saddle structure only in a neighborhood of the origin.

pair of angles opposite to an interior edge is bounded by and if the angle opposite to every edge on the boundary is less than or equal to=2. Let.'a/a2Ndenote the standard nodal basis ofV. For later purposes, we note that ifT is weakly acute then [5]

Ki;j WD Z

˝r'ai r'aj dx60 for allai ¤aj 2N; (3.1)

whereN D fa1; : : : ; aNg denotes the set of nodes in T. For completeness we include now a monotonicity estimate due to Bartels [5] for finite element methods, following the seminal work of Alouges [1].

LEMMA3.1 (monotonicity) LetT be weakly acute, and leten2V³be such thatjen.a/j>1for all a2N, and definen2V³by settingn.a/Den.a/=jen.a/jfor alla2N. Then

krnk6krenk: (3.2)

Proof. Let.'ai/ai2Ndenote the nodal basis ofV. Besides (3.1), the symmetric matrix.Ki;j/^N_i;jD1

satisfiesPN

jD1Ki;j D0owing toPN

jD1'aj D1. We observe the relations jjrnjj²D

XN

i;jD1

Ki;jn.ai/n.aj/

D 1 2

XN

i;jD1

Ki;jn.ai/ n.aj/ n.ai/ C 1

2 XN

i;jD1

Ki;jn.aj/ n.ai/ n.aj/

D 1 2

XN

i;jD1

Ki;j

ˇˇn.a_i/ n.aj/ˇˇ²:

The assertion is proved ifjn.ai/ n.aj/j² 6 jen.ai/ en.aj/j²for alli; j D 1; ; N. Hence, it suffices to showˇˇ_jaj^a _jbj^bˇˇ6ˇˇa bˇˇ, fora; b 2R³withjaj;jbj>1. This follows from the Lipschitz continuity ofS² W fx2R³W jxj>1g !S²; x 7!x=jxj.

For a fixed time-step size > 0lettj Djfor allj >0. Givenq2ŒV³we define the space of tangential updates with respect to the sphere for a given vector fieldqwithjq.a/j D1for alla2N by

FŒqD˚

r 2ŒV₀³Wr.a/q.a/D0 for alla2N : (3.3) Since we use time-independent boundary conditions we may assume that we are given approximationsn⁰ 2 ŒV³ andu⁰ 2 VofnD anduD withjn⁰.a/j D 1for alla 2 N. Moreover we replace the additional variablezQ D uCıdivnp, which has a Dirichlet boundary valueg, byz D Qz g, which has vanishing trace. Givenn⁰ 2 ŒV³ andu⁰ 2 V, we letz⁰ 2 V₀ be an approximation toz.0/defined as

.z⁰; y/D .g; y/ .ru⁰;ry/ ı.n_p⁰;ry/ for ally2V₀; (3.4) and observe that the right-hand side in this equality defines a continuous linear form onV. Since the L²inner product is a norm onV₀ (with zero Dirichlet conditions), existence of a unique solution z⁰follows from the Lax–Milgram lemma.

In the numerical analysis of our proposed scheme we will need to controlkrz⁰k. For this we assume for simplicity that u⁰ˇˇ

@˝ D 0. Then, we define the discrete Laplacian ⁰ with zero boundary values for a finite element functionvto be the unique element⁰v2V₀that satisfies

.⁰v; w/D .rv;rw/ for allw2V₀ and let˘⁰denote theL²projection ontoV₀. We then have that

z⁰D ˘⁰.g ıdivn_p⁰/C⁰u⁰ and

krz⁰k6kr

˘⁰.g ıdivn_p⁰/C⁰u⁰

k: (3.5)

The assumptionu⁰ˇˇ_@˝ D 0can be avoided by appropriately splittingg D u.0/ˇˇ_@˝ C Qg and replacinggbygQin the foregoing discussion.

3.2 Difference quotients

We use two definitions of difference quotients in time, i.e., for any sequence.u^j/and for a sequence .n^j/which is obtained by a post processing of a sequence.en^j/we write

dtu^j D 1

u^j u^{j 1}

; edtn^j D 1

en^j n^{j 1} : For allj 2N,j >1, we have

.dtu^j; u^j/ D 1 2

ku^jk² ku^{j 1}k²

C

2kdtu^jk²D 1

2 dtku^jk²C

2kdtu^jk²; (3.6) .edtn^j;en^j/ D 1

2

ken^jk² kn^{j 1}k² C

2kedtn^jk²D 1

2edtkn^jk²C

2kedtn^jk²: (3.7) We recall a useful estimate for the discrepancyen^j n^j, and include a short proof [6, Proposition 4.2].

LEMMA3.2 (Discrepancyken^j n^jk) For the sequencesfen^jg^j>0andfn^jg^j>0constructed by the numerical scheme in Section3.3below, that is, satisfying in particular the orthogonality condition .n^{j 1};nQ^j n^{j 1}/D 0for allj >0, where we setn ¹ D0, the following estimate holds for all j >1,

ken^j n^jk²60⁴kedtn^jk²kredtn^jk²: (3.8) A possible numerical value for0is0D25.

Proof. For all nodesa2N, we can write ˇˇen^j.a/ n^j.a/ˇˇD ˇˇ

ˇˇen^j.a/ en^j.a/

jen^j.a/j ˇˇ

ˇˇD jen^j.a/j 1 :

Sinceen^j.a/Dn^{j 1}.a/Cedtn^j.a/we obtain from orthogonality, the normalization ofn^{j 1}, and the estimate.1Cx²/¹⁼²61C¹₂x²that

ˇˇen^j.a/ n^j.a/ˇˇ6ˇˇn^{j 1}.a/j²C²jedtn^j.a/ˇˇ²1=2

16 ² 2

ˇˇedtn^j.a/ˇˇ²:

In view of standard estimates for mass-lumping we get the inequality jjen^j n^jjj²_L²_{.T /} 6 ⁴jjedtn^jjj_L⁴⁴_{.T /}. We take the sum over all triangles, recall the fixed boundary values forn^j and use the multiplicative interpolation estimate [25] to obtain

ken^j n^jk²_L²_.˝/6⁴kedtn^jk⁴_L⁴_.˝/625⁴kedtn^jk_L²²_.˝/kredtn^jk²_L²_.˝/: This estimate concludes the proof.

We finally recall a variant of discrete integration in time: for16j 6J we have dt.u^jv^j/D dtu^jv^j Cu^{j 1}dtv^j, whence for16`6J

u^Jv^J u^{` 1</sup>v<sup>` 1}D XJ

jD`

dt.u^jv^j/D XJ

jD`

dtu^jv^j C XJ

jD`

u^{j 1}dtv^j: (3.9)

3.3 Numerical scheme

We propose a semi-implicit method for (1.4) in which the computation of the director field is naturally decoupled from the calculation ofu^j andz^j. We seten⁰ D n⁰ and seek forj > 1and givenu^{j 1},en^{j 1},n^{j 1}functions

u^j 2V₀; z^j 2V₀; edtn^j 2FŒn^{j 1} such that

.z^j; y/C.ru^j;ry/Cı.en_p^{j 1};ry/ D .g; y/ for ally2V₀; (3.10) .dtu^j; v/ .rz^j;rv/ D.rg;rv/ for allv2V₀; (3.11) .edtn^j; m/ ı.rz^j; mp/C ren^j;rm

Dı.rg; mp/ for allm2FŒn^{j 1}; (3.12) whereen^j Dn^{j 1}Cedtn^j. Now set

n^j.a/D en^j.a/

jen^j.a/j D n^{j 1}.a/Cedtn^j.a/

jn^{j 1}.a/Cedtn^j.a/j; for alla2N:

We remark that the system (3.10)–(3.11) has the structure of a saddle-point problem, that is similar to a hybrid formulation of the bilaplacian (with penalty term), i.e., (3.10)–(3.11) can be rewritten as

.z^j; y/C.ru^j;ry/D .g; y/ ı.en_p^{j 1};ry/ for ally2V₀; .rz^j;rv/ ¹.u^j; v/D ¹.u^{j 1}; v/ .rg;rv/ for allv2V₀:

Owing to the essential boundary conditions imposed onz^j andu^j aP1 P1 discretization is stable. The Lax–Milgram lemma implies the unique solvability of (3.12) on the non-empty linear spaceFŒn^{j 1}.

3.4 Stability analysis

In this section we derive energy estimates for the solutions of the finite element discretization. In particular, we verify bounds (uniform in) for the following quantities:

A.J / D 1

4krz^Jk²C1 2

XJ

jD1

kdtru^jk²C 2

XJ

jD1

kdtrz^jk²;

B.J / D 1

2kz^Jk²C 1

2krn^Jk²C 1 2

XJ

jD1

kdtu^jk²C kedtn^jk²

C 2

XJ

jD1

kdtz^jk²C kedtrn^jk²

;

for allJ > 1. The quantities which are quadratic in do not provide uniform estimates for the solutions. However, they are needed in order to control various terms that appear on the right- hand sides of the following estimates. We state all bounds forj > 1and use the conventions that

Pk

jD`.: : : / D 0 and supjD`;:::;k.: : : / D 0 wheneverk < `. We seten ¹ D en⁰ D n⁰, whence dten⁰D0, and recall the definition (3.4) ofz⁰2V⁰. We assume

uD2H³.˝/; nD 2H².˝/; g2H².˝/: (3.13) Combined with (3.4), this implies thatkrz⁰kis uniformly bounded with respect toh.

LEMMA3.3 (first strong estimate) If (3.13) holds, then there exists a constant0such that forJ >1 A.J /D 1

4krz^Jk²C1 2

XJ

jD1

kdtru^jk²C 2

XJ

jD1

kdtrz^jk²

62₀²ı² sup

jD3;:::;J

kren^{j 2}k²C krn^{j 3}k²X^J

jD3

kedtn^{j 2}k²

Cı² XJ

jD2

kedtn_p^{j 1}k²C krz⁰k²C3

2krgk²: Proof. Using (3.4) and the discrete time derivative of (3.10), we obtain that forj >1

.dtz^j; y/C.dtru^j;ry/Cı.dten_p^{j 1};ry/D0 ; (3.14) becausegis time-independent. The choicesyD dtu^j in (3.14) andvD dtz^j in (3.11) yield

1

2 dtkrz^jk²C

2kdtrz^jk²C kdtru^jk² D ı.dten_p^{j 1};dtru^j/ .rg;dtrz^j/ ; (3.15) because of (3.6). Since (3.12) gives an estimate foredtn^j, it is thus natural to rewrite the right-hand side of (3.15) as follows:

ı.dten_p^{j 1};dtru^j/D ı.edtn_p^{j 1};dtru^j/ ı

.n_p^{j 2} en_p^{j 2};dtru^j/ 61

2kdtru^jk²Cı²kedtn_p^{j 1}k²C ı²

²ken_p^{j 2} n_p^{j 2}k²: Forj D1; 2, the last term vanishes; forj >3we use the estimate (3.8) to infer that

ken_p^{j 2} n_p^{j 2}k²62₀²² kren^{j 2}k²C krn^{j 3}k²

kedtn^{j 2}k²:

We take the sum in the foregoing estimates and multiply by. Using (3.9) and the fact thatg is time-independent, the second term on the right-hand side of (3.15) reduces to

XJ

jD1

.rg;dtrz^j/D .rz^J;rg/C.rz⁰;rg/6 1

4krz^Jk²C1

2krz⁰k²C3 2krgk²: This gives the asserted estimate.

LEMMA3.4 (second strong estimate) Suppose that (3.13) is valid andT is weakly acute, whence T satisfies (3.1). Then the following bound holds forJ >1:

B.J /6 1

2kz^Jk²C1

2kren^Jk²C 1 2

XJ

jD1

kdtu^jk²C kedtn^jk²

C 2

XJ

jD1

kdtz^jk²C kedtrn^jk²

6 1

2kz⁰k²C1

2krn⁰k²Cı²Tkrgk²C T 2kgk² Cı²

XJ

jD2

²kdtrz^jk²Cı²krz^Jk²

C2ı0¹⁼² max

jD1;:::;Jkrz^jk XJ

jD2

kedtn^{j 1}k²C²kedtrn^{j 1}k² :

Proof. We use (3.14) withy Dz^j, (3.11) withvD dtu^j, and (3.12) withmDedtn^j and employ ren^j;edtrn^j

D

2kredtn^jk²C 1 2

kren^jk² krn^{j 1}k²

which is a variant of (3.7), to arrive at 1

2 dtkz^jk²C

2kdtz^jk²C 1

2 kren^jk² krn^{j 1}k² C

2kedtrn^jk²C kedtn^jk²C kdtu^jk² D ı.dten_p^{j 1};rz^j/Cı.rz^j;edtn_p^j/C.dtru^j;rg/Cı.edtn_p^j;rg/ : (3.16) For the first two terms on the right-hand side we have forj >1

ı.dten_p^{j 1};rz^j/Cı.rz^j;edtn_p^j/ D ı

.en_p^{j 1} en_p^{j 2} en_p^j n_p^{j 1}

;rz^j/ Dı .d²_ten_p^j;rz^j/C ı

.en_p^{j 1} n_p^{j 1};rz^j/ :

(3.17)

We substitute (3.17) into (3.16), multiply the resulting expression by, and sum it fromj D 1to J. In view of (3.8) and the fact thaten_p⁰ Dn_p⁰, the second term on the right-hand side becomes

ı XJ

jD2

.en_p^{j 1} n_p^{j 1};rz^j/ 6ı XJ

jD2

ken_p^{j 1} n_p^{j 1}k krz^jk

6ı max

jD2;:::;Jkrz^jk XJ

jD2

0²kedtn^{j 1}k kedtrn^{j 1}k

6 ı

20¹⁼² max

jD2;:::;Jkrz^jk XJ

jD2

kedtn^{j 1}k²C²kedtrn^{j 1}k² : (3.18)

Instead, for the first term on the right-hand side we use the partial summation formula (3.9), to deduce

ı XJ

jD1

.d²_ten_p^j;rz^j/ Dı XJ

jD1

.dt.dten_p^j/;rz^j/

D ı XJ

jD1

.dten_p^{j 1};dtrz^j/ ı .dten_p⁰;rz⁰/Cı .dten_p^J;rz^J/

D ı XJ

jD2

.edtn_p^{j 1};dtrz^j/ ı XJ

jD2

.n_p^{j 2} en_p^{j 2}; dtrz^j/ Cı .edtn_p^J;rz^J/Cı.n_p^{J 1} en_p^{J 1};rz^J/ DICIICIIICIV:

(3.19) We now examine termsI–IVseparately. ForIandIIIwe simply use Cauchy–Schwarz to write

I 6 1 4

XJ

jD2

kedtn_p^{j 1}k²Cı²² XJ

jD2

kdtrz^jk²; III6

4kedtn_p^Jk²Cı²krz^Jk²: Sincekdtrz^jk62 max

16j6Jkrz^jk, we proceed as in (3.19) to find II6ı

XJ

jD2

kn_p^{j 2} en_p^{j 2}k kdtrz^jk

6ı0¹⁼² max

jD1;:::;Jkrz^jk XJ

jD1

kedtn^{j 2}k²C²kedtrn^{j 2}k and

IV6ıkn_p^{J 1} en_p^{J 1}k krz^Jk6 ı

20¹⁼²krz^Jk

kedtn^{J 1}k²C²kedtrn^{J 1}k : It remains to estimate the terms involvingg in the summation of (3.16). For the first term we use (3.13) and integrate by parts to get

XJ

jD1

.rg;rdtu^j/D XJ

jD1

.g;dtu^j/6 1 2

XJ

jD1

kdtu^jk²C T 2kgk²; whereT >J is the final time. For the second term we simply use Cauchy–Schwarz

ı XJ

jD1

.rg;edtn_p^j/6 1 4

XJ

jD1

kedtn^jk²Cı²Tkrgk²:

We finally combine the foregoing estimates and use the monotonicitykrn^jk 6 kren^jkforj D 1; 2; : : : ; J 1, established in (3.2), as well asen⁰Dn⁰, to obtain the asserted estimate.

THEOREM3.5 (a priori estimates) Let (3.13) hold andT be weakly acute. We define E.J /D 1

2kz^Jk²C1

2krn^Jk²; Z0D krz⁰k²C 3 2krgk²; G Dı²Tkrgk²C T

2kgk²; and set for arbitrary0 < " < 1

BDE.0/ CGC2^"; AD16₀²ı²B²C2ı²BCZ0: Suppose finally that

^{1 "}6 1

4ı²A; ^{1=2 "} 6 1 8ı0A¹⁼²B

: (3.20)

Then for allJ >1we have A.J / D 1

2krz^Jk²C 2

XJ

jD1

kdtrz^jk²C 1 2

XJ

jD1

kdtru^jk²6A

B.J / DE.J /C XJ

jD1

kdtu^jk²C 1

2kedtn^jk² C

2 XJ

jD1

kdtz^jk²C kedtrn^jk² 6B :

Proof. We proceed by induction. ForJ D1the estimates of Lemmas3.3and3.4imply A.1/6Z06A; B.1/6E.0/6B:

Suppose now that the assertion has been verified forJ 1>1, i.e., forj D1; : : : ; J 1we have A.j /6AandB.j /6B. Lemma3.3and the definition ofAimply

A.J /616₀²ı² sup

jD0;:::;J 1

B.j /B.J 2/C2ı²B.J 1/CZ06A : Similarly, Lemma3.4gives the upper bound

B.J /6E.0/CGC4ı²A.J /C8ı0¹⁼² max

jD1;:::;JA.j /¹⁼²B.J 1/ : The induction hypothesisB.J 1/6Band the conditions onlead to

B.J /6E.0/CGC2^"DB which proves the assertion of the theorem.

Remark 3.1 (regularity ofg). The H²-regularity of g, assumed in (3.13), can be weakened to g2H¹.˝/at the expense of a more technical proof of Theorem3.5. In fact, we avoid integration by parts in the first term involvinggin Lemma3.4and instead write

XJ

jD1

.rg;dtru^j/6 2

XJ

jD1

krdtu^jk²C 1 2

XJ

jD1

krgk²

for > 0arbitrarily small. The recursion forB.J /in the proof of Theorem3.5now becomes B.J /6E.0/CGCA.J /C4ı²A.J /C8ı0¹⁼²A.J /¹⁼²B.J 1/:

Replacing this bound into the recursion ofA.J /allows us to absorb the termA.J /forsufficiently small.

3.5 Weak solution of(1.4)

The a priori estimates of Theorem3.5allow us to establish the existence of a weak solution of the continuousL²flow that satisfies an energy inequality. We now state precisely the notion of solution already introduced in (1.4) and refer the reader to [6] for details about passing to the limit.

DEFINITION3.6 (weak solution) Let˝ R²be a bounded and convex Lipschitz domain and fix T > 0. We call a pair.u; n/a weak solution of (1.4) in the time intervalI D.0; T /if the following assertions are true:

(i) n2H¹.IIL².˝IR²//\L¹.IIH¹.˝IR²//, u2H¹.IIL².˝//\L¹.IIH².˝//;

(ii) jn.t; x/j D1for almost every.t; x/2I˝;

(iii) n.0;/DnD,u.0;/DuDwithuD2H³.˝/andnD2H².˝/;

(iv) n.t;/j^@˝DnDandu.t;/j^@˝ DuDin the sense of traces for almost everyt 2I;

(v) uCıdivnp 2L².IIH¹.˝//and satisfies for a.e.t 2I that.uCıdivnp/ˇˇ Dgwith g2H².˝/given;

(vi) for all.m; v/2 L².IIH₀¹.˝IR²//L² IIH².˝/\H₀¹.˝/

satisfyingmn D0almost everywhere inI ˝we have

Z

I

˚.@tu; v/C.uCıdivnp; v/ dt Z

@˝

g@vdS D0 ; Z

I

˚.@tn; m/C.uCıdivnp; ıdivmp/C.rn;rm/ dt D0 :

Remark3.2. We note that the formulation of Theorem3.5allows one to deduce the existence of a solution.u; n/that satisfies the energy inequality

1

2kuCıdivnpk²C1

2krnk²C Z T

0 k@tuk²C 1

2k@tnk² dt 6 1

2kuDCıdiv.nD/pk²C 1

2krnDk²Cı²Tkrgk²C T

2 kgk²: We refer the reader to [31] for related existence theories in the context of the harmonic map heat flow.

4. Numerical experiments for graphs

In this section we report on various numerical experiments for biomembranes carried out with the scheme devised and analyzed in the previous sections. Since we want to illustrate the interaction of defects and shape we consider the case of a membrane in the gel phase where the director field prefers to have a fixed angle with respect to the normal to the surface, say=2for convenience.

As in Section 2 the director field n has unit length but is allowed to develop an out-of-plane component to accommodate for topological defects; we omit the indexp throughout this section for the tangential partnp ofn. We thus augment the system of equations discussed in Section3.4 by the term" ².en^j₃; m3/, where the subscript3refers to the third, or out-of-plane, component of a vectorfield, in (3.12), i.e., for the evolution of the director field we employ the equation

.edtn^j; m/ ı.rz^j; mp/C ren^j;rm

C" ².en^j₃; m3/D0 : This modification corresponds to the additional penalty term

1 2"²

Z

˝jn3j²dx in the energy, i.e., our energy functional is

E.u; z; n/D 1 2

Z

˝jzj²dxC1 2

Z

˝jrnj²dxC 1 2"²

Z

˝jn3j²dx

subject to the relationz D uCıdiv.n1; n2/, the pointwise constraintjnj D 1, and Dirichlet boundary conditions foru, z, and n. We remark that the inclusion of an implicit treatment of the convex penalty term in the stability analysis for the numerical scheme in Section3poses no difficulties.

The goal of this section is to explore the qualitative behavior of the evolution for specific initial conditions with defects. Here, the terminology of a defect refers to a singularity in the renormalized planar part of the director field which is also called vortex. This evolution typically shows an initial phase with a significant change of the shape in order for the system to adjust to the given initial and boundary values which is followed by a slower evolution towards an equilibrium shape. In the figures we display typical intermediate shapes and states which are close to an equilibrium. In our simulations the domain˝and the parametersı,T, and"are given by

˝D. 1=2; 1=2/²; ıD1; T D1; "D10 ¹:

We denote by.r; /the usual polar coordinates inR²(with respect to the origin). The function'"

which is used in the extension of a function given on@˝to˝ is equal to'".r/D tanh.r="/. The initial values are always chosen to be

u0D0; gDI_hŒdivn0

for different choices of n0 and where I_h is the nodal interpolation operator. The sequence of triangulationsT<sub>`</sub>is generated by`uniform refinements (division of each triangle into four congruent ones) of the initial triangulationT₀ of˝ which consists of two triangles obtained by dividing˝ along the diagonalx1 D x2. Hence the mesh-sizeh` is given byh` Dp

22 <sup>`</sup>. Moreover we used ` Dh`=.8p

2/as time-step size.

4.1 Positive degree-one defect

We choose boundary conditions which correspond to a defect of degree one, i.e., n0j^@˝.x1; x2/Dn0j^@˝.rcos; rsin/D cos./;sin./; 0

D eⁱ ; 0

; .x1; x2/2@˝