Hydrodynamic interactions in polar liquid crystals on evolving surfaces
HHydrodynamic interactions in polar liquid crystals on evolving surfaces
Ingo Nitschke and Sebastian Reuther
Institute of Scientific Computing,Technische Universit¨at Dresden,Germany
Axel Voigt
Institute of Scientific Computing,Technische Universit¨at Dresden and Dresden Center for Computational MaterialsScience (DCMS) and Center for Systems Biology Dresden (CSBD), Dresden,Germany
We consider the derivation and numerical solution of the flow of passive and active polarliquid crystals, whose molecular orientation is subjected to a tangential anchoring on anevolving curved surface. The underlying passive model is a simplified surface Ericksen-Leslie model, which is derived as a thin-film limit of the corresponding three-dimensionalequations with appropriate boundary conditions. A finite element discretization is con-sidered and the effect of hydrodynamics on the interplay of topology, geometric prop-erties and defect dynamics is studied for this model on various stationary and evolvingsurfaces. Additionally, we consider an active model. We propose a surface formulationfor an active polar viscous gel and exemplarily demonstrate the effect of the underlyingcurvature on the location of topological defects on a torus.
PACS numbers: 61.30.Jf, 61.30.Hn, 47.50.Cd, 47.11.Fg
I. INTRODUCTION
Liquid crystals (LCs) are partially ordered materialsthat combine the fluidity of liquids with the orientationalorder of crystalline solids (Chaikin and Lubensky, 1995;de Gennes and Prost, 1993). Topological defects are a keyfeature of LCs if considered under external constraints.In particular on curved surfaces these defects are impor-tant and have been intensively studied on a sphere (Bates et al. , 2010; Dhakal et al. , 2012; Dzubiella et al. , 2000;Koning et al. , 2013; Napoli and Vergori, 2012; Shin et al. ,2008) and under more complicated constraints (Martinez et al. , 2014; Prinsen and van der Schoot, 2003; Stark,2001). LCs on curved surfaces can be realized on vari-ous levels. One possibility is to prepare a double emul-sion of two concentric droplets (Fern´andez-Nieves et al. ,2007) for which the intervening shell is filled with molec-ular or colloidal LCs which show planar anchoring at thetwo curved interfaces (Liang et al. , 2013, 2011; Lopez-Leon et al. , 2012). Also air bubbles covered by microrodshave been prepared and studied in real-space (Zhou et al. ,2008). Moreover, topological defects for charged colloidalspheres confined on a sphere were experimentally studied(Guerra et al. , 2018). Ellipsoidal colloids bound to curvedfluid-fluid interfaces with negative Gaussian curvature(Liu et al. , 2018) and spherical droplets covered with as-pherical surfactants (Yang et al. , 2018) were explored.Even living and motile “particles” like cells (Bade et al. ,2017) and suspensions of microtubules and kinesin (Ellis et al. , 2018; Keber et al. , 2014) were recently studied onsurfaces with non-constant curvature. In all these studies a tight coupling between topology, geometric propertiesand defect dynamics is observed. In equilibrium defectsare positioned according to geometric properties of thesurface (Kralj et al. , 2011; Lubensky and Prost, 1992;Nelson, 2002). Creation and annihilation of defects canresult from geometric interaction, leading to different re-alizations of the Poincar´e-Hopf theorem on topologicallyequivalent but geometrically different surfaces (Nestler et al. , 2018). Also changes in the phase diagram can beinduced by the geometry, e. g. allowing for coexistenceof isotropic and nematic phases in surface LCs (Nitschke et al. , 2018). In active systems the observed phenomenaare even richer, including, e. g., oscillating defect patterns(Alaimo et al. , 2017; Keber et al. , 2014) and circulatingband structures (Sknepnek and Henkes, 2015). The ef-fect of hydrodynamics on these phenomena is more orless unexplored.Most of the theoretical studies of these phenomena useparticle methods. Despite the interest in such methods acontinuous description would be more essential for pre-dicting and understanding the macroscopic relation be-tween type and position of the defects and geometricproperties of the surface. Also the influence of hydro-dynamics and dynamic shape changes on these relationswould be much more appropriate to study within a con-tinuous approach. However, a coherent model, which ac-counts for the complex interplay between topology, ge-ometry, defect interactions, hydrodynamics and shapechanges, is still lacking. In (Napoli and Vergori, 2016), anattempt in this direction is proposed but for a fixed sur-face. We here extend this approach and propose a min- a r X i v : . [ phy s i c s . c o m p - ph ] J a n imal continuous surface hydrodynamic LC model, whichcontains the evolution of the surface, tangential polar or-dering and surface hydrodynamics. The passive modelis derived as a thin-film limit of the simplified Ericksen-Leslie model (Lin and Liu, 2000). We describe a numeri-cal approach to solve this model on general surfaces anddemonstrate by simulations various expected and someunexpected phenomena on ellipsoidal and toroidal sur-faces. These phenomena result from the tight couplingof the geometry with the fluid velocity and the directorfield. However, a full exploration of the rich nonlinearphenomena resulting from these relations goes beyondthe scope of the paper. This also holds for the extensionto active systems. The proposed model of surface activepolar viscous gels follows as a thin-film limit of a three-dimensional active polar viscous gel model, which com-bines a more general Ericksen-Leslie model with activecomponents (Kruse et al. , 2004; Simha and Ramaswamy,2002). The model can be derived and numerically solvedusing the same concepts. We here only formulate themodel and exemplarily demonstrate numerically the ef-fect of the underlying curvature on the location of topo-logical defects in an active system. Throughout the wholepaper we consider the evolution of the surface to be pre-scribed and the surface to be decoupled from the sur-rounding bulk phases in order to highlight the surfacehydrodynamics and its coupling with topological and ge-ometric effects. II. THE ERICKSEN-LESLIE MODEL
The Ericksen-Leslie model (Ericksen, 1961, 1976;Leslie, 1968) is an established model for LCs, whose re-laxation dynamics are affected by hydrodynamics. In(Lin and Liu, 2000) a simplified model was introducedand analyzed. This system already retains the mainproperties of the original Ericksen-Leslie model (Hu andWu, 2013; Huang et al. , 2014; Lin et al. , 2010; Wang et al. , 2013) and will be considered as a starting modelto derive a surface hydrodynamic LC model by meansof a thin-film limit, see appendix B. The resulting sim-plified surface Ericksen-Leslie model (cf. eqs. (B30) –(B32)) reads π S ∂ t v + ∇ S v v = v n ( B v + ∇ S v n ) − ∇ S p S − ν ∆ dR v + ν (2 κ v + ∇ S ( v n H ) − S ( v n B )) − λ div S σ E S (1)div S v = v n H (2) π S ∂ t p + ∇ S v p = η (cid:0) ∆ DG p − B p (cid:1) − ω n (cid:16) (cid:107) p (cid:107) S − (cid:17) p (3)where v ( t ) ∈ T S ( t ) denotes the tangential surface veloc-ity, p ( t ) ∈ T S ( t ) the tangential director field, represent-ing the averaged molecular orientation, p S ( x , t ) ∈ R the surface pressure and σ E S = ( ∇ S p ) T ∇ S p + ( B p ) ⊗ ( B p )the extrinsic surface Ericksen stress tensor. The model isdefined on a compact smooth Riemannian surface S ( t ).We consider initial conditions v ( x , t = 0) = v ( x ) ∈ T x S (0) and p ( x , t = 0) = p ( x ) ∈ T x S (0). The posi-tive constants ν , λ and η denote the fluid viscosity, thecompetition between kinetic and elastic potential energyand the elastic relaxation time for the molecular orien-tation field, respectively. κ is the Gaussian curvature, H the mean curvature, B the shape operator, v n a pre-scribed normal velocity of the surface and ω n a penaliza-tion parameter to enforce (cid:107) p (cid:107) = 1 weakly. T x S ( t ) is thetangent space on x ∈ S ( t ), T S ( t ) = (cid:116) x ∈S ( t ) T x S ( t ) thetangent bundle, π S the projection to the tangential spacew. r. t. the surface S ( t ) and ∇ S v , ∇ S , div S , ∆ dR as wellas ∆ DG the covariant directional derivative, covariantgradient, surface divergence, Laplace-deRham operatorand Bochner Laplacian, respectively. The system com-bines an incompressible surface Navier-Stokes equation(Jankuhn et al. , 2017; Miura, 2017; Reuther and Voigt,2015; Yavari et al. , 2016) with a weak surface Frank-Oseen model (Nestler et al. , 2018) on an evolving surface.For a general discussion on transport of vector-valuedquantities on evolving surfaces we refer to (Nitschke andVoigt, 2019). The used formulation with the projectionoperator π S requires the presence of an embedding space,which is R in our case, see appendix B for details.For λ = 0 eqs. (1) and (2) are the surface Navier-Stokes equation for an incompressible surface fluid on anevolving surface. These equations can be obtained asa thin-film limit of the three-dimensional Navier-Stokesequation in an evolving domain (Miura, 2017) or by avariational derivation (Yavari et al. , 2016). If only a sta-tionary surface is considered, i. e. v n = 0, the equationsreduce to the incompressible surface Navier-Stokes equa-tion as considered in (Arroyo and DeSimone, 2009; Ebinand Marsden, 1970; Mitrea and Taylor, 2001; Nitschke et al. , 2017, 2012; Reuther and Voigt, 2015, 2018). Com-pared with its counterpart in flat space, not only theoperators are replaced by the corresponding surface op-erators, also an additional contribution from the Gaus-sian curvature arises. This additional term results fromthe surface divergence of the surface strain rate tensor,see (Arroyo and DeSimone, 2009; Jankuhn et al. , 2017).The unusual sign results from the definition of the surfaceLaplace-deRham operator (Abraham et al. , 1988). Eq.(3) with v = , v n = 0 and the Laplace-deRham operator ∆ dR instead of the Bochner Laplacian ∆ DG has been de-rived as a thin-film limit in (Nestler et al. , 2018) and mod-els the L -gradient flow of a weak surface Frank-Oseenenergy. The different operators result from different one-constant approximations in the Frank-Oseen energy, seeappendix B for details. Again an additional geometricterm enters in this equation if compared with the corre-sponding model in flat space. The term with the shapeoperator B results from the influence of the embedding(Napoli and Vergori, 2012; Nestler et al. , 2018; Segatti et al. , 2014). The coupled system eqs. (1) - (3) with v n = 0 can be considered as the surface counterpart ofthe model in (Lin and Liu, 2000). Related surface modelshave been proposed and analyzed in (Napoli and Vergori,2016; Shkoller, 2002). The model in (Shkoller, 2002) isderived from a variational principle on a stationary sur-face and thus only contains intrinsic terms. It differs fromeqs. (1) - (3) with v n = 0 by the extrinsic term B p andthe extrinsic contribution in the surface Ericksen stresstensor ( B p ) ⊗ ( B p ). The model in (Napoli and Vergori,2016) coincides with our formulation with v n = 0 if aspecific parameter set is considered, see also (Napoli andVergori, 2010). However, note that in their notation thesymbol ˜ ∇ S denotes the surface gradient operator, whilewe use ∇ S as the covariant gradient operator. Both arerelated to each other by ∇ S p + ν ⊗ B p = ˜ ∇ S p , where ν denotes the surface normal. III. NUMERICAL METHOD
Eqs. (1) - (3) are a system of vector-valued surfacePDEs. Numerical approaches have been developed forsuch equations on general surfaces only recently, see (Ol-shanskii et al. , 2018; Reuther and Voigt, 2018) for thesurface (Navier-)Stokes equation, (Nestler et al. , 2018)for the surface Frank-Oseen model and (Hansbo et al. ,2016) for a surface vector-Laplace equation. Earlier ap-proaches using vector spherical harmonics, e. g. (Free-den and Schreiner, 2008; Nestler et al. , 2018; Praeto-rius et al. , 2018), are restricted to a sphere or radialmanifold shapes (Gross and Atzberger, 2018) and ap-proaches which rewrite the surface Navier-Stokes equa-tion in a surface vorticity-stream function formulation(Mickelin et al. , 2018; Nitschke et al. , 2012; Reusken,2018; Reuther and Voigt, 2015) are limited to surfaceswith genus g ( S ) = 0, see (Nitschke et al. , 2017; Reutherand Voigt, 2018) for details. For the numerical solution ofeqs. (1) - (3) we combine the methods in (Nestler et al. ,2018; Reuther and Voigt, 2018) in an operator splittingapproach. The idea behind these methods is to extendthe variational space from vectors in T S to vectors in R ,while penalizing the normal components. This allows tosplit the vector-valued surface PDE into a set of cou-pled scalar-valued surface PDEs for each component forwhich established numerical methods are available, seethe review (Dziuk and Elliott, 2013).The corresponding extended problem to eqs. (1) - (3) reads π S ∂ t (cid:98) v + ∇ S (cid:98) v (cid:98) v = v n ( B (cid:98) v + ∇ S v n ) − ∇ S p S − ν (cid:98) ∆ dR (cid:98) v + ν (2 κ (cid:98) v + ∇ S ( v n H ) − S ( v n B )) − λ div S (cid:98) σ E S − α v ( (cid:98) v · ν ) ν (4)div S (cid:98) v = v n H (5) π S ∂ t (cid:98) p + ∇ S (cid:98) v (cid:98) p = η (cid:16) (cid:98) ∆ DG (cid:98) p − B (cid:98) p (cid:17) − ω n (cid:16) (cid:107) (cid:98) p (cid:107) − (cid:17) (cid:98) p − α p ( ν · (cid:98) p ) ν (6)with (cid:98) v = (cid:98) v x e x + (cid:98) v y e y + (cid:98) v z e z , (cid:98) p = (cid:98) p x e x + (cid:98) p y e y + (cid:98) p z e z ∈ R and (cid:98) σ E S = ( ∇ S (cid:98) p ) T ∇ S (cid:98) p + ( B (cid:98) p ) ⊗ ( B (cid:98) p ). We furtheruse div S (cid:98) v = ∇ · (cid:98) v − ν · ( ∇ (cid:98) v · ν ), rot S (cid:98) v = − div S ( ν × (cid:98) v )and (cid:98) ∆ DG (cid:98) p = div S ∇ S (cid:98) p and (cid:98) ∆ dR (cid:98) v = − (rot S rot S (cid:98) v −∇ S ( v n H )) since div S (cid:98) v = − v n H . The normal compo-nents (cid:98) v · ν and (cid:98) p · ν are penalized by the additional terms α v ( ν · (cid:98) v ) ν and α p ( ν · (cid:98) p ) ν with penalization parameters α v and α p . For convergence results in α v and α p for the sur-face Navier-Stokes and the surface Frank-Oseen problemwe refer to (Nestler et al. , 2018; Reuther and Voigt, 2018).Without these penalization terms the system of equations(4) - (6) is an under-determined problem, since the vec-tor fields are considered in R and therefore the normalcomponents are completely arbitrary, see (Nestler et al. ,2018; Reuther and Voigt, 2018) for details. Eqs. (4) - (6)can now be solved for each component (cid:98) v x , (cid:98) v y , (cid:98) v z , (cid:98) p x , (cid:98) p y , (cid:98) p z and p S using standard approaches for scalar-valuedproblems on surfaces, such as the surface finite elementmethod (Dziuk and Elliott, 2007a,b, 2013), level set ap-proaches (Bertalmio et al. , 2001; Dziuk and Elliott, 2008;Greer et al. , 2006; St¨ocker and Voigt, 2008) or diffuse in-terface approximations (R¨atz and Voigt, 2006). We con-sider a simple operator splitting approach and solve eqs.(4) - (5) and eq. (6) iteratively in each time step, employ-ing the same surface finite element discretizations as in(Nestler et al. , 2018; Reuther and Voigt, 2018). A semi-implicit Euler discretization in time is used. Thereby,the nonlinear transport term in eq. (4) and the norm-1 penalization term in eq. (6) are linearized in time bya Taylor-1 expansion and the transport term in eq. (6)as well as the term including the surface Ericksen stresstensor in eq. (4) are coupling terms in the operatorsplitting scheme. Additionally, we employ an adaptivetime-stepping scheme which is based on the combinationof changes in the surface Frank-Oseen energy and theCourant-Friedrichs-Lewy (CFL) condition. For more de-tails we refer to appendix A. The resulting discrete equa-tions are implemented in the FEM-toolbox AMDiS (Veyand Voigt, 2007; Witkowski et al. , 2015), where we addi-tionally use a domain decomposition ansatz to efficientlydistribute the workload on many cores systems. FIG. 1 Top: Evolution of the director field (cid:98) p on a stationaryellipsoid of the dry case (top row) and the wet case (bottomrow) for t = 1, 2, 3, 4, 6 (left to right). Bottom: SurfaceFrank-Oseen energy F P and surface kinetic energy F kin vs.time t . IV. RESULTS
In the following simulations we use λ = 0 . α v = 10 , ω n = 10 and α p = 10 , where all parameters are treatedas nondimensional. We compare the solution of eqs. (4)- (6) (the so-called wet case) and the solution of eq. (6)with (cid:98) v = 0 (the so-called dry case). To highlight thedifferences we take the surface Frank-Oseen energy F P and the surface kinetic energy F kin into account, whichread in the extended form incorporating the penalizationterm, F P := (cid:90) S η (cid:16) (cid:107)∇ S (cid:98) p (cid:107) + ( B (cid:98) p ) (cid:17) + ω n (cid:0) (cid:107) (cid:98) p (cid:107) − (cid:1) d S + α p (cid:90) S ( (cid:98) p · ν ) d SF kin := 12 (cid:90) S (cid:98) v d S .First, we consider eqs. (4) - (6) on a stationary,i. e. v n = 0, ellipsoidal shape with major axes param-eters (0 . , . , . (cid:98) p = ∇ S ψ / (cid:107)∇ S ψ (cid:107) with ψ = x /
10 + x + x /
10 and x = ( x , x , x ) T the Euclidean coordinate vector. Thelatter generates a vector field with two +1 defects – to bemore precise a source and a sink defect – and an out-of-equilibrium solution. Furthermore, we use ν = 2, η = 0 . h m = 1 . · − , where h m denotes the maximummesh size. Figure 1 shows the influence of the hydrody-namics on the dynamical evolution of the director field.The two defects, which fulfill the Poincar´e-Hopf theorem,evolve towards the geometrically favorable positions ofhigh Gaussian curvature, the director field aligns withthe minimal curvature lines of the geometry and as inflat space the hydrodynamics enhances the evolution to-wards the equilibrium configuration, which coincides forthe dry and the wet case.In the next example we consider a stationary toruswith major radius R = 2, minor radius r = 0 . x -axis as symmetry axis. Again we use the trivial solu-tion as initial condition for the velocity (cid:98) v and a random(normalized) vector field for the director field (cid:98) p . Here,we use the simulation parameters ν = 1 and η = 0 .
4. Themaximum mesh size is fixed at h m = 2 . · − . All otherparameters are equal to that used in Figure 1. In Figure 2we focus on the annihilation of defects in one realization.Figure 2 (top) shows the evolution of the director field (cid:98) p for the dry and the wet case. Again in the wet case thedynamics is enhanced, which is quantified by the strongeroverall decay of the surface Frank-Oseen energy, cf. Fig-ure 2 (middle). Additionally, in Figure 2 (bottom) thecorresponding flow field (cid:98) v is shown for the considered an-nihilation of a source (+1) and a saddle ( −
1) defect in thedirector field (cid:98) p . After all defects are annihilated, whichagain is in accordance with the Poincar´e-Hopf theorem,the velocity field (cid:98) v decays to zero and the director field (cid:98) p aligns with the minimal curvature lines of the geometry.The reached equilibrium configuration coincides for boththe dry and the wet case.While in the two previous examples the expected min-imal energy configuration was reached, we now con-sider an initial condition for which only a local mini-mum can be reached. We use (cid:98) p = ∇ S ψ / (cid:107)∇ S ψ (cid:107) with ψ = exp (cid:0) − ( x − m ) / (cid:1) and m = ( R, , r ) T as initialcondition for the director field. This produces two ± π/ ν = 1, η = 0 . h m = 2 . · − .In a flat geometry with zero curvature these two pairswould annihilate. However, due to the geometric inter-action in the present case resulting from the difference ofthe Gaussian curvature inside and outside of the torus,the reached nontrivial defect configuration is stable andthe two ± − FIG. 2 Top: Evolution of the director field (cid:98) p on a torus of thedry case (top row) and the wet case (bottom row) for t = 0 . .
5, 6 . F P and surface kinetic energy F kin vs. time t . Bottom:Velocity field (cid:98) v for the annihilation of a source (+1, left) anda saddle ( −
1, right) defect in the director field (cid:98) p (red dots)for t = 3, 3 .
75, 3 .
81, 4 .
05, 4 .
5, 5 . are attracted to regions with negative Gaussian curva-ture, i. e. the inner of the torus, and +1 defects are at-tracted to regions with positive Gaussian curvature, i. e.the outer of the torus, see Figure 3. The reached con-figuration, is a local minimum with a significantly largersurface Frank-Oseen energy F P as the defect-free config-uration. In this example we did not find any significantdifference between the dry and the wet case, when thezero initial condition for the velocity (cid:98) v is used. However,if we use a Killing vector field for the velocity as initialcondition, i. e. (cid:98) v = 1 / − x , x , T , cf. (Nitschke et al. , FIG. 3 Top: Schematic defect positions of the initial condi-tion (left) and the final configuration (right) on a torus withthe analytical initial condition for the director field (cid:98) p and zeroinitial condition for the velocity field (cid:98) v . Red dots are indicat-ing +1 defects (sources or sinks) and blue dots are indicating − (cid:98) p for t = 1, 5, 25 (left to right).FIG. 4 Surface Frank-Oseen energy F P and surface kineticenergy F kin vs. time t for the analytical initial conditionfor the director field (cid:98) p and the killing vector field as initialcondition for the velocity (cid:98) v . (cid:98) v = 0, withthe rotation angle depending on the strength of the initialvelocity and the viscosity.We now let the ellipsoid from Figure 1 evolve by pre-scribing the normal velocity v n , such that the ellipsoidchanges to a sphere and afterwards to an ellipsoid with adifferent axis orientation and vice versa to obtain a shapeoscillation. The surface area remains constant during theevolution. Figure 5 shows schematically the evolution ofthe geometry and the axes parameters for one period ofoscillation. We use the same simulation parameters andinitial conditions as considered in Figure 1. The evolu- FIG. 5 Left: Schematic description of the ellipsoid evolutionfor a half period of oscillation. Descending gray scale indi-cates increasing time. The motion in the second half of theoscillation is reversed, respectively. Right: Major axes pa-rameters for the ellipsoid over a full period of oscillation. Thetime of one oscillation period is considered to be T = 160.The major axes parameters are chosen such that the surfacearea of the ellipsoid is conserved over time. tion of the director field (cid:98) p is shown in Figure 6, againfor the dry (top) and the wet (bottom) case. The defectpositions again reallocate at their geometrically favor-able position. However, due to the change in the geome-try the time scale for the reallocation competes with thetime-scale for the shape changes. The enhanced evolutiontowards the minimal energy configuration with hydrody-namics becomes even more significant in these situations.Already slight modifications of the geometry are enoughto push the defect after crossing the sphere configurations(with no preferred defect position) to the energeticallyfavorable state. In the dry case there is a strong delayand much stronger shape changes are needed to push thedefect to the energetically favorable position. First anenergy barrier for reallocating the defect position has tobe overcome, which is shown by the further increase ofthe red line after the blue line has already dropped aftercrossing the sphere configuration in Figure 6 (middle).The parameters and the initial condition are further cho-sen in such a way that the defects in the dry case notquite reach the position at the poles if the shape evolu-tion crosses the sphere. In the wet case they have movedbeyond. This results in a constant orientation in the drycase and a flipping of the orientation of the director fieldin the wet case in each oscillation. The final configu-ration in Figure 6 after completing one oscillation cycleis energetically equivalent for the dry and the wet caseeven if the orientation of the director field (cid:98) p differs, seealso the video in the supplementary material. This be-havior clearly depends on the used parameters. However,it also demonstrates the strong influence hydrodynamicsmight have in such highly nonlinear systems, where thetopology, geometric properties and defect dynamics arestrongly coupled.These examples together with the demonstrated en-ergy reduction by creation of additional defects in ge- FIG. 6 Top: Evolution of the director field (cid:98) p on an evolvingellipsoid of the dry case (top row) and the wet case (bottomrow) for t = 3, 40, 95, 145, 160 (left to right). Middle: SurfaceFrank-Oseen energy F P and surface kinetic energy F kin vs.time t for the first period of oscillation. Bottom: SurfaceFrank-Oseen energy F P and surface kinetic energy F kin vs.time t over five periods of oscillation. ometrically favored positions in (Nestler et al. , 2018),which is expected to hold also for the wet case, leads toa very rich phase space, considering geometric and ma-terial properties, whose exploration is beyond the scopeof this paper. FIG. 7 Surface Frank-Oseen energy F P and surface kineticenergy F kin vs. time t for the simulation with the dampedoscillation of the defects around the minimal defect configu-ration. V. DISCUSSION
Eqs. (1) - (3) have been derived as a thin-film limitof a three-dimensional simplified Ericksen-Leslie model,see appendix B. In (Shkoller, 2002) a similar modelwas proposed, which differs from eqs. (1) – (3) with v n = 0 in the extrinsic contributions. Especially thesurface Ericksen stress tensor is considered to be σ E S =( ∇ S p ) T ∇ S p . To show the strong difference between thisintrinsic and the extrinsic surface Ericksen stress tensor σ E S = ( ∇ S p ) T ∇ S p + ( B p ) ⊗ ( B p ) considered here andin (Napoli and Vergori, 2016), we come back to the sta-tionary ellipsoid in Figure 1. We use slightly differentparameters, i. e. ν = 0 . η = 0 . λ = 1, which leadto a damped oscillation of the defects around the ener-getically favorable positions before they reach the finalstate configuration as in Figure 1. In Figure 7 the dif-ferences in the time evolution of the surface Frank-Oseenenergy as well as the surface kinetic energy are shownfor both cases the intrinsic and extrinsic surface Erick-sen stress tensor. The influence of the hydrodynamics ismuch stronger for the extrinsic surface Ericksen stress.Together with the example in Figure 6 such differencesin the dynamics might have a huge impact on the overallevolution if also shape changes are considered.All results so far are for the simplified surface Ericksen-Leslie model. However, appendix B provides all necessarytools to do the thin-film analysis also for more compli-cated systems, such as more general Ericksen-Leslie mod-els or active versions of them. Here, we provide the for-mulation for a surface active polar viscous gel, see (Kruse et al. , 2004; Simha and Ramaswamy, 2002) and (Marth et al. , 2015; Tjhung et al. , 2012) for the considered three-dimensional formulation, which correspond to eqs. (B1) - (B3) with boundary conditions (B4) - (B6), i. e. ∂ t V + ∇ U V = −∇ P Ω h + ν ∆ V + div σ A − λ div σ E div V = 0 ∂ t P + ∇ U P = H + α DP + Ω P with σ A = 12 ( P ⊗ H − H ⊗ P ) − α P ⊗ H + H ⊗ P )+ β P ⊗ PH = η ∆ P − ω n (cid:16) (cid:107) P (cid:107) h − (cid:17) PD = 12 (cid:16) ∇ V + ( ∇ V ) T (cid:17) Ω = 12 (cid:16) ∇ V − ( ∇ V ) T (cid:17) and α, β ∈ R . The Navier-Stokes equation now containsadditional distortion and active stresses, combined in σ A ,while in the director field equation additional contribu-tions from the strain rate tensor D and the vorticity ten-sor Ω arise. The corresponding thin-film limit reads π S ∂ t v + ∇ S v v = v n ( B v + ∇ S v n ) + ν (cid:0) − ∆ dR v + 2 κ v (cid:1) + ν ( ∇ S ( v n H ) − S ( v n B )) − ∇ S p S + div S σ A S − λ div S σ E S − − α p T B ( H v + ∇ S v n )) p + 1 + α p T B v ) B p (7)div S v = v n H (8) π S ∂ t p + ∇ S v p = h + α D S p + Ω S p − αv n B p (9)with σ A S = 12 ( p ⊗ h − h ⊗ p ) − α p ⊗ h + h ⊗ p )+ β p ⊗ ph = η (cid:0) ∆ DG p − B p (cid:1) − ω n (cid:16) (cid:107) p (cid:107) S − (cid:17) p D S = 12 (cid:16) ∇ S v + ( ∇ S v ) T (cid:17) Ω S = 12 (cid:16) ∇ S v − ( ∇ S v ) T (cid:17) . Besides the corresponding surface operators and the ad-ditional geometric coupling terms with the shape opera-tor B and the mean curvature H we also obtain an ex-plicit appearance of the normal velocity v n in the directorfield equation. Overall the additional terms in the moregeneral Ericksen-Leslie model lead to an even tighter cou-pling between geometric properties and dynamics. Thedescribed numerical approach, see appendix A, can beadapted to also solve the surface active polar viscous gelmodel proposed in eqs. (7) - (9). Figure 8 shows resultson a torus. We use the same torus and the same initial FIG. 8 Top: Snapshot of the director field (cid:98) p of the surfaceactive polar viscous gel model on the torus for t = 14 .
3. Bot-tom: Number of defects per area for the inner ( κ >
0) andouter ( κ <
0) region of the torus vs. time t . conditions for the velocity field (cid:98) v and the director field (cid:98) p as considered in Figure 2. The parameters are adapted to α = 1 . β = 20, λ = 0 . ν = 200, η = 0 .
01 and ω n = 5,while all other parameters remain unchanged. Figure 8(top) shows a snapshot of the director field, see also thevideo in the supplementary material. The significantlyreduced parameter η yields smeared out defects and pro-motes the annihilation and creation of new defects. InFigure 8 (bottom) the number of defects per area againstthe time is plotted. Thereby, we distinguish between theinner ( κ <
0) and the outer ( κ >
0) region of the torusand observe slightly more defects per area in the innerpart. This might be due to the stronger geometric force(resulting from a higher absolute value of the Gaussiancurvature in the inner region), the continuous creationand annihilation of defects as well as the fact that de-fects of opposite topological charge are attracted to eachother. As in the passive case, a detailed analysis of suchphenomena has to be discussed elsewhere.
ACKNOWLEDGMENTS
This work was financially supported by the GermanResearch Foundation (DFG) through project Vo899-19.We used computing resources provided by J¨ulich Super-computing Centre within project HDR06.
Appendix A: Numerics
To efficiently solve the surface Navier-Stokes equationin (Reuther and Voigt, 2018) a heavy assembly workloadwas avoided by applying ν × to eq. (4) and consider-ing the rotated velocity field (cid:98) w := ν × (cid:98) v . Since onlythe Rot S rot S ( · ) operator occurs as second order oper-ator we here can also use the same idea to reduce theassembly costs. Thus, the rotated version of eqs. (4) and(5) with tangential penalization of the rotated velocity (cid:98) w now reads π S ∂ t (cid:98) w + ∇ (cid:98) w (cid:98) w = − rot S p + ν ( ∇ S div S (cid:98) w + 2 κ (cid:98) w )+ f g + f (cid:98) w − λ ν × div S (cid:98) σ E S (A1)rot S (cid:98) w = v n H (A2)where we used for convenience the abbreviations f g := v n rot S v n + 2 ν ( H rot S v n − ν × ( B∇ S v n )) f (cid:98) w := − v n ν × ( B ( ν × (cid:98) w )) − α v ( (cid:98) w · ν ) ν and the alternative form of the viscous terms proposedin (B23). Time discretization
For the discretization in time we again use the sameapproach proposed in (Reuther and Voigt, 2018). Letthe time interval [0 , t end ] be divided into a sequence ofdiscrete times 0 < t < t < ... with time step width τ m = t m − t m − . Thereby, the superscript denotes thetimestep number. The vector field (cid:98) w m ( x ) correspond tothe respective rotated velocity field (cid:98) w ( x , t m ). All otherquantities follow the same notation. The time deriva-tive is approximated by a standard difference quotientand a Chorin projection method (Chorin, 1968) is ap-plied to eqs. (A1) and (A2). Furthermore, we definethe discrete time derivatives d (cid:98) w τ m := τ m (cid:0) (cid:98) w ∗ − π S (cid:98) w m − (cid:1) and d (cid:98) p τ m := τ m (cid:0)(cid:98) p m − π S (cid:98) p m − (cid:1) , with π S the projectionto the surface at time t m . Thus, we get a time-discreteversion of eqs. (A1), (A2) and (6)d (cid:98) w τ m + ∇ (cid:98) w m − (cid:98) w ∗ = ν ( ∇ S div S (cid:98) w ∗ + 2 κ (cid:98) w ∗ ) + f g + f (cid:98) w ∗ − λ ν × div S (cid:98) σ E S (A3) τ m ∆ S p m = rot S (cid:98) w ∗ − v n H (A4) (cid:98) w m = (cid:98) w ∗ − τ m rot S p m (A5)d (cid:98) p τ m + ∇ (cid:98) v m (cid:98) p m = η (cid:16) (cid:98) ∆ DG (cid:98) p m − B (cid:98) p m (cid:17) − ω n (cid:0) (cid:107) (cid:98) p m − (cid:107) − (cid:1) (cid:98) p m − α p ( ν · (cid:98) p m ) ν (A6)where (cid:98) σ E S is evaluated at the old timestep, i. e. (cid:98) σ E S = (cid:0) ∇ S (cid:98) p m − (cid:1) T ∇ S (cid:98) p m − + (cid:0) B (cid:98) p m − (cid:1) ⊗ (cid:0) B (cid:98) p m − (cid:1) . Note thatfor readability we used a Taylor-0 linearization of thetransport term in (A3) and the norm-1 penalization termin (6). In the simulations from above we performed aTaylor-1 linearization, see (Nitschke et al. , 2017) and(Nestler et al. , 2018) for details. Spatial discretization
The considered extension of the tangential vector fieldsto the Euclidean space allows us to apply the surfacefinite element method (Dziuk and Elliott, 2013) for eachcomponent of the respective vector field. Let S h denotethe interpolation of the surface S ( t m ) at time t = t m suchthat S h := (cid:83) T ∈T T with a conforming triangulation T .Furthermore, we introduce the finite element space V h ( S h ) = (cid:8) v h ∈ C ( S h ) : v h | T ∈ P , ∀ T ∈ T (cid:9) which is used twice as trail and test space and thestandard L scalar product on V h ( S h ), ( α, β ) := (cid:82) S h (cid:104) α, β (cid:105) d S . By using an operator splitting techniquewe decouple the hydrodynamic and the director fieldequation in the following way. First the surface finite el-ement approximation of eqs. (A3), (A4) is solved, whichreads: find (cid:98) w ∗ i , p m ∈ V h ( S h ) s.t. ∀ ξ, η ∈ V h ( S h ) (cid:16) (d (cid:98) w τ m ) i + ∇ (cid:98) w m − (cid:98) w ∗ i − νκ (cid:98) w ∗ i − ( f g + f (cid:98) w ∗ ) i , ξ (cid:17) = − (cid:16) ν div S (cid:98) w ∗ , ( ∇ S ξ ) i (cid:17) − (cid:16) λ (cid:98) σ E S , ∇ S ( ν × ( ξ e i )) (cid:17)(cid:16) τ m ∇ S p m , ∇ S η (cid:17) = (cid:16) v n H − rot S (cid:98) w ∗ , η (cid:17) for i = x, y, z . The resulting vector field (cid:98) w ∗ is then usedto determine (cid:98) w m by the pressure correction step in eq.(A5). The transformation (cid:98) v m = − ν × (cid:98) w m leads to thevelocity field at the new timestep t m . Finally, the surfacefinite element approximation of eq. (A6) is solved, whichreads: find (cid:98) p mi ∈ V h ( S h ) s.t. ∀ ξ ∈ V h ( S h ) (cid:16) (d (cid:98) p τ m ) i + ∇ (cid:98) v m (cid:98) p mi , ξ (cid:17) = (cid:16) η ∇ S (cid:98) p m , ∇ S ( ξ e i ) (cid:17) − (cid:16) η B (cid:98) p mi , ξ (cid:17) − (cid:16) ω n (cid:0) (cid:107) (cid:98) p m − (cid:107) − (cid:1) (cid:98) p mi + α p ( ν · (cid:98) p m ) ν i , ξ (cid:17) for i = x, y, z . Pressure relaxation schemes
In some situations it is useful to modify the Chorin pro-jection scheme (A3), (A4) and (A5). To be more precisethe resulting finite element matrix of the pressure equa-tion (A4) is sometimes ill-conditioned, especially whenthe term rot S (cid:98) w ∗ is big compared to the others. The so-lution of eq. (A4) can be seen as the steady-state solutionof a heat conduction equation where the heat source is determined by the right hand side of eq. (A4). There-fore, we add a relaxation scheme in form of a discretetime derivative on a different timescale to the left handside of eq. (A4), i. e.1 τ ∗ (cid:0) p m +1 ,l +1 − p m +1 ,l (cid:1) − τ m ∆ S p m +1 ,l +1 = − rot S (cid:98) w ∗ + v n H , (A7)where τ ∗ denotes the timestep and l the timestep numberon the different timescale. Instead of solving eq. (A4) theiterative process in eq. (A7) is performed until a steady-state is reached, which is then used in the correction stepin eq. (A5). Appendix B: Thin film limit
We assume a regular moving surface S ( t ) ⊂ R withoutboundaries and a thin film Ω h ( t ) := S ( t ) × [ − h/ , h/ ⊂ R of sufficiently small thickness h , such that thethin film parametrization X ( t, y , y , ξ ) = x ( t, y , y ) + ξ ν ( t, y , y ) is injective for the surface parametrization x ( t, · , · ). Thereby, y and y denote the local surfacecoordinates, ν ( t, · , · ) the surface normal field and ξ ∈ [ − h/ , h/
2] is the local normal coordinate. Since the thinfilm is moving according to the surface, the parametriza-tion X is not unique, which arises from the choice ofan observer within the thin film. For a pure Eulerianobserver, i. e. for the observer velocity W = ∂ t X = 0,we are not able to formulate proper intrinsic physicsat the surface S for h →
0. To overcome this is-sue, we choose a transversal observer as the surface ob-server parametrization x , i. e. Eulerian in the tangen-tial space and Lagrangian in normal direction and hence ∂ t x = ∂ t X (cid:12)(cid:12) S = v n ν , where v n = W ξ (cid:12)(cid:12) S is the normalsurface velocity of S . To ensure a constant thickness h of the thin film and that the surface S is located in themiddle of the thin film over time, we stipulating the sametransversal behavior for both boundary surfaces. Thus,we get ∂ t X (cid:12)(cid:12) ∂ Ω h = v n ν ± h ∂ t ν = W ξ (cid:12)(cid:12) ∂ Ω h ν and there-fore W ξ (cid:12)(cid:12) ∂ Ω h = v n , since ∂ t ν is always tangential on theboundaries ∂ Ω h .For notational compactness of tensor algebra we usethe thin film calculus presented in (Nitschke et al. , 2018)and (Nestler et al. , 2018) (Appendix) which is basedon Ricci calculus, where lowercase indices i, j, k, . . . de-note components w. r. t. y and y in the surface co-ordinate system and uppercase indices I, J, K, . . . de-note components w. r. t. y , y and ξ in the extendedthree dimensional thin film coordinate system. Met-ric quantities at the surface S are the metric tensor g ij = (cid:104) ∂ i x , ∂ j x (cid:105) R (first fundamental form), the shape op-erator B ij = − (cid:104) ∂ i x , ∂ j ν (cid:105) R (second fundamental form),its square (cid:2) B (cid:3) ij = B ik B kj = (cid:104) ∂ i ν , ∂ j ν (cid:105) R (third fun-damental form), the mean curvature H = tr B = B ii ,0the Gaussian curvature κ = det B (cid:93) = det (cid:8) B ij (cid:9) and theChristoffel symbols Γ kij = g kl ( ∂ i g jl + ∂ j g il − ∂ l g ij ) for co-variant differentiating (e. g. [ ∇ S p ] ik = p i | k = ∂ k p i +Γ ikj p j for a contravariant vector field p ∈ T S = T S ). In thethin film Ω h , the metric tensor G IJ = (cid:104) ∂ I X , ∂ J X (cid:105) R andthe Christoffel symbols L KIJ = G KL ( ∂ I G JL + ∂ J G IL − ∂ L G IJ ) for covariant differentiating (e. g. [ ∇ P ] IK = P I ; K = ∂ K P I + L IKJ P J for a contravariant vector field P ∈ T Ω h = T Ω h ) can be developed at the surface by G ij = g ij − ξB ij + ξ (cid:2) B (cid:3) ij , G ξξ = 1, G ξi = G iξ = 0, G ij = g ij + O ( ξ ), G ξξ = 1, G ξi = G iξ = 0, L kij =Γ kij + O ( ξ ), L ξij = B ij + O ( ξ ), L kiξ = L kξi = − B ki + O ( ξ )and L ξIξ = L ξξI = L Kξξ = 0, see (Nitschke et al. , 2018) and(Nestler et al. , 2018) for details.Our starting point is the simplified local three dimen-sional Ericksen-Leslie model (Lin and Liu, 2000), i. e. ∂ t V + ∇ U V = −∇ P Ω h + ν ∆ V − λ div σ E (B1)div V = 0 (B2) ∂ t P + ∇ U P = η ∆ P − ω n (cid:16) (cid:107) P (cid:107) h − (cid:17) P (B3)in Ω h × R + with fluid velocity V ∈ T Ω h , relative fluidvelocity U = V − W ∈ T Ω h , with respect to the observervelocity W = ∂ t X , director field P ∈ T Ω h , pressure P Ω h ,Ericksen stress tensor σ E = ( ∇ P ) T ∇ P , fluid viscosity ν , competition between kinetic and elastic potential en-ergy λ and elastic relaxation time for the molecular ori-entation field η . Besides initial conditions, we considerhomogeneous Dirichlet boundary conditions for the nor-mal components and Neumann boundary conditions forthe tangential components of the director and homoge-neous Navier boundary conditions for the velocity field,i. e. (cid:104) P , ν (cid:105) Ω h = P ξ = 0 (B4) ∇ ν (cid:16) π (cid:91)∂ Ω h P (cid:17) (cid:91) = { P i ; ξ } = 0 (B5) π (cid:91)∂ Ω h ( ν · L V G ) = { V i ; ξ + V ξ ; i } = 0 (B6)at the boundaries ∂ Ω h in its covariant form. Thereby, L V G denotes the viscous stress tensor in terms of theLie derivative L and π (cid:91)∂ Ω h : T Ω h (cid:12)(cid:12) ∂ Ω h → T ∂ Ω h is theorthogonal covariant projection into the covariant bound-ary tangential space. Note that it holds (cid:104) V , ν (cid:105) Ω h = V ξ = W ξ = v n on ∂ Ω h , which follows from the special choice of thetransversal observer from above.As proposed in (Nitschke et al. , 2018) the boundaryquantities are continuable to the surface S by Taylor ex- pansions at the boundaries, which, e. g., results in P ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) P ξ ; ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) (B7) P i ; ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) P i ; ξ ; ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) (B8) V ξ ; ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) V ξ ; ξ ; ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) (B9) V i ; ξ (cid:12)(cid:12) S + V ξ ; i (cid:12)(cid:12) S = O (cid:0) h (cid:1) (B10) V i ; ξ ; ξ (cid:12)(cid:12) S + V ξ ; i ; ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) . (B11)Note that the right identity of eq. (B9) is achieved byusing V ξ (cid:12)(cid:12) X ( ξ = − h ) = V ξ (cid:12)(cid:12) X ( ξ = h ) = V ξ (cid:12)(cid:12) X ( ξ =0) = v n , the re-lated second order difference quotient and V ξ ; ξ ; ξ = ∂ ξ V ξ ,since vanishing Christoffel symbols L Kξξ .With all these tools from above, we are able to realize athin film limit of eqs. (B1) – (B3) for h →
0, consistently,by Taylor expansion of the equations at the surface. Thethin film Ω h is a flat Riemannian manifold and thereforethe Riemannian curvature tensor vanish, i. e. covariantderivatives commute, and with the continuity equation(B2) we obtain[ ∆ V ] I = [div ∇ V + ∇ div V ] I = V ; KI ; K + V K ; K ; I = V ; KI ; K + V K ; I ; K = [div L V G ] I . (B12)Hence, we have to develop three divergence terms of 2-tensors at the surface in eqs. (B1) – (B3), namely div T (cid:12)(cid:12) S for T being either L V G , σ E or ∇ P . By using eqs. (B7)– (B11) it holds T iξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) , ∂ ξ T iξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) , T iξ ; ξ (cid:12)(cid:12) S = O (cid:0) h (cid:1) . (B13)The covariant tangential components of the divergenceat the surface are[div T ] i (cid:12)(cid:12)(cid:12) S = T Ki ; K (cid:12)(cid:12)(cid:12) S = (cid:16) ∂ K T Ki − L JKi T KJ + L KKL T Li (cid:17) (cid:12)(cid:12)(cid:12) S ,where it holds by using eq. (B13) ∂ K T Ki (cid:12)(cid:12)(cid:12) S = ∂ k T ki (cid:12)(cid:12)(cid:12) S + O (cid:0) h (cid:1) , (B14) L JKi T KJ (cid:12)(cid:12)(cid:12) S = (cid:16) Γ jki T kj + B ik T kξ (cid:17) (cid:12)(cid:12)(cid:12) S + O (cid:0) h (cid:1) , (B15) L KKL T Li (cid:12)(cid:12)(cid:12) S = Γ kkl T li (cid:12)(cid:12)(cid:12) S + O (cid:0) h (cid:1) . (B16)Adding this up and take the metric compatibility of ∇ S into account, we obtain[div T ] i (cid:12)(cid:12)(cid:12) S = (cid:16) T ki (cid:12)(cid:12)(cid:12) S (cid:17) | k − B ik T kξ (cid:12)(cid:12)(cid:12) S + O (cid:0) h (cid:1) = g kl (cid:18)(cid:16) T il (cid:12)(cid:12)(cid:12) S (cid:17) | k − B ik T ξl (cid:12)(cid:12)(cid:12) S (cid:19) + O (cid:0) h (cid:1) . (B17)1Note that all normal derivatives vanished here, i. e. thereis no need for a higher order expansion in ξ of the thin filmChristoffel symbols as a consequence of the used bound-ary conditions. To substantiate the tensor T ∈ T (2) Ω h ,we first observe that V i ; j (cid:12)(cid:12)(cid:12) S = (cid:16) ∂ j V i − L kji V k − L ξji V ξ (cid:17) (cid:12)(cid:12)(cid:12) S = v i | j − v n B ij , (B18) P i ; j (cid:12)(cid:12)(cid:12) S = (cid:16) ∂ j P i − L kji P k − L ξji P ξ (cid:17) (cid:12)(cid:12)(cid:12) S = p i | j + O (cid:0) h (cid:1) , (B19) P ξ ; j (cid:12)(cid:12)(cid:12) S = (cid:16) ∂ j P ξ − L kjξ P k (cid:17) (cid:12)(cid:12)(cid:12) S = B kj p k + O (cid:0) h (cid:1) , (B20)where v i := V i (cid:12)(cid:12) S and p i := P i (cid:12)(cid:12) S , i. e. v := π S V (cid:12)(cid:12) S ∈ T S and in contravariant form p := π S P (cid:12)(cid:12) S ∈ T S . Wefurther obtain[ L V G ] il (cid:12)(cid:12)(cid:12) S = v i | l + v l | i − v n B il = [ L v g − v n B ] il , (cid:2) σ E (cid:3) il (cid:12)(cid:12)(cid:12) S = (cid:104) ( ∇ P ) T ∇ P (cid:105) il (cid:12)(cid:12)(cid:12) S = (cid:0) G jk p j ; i p k ; l + p ξ ; i p ξ ; l (cid:1) (cid:12)(cid:12)(cid:12) S = (cid:104) ( ∇ S p ) T ∇ S p + ( B p ) ⊗ ( B p ) (cid:105) il + O (cid:0) h (cid:1) which we define as the extrinsic surface Ericksen stresstensor σ E S := ( ∇ S p ) T ∇ S p + ( B p ) ⊗ ( B p ), and finally get π S ( ∆ V ) (cid:12)(cid:12)(cid:12) S = div S ( L v g − v n B ) + O (cid:0) h (cid:1) , (B21) π S ( ∆ P ) (cid:12)(cid:12)(cid:12) S = ∆ DG p − B p + O (cid:0) h (cid:1) , π S div σ E (cid:12)(cid:12)(cid:12) S = div S σ E S + O (cid:0) h (cid:1) ,where eq. (B17) was used and ∆ DG := div S ◦∇ S de-notes the Bochner Laplacian. Evaluating eq. (B2) at thesurface and using the boundary condition (B9) and theidentity in eq. (B18) gives0 = div V (cid:12)(cid:12)(cid:12) S = (cid:0) G ij V i ; j + V ξ : ξ (cid:1) (cid:12)(cid:12)(cid:12) S = div S v − v n H + O (cid:0) h (cid:1) . (B22)Furthermore, we introduce the two curl operators rot S : T (1) S → T S for vector fields and Rot S : T S → T (1) S for scalar fields, see (Nestler et al. , 2018) and (Nitschke et al. , 2017) for their definitions. Therefore, rewriting eq.(B21) yields π S ( ∆ V ) (cid:12)(cid:12)(cid:12) S = Rot S rot S v + 2 ( κ v + H∇ S v n ) − B∇ S v n + O (cid:0) h (cid:1) = − ∆ dR v + 2 κ v + ∇ S ( v n H ) − S ( v n B ) + O (cid:0) h (cid:1) , (B23) as a consequence of eq. (B22) and the Weizenb¨ock ma-chinery, i. e. interchanging covariant derivatives w. r. t.the Riemannian curvature tensor of S , cf. (Arroyo andDeSimone, 2009). With (cid:107) P (cid:107) h (cid:12)(cid:12)(cid:12) S = (cid:16) G ij P i P j + ( P ξ ) (cid:17) (cid:12)(cid:12)(cid:12) S = (cid:107) p (cid:107) S + O (cid:0) h (cid:1) we get for the remaining terms on the right hand sides ofeqs. (B3) and (B1) π S (cid:16)(cid:16) (cid:107) P (cid:107) h − (cid:17) P (cid:17) (cid:12)(cid:12)(cid:12) S = (cid:16) (cid:107) p (cid:107) S − (cid:17) p + O (cid:0) h (cid:1) , π S ( ∇ P Ω h ) (cid:12)(cid:12)(cid:12) S = ∇ S p S ,with the surface pressure p S := P Ω h (cid:12)(cid:12) S . Note that it holdsfor the relative velocity at the surface U (cid:12)(cid:12) S = π S V (cid:12)(cid:12) S , i. e. U ξ (cid:12)(cid:12) S = 0 and U i (cid:12)(cid:12) S = v i , on the left hand sides of eqs.(B1) and (B3). Hence, eqs. (B18) and (B19) now read[ ∇ U V ] i (cid:12)(cid:12)(cid:12) S = v k V i ; k (cid:12)(cid:12)(cid:12) S = (cid:2) ∇ S v v − v n B v (cid:3) i [ ∇ U P ] i (cid:12)(cid:12)(cid:12) S = v k P i ; k (cid:12)(cid:12)(cid:12) S = (cid:2) ∇ S v p (cid:3) i + O (cid:0) h (cid:1) .By using ∂ t X = W = W ξ ν , W ξ (cid:12)(cid:12) S = v n and V = V k ∂ k X + V ξ ν (analogously for P ) we obtain for the par-tial time derivatives[ ∂ t V ] i (cid:12)(cid:12)(cid:12) S = (cid:104) ∂ t V , ∂ i X (cid:105) Ω h (cid:12)(cid:12)(cid:12) S = (cid:10) (cid:0) ∂ t V k (cid:1) ∂ k X + V k ∂ k W + ( ∂ t V ξ ) ν + V ξ ∂ t ν , ∂ i X (cid:11) Ω h (cid:12)(cid:12)(cid:12) S = g ik ∂ t v k + v k (cid:104) ( ∂ k W ξ ) ν + v n ∂ k ν , ∂ i X (cid:105) Ω h (cid:12)(cid:12)(cid:12) S − v n (cid:104) ν , ∂ i W (cid:105) Ω h (cid:12)(cid:12)(cid:12) S = g ik ∂ t v k − v n (cid:0) B ik v k + ∂ i v n (cid:1) (B24)and [ ∂ t P ] i (cid:12)(cid:12)(cid:12) S = (cid:104) ∂ t P , ∂ i X (cid:105) Ω h (cid:12)(cid:12)(cid:12) S = (cid:10)(cid:0) ∂ t P k (cid:1) ∂ k X + P k ∂ k W , ∂ i X (cid:11) Ω h (cid:12)(cid:12)(cid:12) S + O (cid:0) h (cid:1) = g ik ∂ t p k − v n B ik p k + O (cid:0) h (cid:1) . (B25)Therefore, we have π S [ ∂ t V + ∇ U V ] (cid:12)(cid:12)(cid:12) S = (cid:0) ∂ t v i (cid:1) ∂ i x + ∇ S v v − v n (2 B v + ∇ S v n ) (B26)for the tangential part of the fluid acceleration in eq.(B1). The same term was proposed in (Yavari et al. ,2016) by variation of the kinetic energy of a moving man-ifold in the context of Lagrangian field theory. Moreover,2we also find this acceleration term in (Nitschke and Voigt,2019), where a covariant material derivative is derived interms of covariant tensor transport through a three di-mensional moving spacetime embedded in a four dimen-sional absolute space. In this context eq. (B26) can beobtained by taking the spatial part of the covariant ma-terial derivative for the special case of velocity fields anda transversal observer. Evaluating the transport term forthe director field P at the surface yields π S [ ∂ t P + ∇ U P ] (cid:12)(cid:12)(cid:12) S = (cid:0) ∂ t p i (cid:1) ∂ i x + ∇ S v p − v n B p ,which can also be found in (Nitschke and Voigt, 2019),but as the spatial part of the covariant material derivativeof a so-called instantaneous vector field from a point ofview of a transversal observer. Finally, under boundaryconditions (B4) – (B6) and h →
0, eq. (B2) and thetangential parts of eqs. (B1) and (B3) reduces to (cid:0) ∂ t v i (cid:1) ∂ i x + ∇ S v v − v n (2 B v + ∇ S v n )= ν (cid:0) − ∆ dR v + 2 κ v + ∇ S ( v n H ) − S ( v n B ) (cid:1) − ∇ S p S − λ div S σ E S (B27)div S v = v n H (B28) (cid:0) ∂ t p i (cid:1) ∂ i x + ∇ S v p − v n B p = η (cid:0) ∆ DG p − B p (cid:1) − ω n (cid:16) (cid:107) p (cid:107) S − (cid:17) p (B29)in S × R + . This system of PDEs has full rank, i. e.it contains five independent coupled equations with fivedegree of freedoms v , v , p , p , depending on an ar-bitrary choice of local coordinates, and p S . A full dis-cussion about the normal parts of eqs. (B1) and (B3)does not belong to this paper. Nevertheless, the nor-mal parts would give us two additional equations, con-sistently w. r. t. h , and two new free scalar valued quan-tities γ , γ ∈ T S . With our boundary conditions andassumption, w. r. t. the moving thin film geometry, thiswould be γ := ∂ ξ P Ω h (cid:12)(cid:12) S and γ := P ξ ; ξ ; ξ (cid:12)(cid:12) S = ∂ ξ P ξ (cid:12)(cid:12) S .Both degrees of freedom would occur as zero order differ-ential terms, i. e. the upcoming normal equations wouldnot have any influence to eqs. (B27) – (B29). Therefore,the thin film limit of the normal part equations can beomitted as long as we are not interested in the quantities γ and γ .The partial time derivatives in eqs. (B27) and (B29)are realized only at the contravariant vector proxies of v and p w. r. t. locally defined charts at the surface. Thereason for the absence of an intrinsic covariant vector op-erator notation for the time derivative (similar to ∇ S ) isthat the time t is not a coordinate of a moving space in apure spatial perspective, especially for a moving surfacewith v n (cid:54) = 0. Unfortunately, most of the numerical toolsfor solving surface PDEs do not work with a locally de-fined vector basis. They mimic vector-valued problemsas a system of scalar-valued problems under the assump-tion of Euclidean coordinates. This means for the surface problem (B27) – (B29) that the tangential velocity fieldand the director field are considered to be vector fields in R , which results in an under-determined problem. Thetwo additional degrees of freedom can be handled in dif-ferent ways, e. g. by penalty methods (Nestler et al. , 2018;Reuther and Voigt, 2018) or by using Lagrange multipli-ers (Jankuhn et al. , 2017). The terms ∂ t v and ∂ t p makescertainly sense, if we consider v , p ∈ T R , but note thatin general ∂ t v as well as ∂ t p are no longer part of the tan-gential space of the surface S . Nevertheless, we can useonly the tangential part of ∂ t v = (cid:0) ∂ t v j (cid:1) ∂ j x + v j ∂ j ∂ t x for a transversal observer, i. e.[ ∂ t v ] i = (cid:104) ∂ t v , ∂ i x (cid:105) R = g ij ∂ t v j − v n B ij v j .Analogously, the same holds for ∂ t p . Finally, we obtainby rewriting eqs. (B27) – (B29) π S ∂ t v + ∇ S v v − v n ( B v + ∇ S v n )= ν (cid:0) − ∆ dR v + 2 κ v + ∇ S ( v n H ) − S ( v n B ) (cid:1) − ∇ S p S − λ div S σ E S (B30)div S v = v n H (B31) π S ∂ t p + ∇ S v p = η (cid:0) ∆ DG p − B p (cid:1) − ω n (cid:16) (cid:107) p (cid:107) S − (cid:17) p (B32)in S × R + , if v , p ∈ T S ⊂ T R .Eq. (B3) together with the boundary condition (B5)is the L -gradient flow along the material motion tominimize the Frank-Oseen energy functional with ma-terial constants η = K = K = K and K = 0,see (Mori et al. , 1999). In the thin film limit, thisleads to the minimization of the surface Frank-Oseenenergy η (cid:82) S (cid:107)∇ S p (cid:107) S d S . This situation differs from(Nestler et al. , 2018), where the one-constant approx-imation η = K = K = K = − K was as-sumed, which leads to minimizing the distortion en-ergy η (cid:82) S (rot S p ) + (div S p ) d S . However, for the case ω n → ∞ , where (cid:107) p (cid:107) S = 1 a. e., both energies only differby a constant value K (cid:82) S κ d S = πKχ ( S ), where χ ( S )denotes the Euler characteristic. Thus, the minimizers ofboth energies are equal. REFERENCES
Abraham, R, J. E. Marsden, and T. S. Ratiu (1988),
Man-ifolds, Tensor Analysis, and Applications , Applied Mathe-matical Sciences No. 75 (Springer).Alaimo, F, C. K¨ohler, and A. Voigt (2017), “Curvature con-trolled defect dynamics in topological active nematics,” Sci.Rep. , 5211.Arroyo, M, and A. DeSimone (2009), “Relaxation dynamicsof fluid membranes,” Phys. Rev. E , 031915.Bade, N D, R. Kamien, R. Assoian, and K. Stebe (2017),“Patterned Cell Alignment in Response to Macroscale Cur-vature,” Biophys. J. (3), 536a. Bates, M A, G. Skaˇcej, and C. Zannoni (2010), “Defectsand ordering in nematic coatings on uniaxial and biaxialcolloids,” Soft Matter (3), 655–663.Bertalmio, M, L.-T. Cheng, S. Osher, and G. Sapiro (2001),“Variational Problems and Partial Differential Equationson Implicit Surfaces,” J. Comput. Phys. , 759–780.Chaikin, P M, and T. C. Lubensky (1995), Principles ofcondensed matter physics (Cambridge University Press).Chorin, A J (1968), “Numerical solution of the Navier-Stokesequations,” Math. Comp. , 745–762.Dhakal, S, F. J. Solis, and M. Olvera de la Cruz (2012), “Ne-matic liquid crystals on spherical surfaces: Control of defectconfigurations by temperature, density, and rod shape,”Phys. Rev. E , 011709.Dziuk, G, and C. M. Elliott (2007a), “Finite elements onevolving surfaces,” IMA J. Numer. Anal. , 262–292.Dziuk, G, and C. M. Elliott (2007b), “Surface finite elementsfor parabolic equations,” J. Comput. Math. , 385–407.Dziuk, G, and C. M. Elliott (2008), “Eulerian finite elementmethod for parabolic PDEs on implicit surfaces,” InterfacesFree Bound. (1), 119–138.Dziuk, G, and C. M. Elliott (2013), “Finite element methodsfor surface PDEs,” Acta Numer. , 289–396.Dzubiella, J, M. Schmidt, and H. L¨owen (2000), “Topologicaldefects in nematic droplets of hard spherocylinders,” Phys.Rev. E , 5081–5091.Ebin, D G, and J. Marsden (1970), “Groups of Diffeomor-phisms and the Motion of an Incompressible Fluid,” Ann.Math. , 102–163.Ellis, P W, D. J. G. Pearce, Y.-W. Chang, G. Goldsztein,L. Giomi, and A. Fern´andez-Nieves (2018), “Curvature-induced defect unbinding and dynamics in active nematictoroids,” Nat. Phys. , 85–90.Ericksen, J L (1961), “Conservation Laws for Liquid Crys-tals,” T. Soc. Rheo. (1), 23–34.Ericksen, J L (1976), “Equilibrium Theory of Liquid Crys-tals,” in Advances in Liquid Crystals , edited by G. H.Brown (Elsevier) pp. 233–298.Fern´andez-Nieves, A, V. Vitelli, A. S. Utada, D. R. Link,M. M´arquez, D. R. Nelson, and D. A. Weitz (2007), “NovelDefect Structures in Nematic Liquid Crystal Shells,” Phys.Rev. Lett. , 157801.Freeden, W, and M. Schreiner (2008), Spherical Functions ofMathematical Geosciences: A Scalar, Vectorial, and Ten-sorial Setup , Advances in Geophysical and Environmen-tal Mechanics and Mathematics (Springer, Berlin, Heidel-berg).de Gennes, P G, and J. Prost (1993),
The Physics of LiquidCrystals (Second Edition, Oxford Science Publications).Greer, J B, A. L. Bertozzi, and G. Sapiro (2006), “Fourthorder partial differential equations on general geometries,”J. Comput. Phys. , 216–246.Gross, B J, and P. J. Atzberger (2018), “Hydrodynamic flowson curved surfaces: Spectral numerical methods for radialmanifold shapes,” J. Comput. Phys. , 663–689.Guerra, R E, C. P. Kelleher, A. D. Hollingsworth, and P. M.Chaikin (2018), “Freezing on a sphere,” Nature , 346.Hansbo, P, M. G. Larson, and K. Larsson (2016), “Anal-ysis of Finite Element Methods for Vector Laplacians onSurfaces,” arXiv:1610.06747.Hu, X, and H. Wu (2013), “Long-time dynamics of the nonho-mogeneous incompressible flow of nematic liquid crystals,”Commun. Math. Sci. (3), 779–806.Huang, J, F. Lin, and C. Wang (2014), “Regularity and Ex- istence of Global Solutions to the Ericksen–Leslie Systemin R2,” Commun. Math. Phys. (2), 805–850.Jankuhn, T, M. A. Olshanskii, and A. Reusken (2017), “In-compressible Fluid Problems on Embedded Surfaces: Mod-eling and Variational Formulations,” IGPM report .Keber, F C, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi,M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch(2014), “Topology and dynamics of active nematic vesi-cles,” Science (6201), 1135–1139.Koning, V, T. Lopez-Leon, A. Fern´andez-Nieves, andV. Vitelli (2013), “Bivalent defect configurations in inho-mogeneous nematic shells,” Soft Matter , 4993–5003.Kralj, S, R. Rosso, and E. G. Virga (2011), “Curvature con-trol of valence on nematic shells,” Soft Matter , 670–683.Kruse, K, J. F. Joanny, F. J¨ulicher, J. Prost, and K. Sekimoto(2004), “Asters, Vortices, and Rotating Spirals in ActiveGels of Polar Filaments,” Phys. Rev. Lett. , 078101.Leslie, F M (1968), “Some constitutive equations for liquidcrystals,” Arch. Ration. Mech. An. (4), 265–283.Liang, H-L, J. H. Noh, R. Zentel, P. Rudquist, and J. Lager-wall (2013), “Tuning the defect configurations in nematicand smectic liquid crystalline shells,” Phil. Trans. Roy. Soc.A (1988), 10.1098/rsta.2012.0258.Liang, H L, S. Schymura, P. Rudquist, and J. Lagerwall(2011), “Nematic-Smectic Transition under Confinement inLiquid Crystalline Colloidal Shells,” Phys. Rev. Lett. ,247801.Lin, F, J. Lin, and C. Wang (2010), “Liquid Crystal Flowsin Two Dimensions,” Arch. Ration. Mech. An. (1),297–336.Lin, F-H, and C. Liu (2000), “Existence of Solutions for theEricksen-Leslie System,” Arch. Ration. Mech. An. (2),135–156.Liu, I B, N. Sharifi-Mood, and K. J. Stebe (2018), “CapillaryAssembly of Colloids: Interactions on Planar and CurvedInterfaces,” Annu. Rev. Condens. Matter Phys. , 283–305.Lopez-Leon, T, A. Fern´andez-Nieves, M. Nobili, and C. Blanc(2012), “Smectic shells,” J. Phys.-Condens. Mat. (28),284122.Lubensky, T C, and J. Prost (1992), “Orientational orderand vesicle shape,” J. de Physique II (3), 371–382.Marth, W, S. Praetorius, and A. Voigt (2015), “A mechanismfor cell motility by active polar gels,” J. Roy. Soc. Interface (107), 10.1098/rsif.2015.0161.Martinez, A, M. Ravnik, B. Lucero, R. Visvanathan,S. Zumer, and I. I. Smalyukh (2014), “Mutually tangledcolloidal knots and induced defect loops in nematic fields,”Nat. Mater. , 258–263.Mickelin, O, J. S(cid:32)lomka, K. J. Burns, D. Lecoanet, G. M. Vasil,L. M. Faria, and J. Dunkel (2018), “Anomalous ChainedTurbulence in Actively Driven Flows on Spheres,” Phys.Rev. Lett. , 164503.Mitrea, M, and M. Taylor (2001), “Navier-Stokes equationson Lipschitz domains in Riemannian manifolds,” Math.Ann. , 955–987.Miura, T-H (2017), “On singular limit equations for incom-pressible fluids in moving thin domains,” arXiv:1703.09698.Mori, H, E. C. Gartland Jr., J. R. Kelly, and P. J. Bos (1999),“Multidimensional Director Modeling Using the Q TensorRepresentation in a Liquid Crystal Cell and Its Applicationto the pi Cell with Patterned Electrodes,” Jpn. J. Appl.Phys. (1R), 135.Napoli, G, and L. Vergori (2010), “Equilibrium of nematicvesicles,” J. Phys. A (44), 445207. Napoli, G, and L. Vergori (2012), “Extrinsic Curvature Ef-fects on Nematic Shells,” Phys. Rev. Lett. , 207803.Napoli, G, and L. Vergori (2016), “Hydrodynamic theory fornematic shells: The interplay among curvature, flow, andalignment,” Phys. Rev. E , 020701.Nelson, D R (2002), “Toward a Tetravalent Chemistry of Col-loids,” Nano Lett. (10), 1125–1129.Nestler, M, I. Nitschke, S. Praetorius, and A. Voigt (2018),“Orientational Order on Surfaces: The Coupling of Topol-ogy, Geometry, and Dynamics,” J. Nonlinear Sci. (1),147–191.Nitschke, I, M. Nestler, S. Praetorius, H. L¨owen, and A. Voigt(2018), “Nematic liquid crystals on curved surfaces – a thinfilm limit,” Proc. Roy. Soc. A , 20170686.Nitschke, I, S. Reuther, and A. Voigt (2017), “Discrete Ex-terior Calculus (EC) for the Surface Navier-Stokes Equa-tion,” in Transport Processes at Fluidic Interfaces , editedby D. Bothe and A. Reusken (Springer) pp. 177–197.Nitschke, I, and A. Voigt (2019), in preperation.Nitschke, I, A. Voigt, and J. Wensch (2012), “A finite elementapproach to incompressible two-phase flow on manifolds,”J. Fluid Mech. , 418–438.Olshanskii, M A, A. Quaini, A. Reusken, and V. Yushutin(2018), “A finite element method for the surface Stokesproblem,” arXiv:1801.06589.Praetorius, S, A. Voigt, R. Wittkowski, and H. L¨owen (2018),“Active crystals on a sphere,” Phys. Rev. E , 052615.Prinsen, P, and P. van der Schoot (2003), “Shape anddirector-field transformation of tactoids,” Phys. Rev. E ,021701.R¨atz, A, and A. Voigt (2006), “PDE’s on surfaces: A diffuseinterface approach,” Commun. Math. Sci. , 575–590.Reusken, A (2018), “Stream Function Formulation of SurfaceStokes Equations,” IGPM report .Reuther, S, and A. Voigt (2015), “The Interplay of Curva-ture and Vortices in Flow on Curved Surfaces,” MultiscaleModel. Sim. (2), 632–643.Reuther, S, and A. Voigt (2018), “Solving the incompressiblesurface Navier-Stokes equation by surface finite elements,”Phys. Fluids (1), 012107.Segatti, A, M. Snarski, and M. Veneroni (2014), “Equilibriumconfigurations of nematic liquid crystals on a torus,” Phys.Rev. E , 012501.Shin, H, M. J. Bowick, and X. Xing (2008), “Topological De-fects in Spherical Nematics,” Phys. Rev. Lett. , 037802.Shkoller, S (2002), “Well-posedness and global attractors forliquid crystals on Riemannian manifolds,” Commun. Part.Diff. Eq. (5-6), 1103–1137.Simha, R A, and S. Ramaswamy (2002), “HydrodynamicFluctuations and Instabilities in Ordered Suspensions ofSelf-Propelled Particles,” Phys. Rev. Lett. , 058101.Sknepnek, R, and S. Henkes (2015), “Active swarms on asphere,” Phys. Rev. E , 022306.Stark, H (2001), “Physics of colloidal dispersions in nematicliquid crystals,” Phys. Rep. , 387–474.St¨ocker, C, and A. Voigt (2008), “Geodesic Evolution Laws -A Level-Set Approach,” SIAM J. Imaging Sci. , 379–399.Tjhung, E, D. Marenduzzo, and M. E. Cates (2012), “Spon-taneous symmetry breaking in active droplets provides ageneric route to motility,” Proc. Natl. Acad. Sci. USA ,12381–12386.Vey, S, and A. Voigt (2007), “AMDiS: Adaptive multidimen-sional simulations,” Comput. Vis. Sci. , 57–67.Wang, W, P. Zhang, and Z. Zhang (2013), “Well-Posedness of the Ericksen–Leslie System,” Arch. Ration. Mech. An. (3), 837–855.Witkowski, T, S. Ling, S. Praetorius, and A. Voigt (2015),“Software concepts and numerical algorithms for a scalableadaptive parallel finite element method,” Adv. Comput.Math. , 1145–1177.Yang, Z, J. Wei, Y. I. Sobolev, and B. A. Grzybowski (2018),“Systems of mechanized and reactive droplets powered bymulti-responsive surfactants,” Nature , 313.Yavari, A, A. Ozakin, and S. Sadik (2016), “Nonlinear Elas-ticity in a Deforming Ambient Space,” J. Nonlinear Sci. (6), 1651–1692.Zhou, W, J. Cao, W. Liu, and S. Stoyanov (2008), “HowRigid Rods Self-Assemble at Curved Surfaces,” Angew.Chem. Int. Ed.48