Analysis and Numerical Approximation of Energy-Variational Solutions to the Ericksen–Leslie Equations

We define the concept of energy-variational solutions for the Ericksen–Leslie equations in three spatial dimensions. This solution concept is finer than dissipative solutions and satisfies the weak-strong uniqueness property. For a certain choice of the regularity weight, the existence of energy-variational solutions implies the existence of measure-valued solutions and for a different choice, we construct an energy-variational solution with the help of an implementable, structure-inheriting space-time discretization. Computational studies are performed in order to provide some evidence of the applicability of the proposed algorithm.


Introduction
Liquid crystals comprise the structural properties of crystals within a fluid.The fluid flow of the nematic phase of liquid crystals can be described by the Ericksen-Leslie system.In this model, the material behaves like a liquid, i.e., no positional order is present, but the molecules exhibit a long-range self-alignment along a direction (see Figure 1).In this way, liquid crystals sustain Figure 1: Alignment of the rod-like molecules in an isotropic liquid, the nematic phase of a liquid crystal, and a solid, [24,Fig. 1.1] anisotropic dynamics such as the polarization of light or the transfer of heat, but at the same time offer the physical flexibility of a fluid, which makes these materials interesting for engineering and sciences [23].
Ericksen [16] and Leslie [31] derived the Ericksen-Leslie equations during their development of an instationary theory of liquid crystals in the 1960s.Let v v v : Ω × [0, T ] → R 3 denote the velocity of the fluid, p p p : Ω × [0, T ] → R its pressure and d d d : Ω × [0, T ] → R 3 the director.We consider the system governed by the equations with µ 4 > 0, µ 5 + µ 6 − λ 2 ≥ 0, µ 1 + λ 2 ≥ 0 in order to ensure the dissipative character of our model.So far, a vast majority of the mathematical work on the Ericksen-Leslie model considers a simplified system with a relaxed unit-norm constraint that is only enforced approximately by adding a Ginzburg-Landau penalization term f f f ε (d) = 1 4ε (|d d d| 2 − 1) 2 to the free energy potential 1 2 |∇d d d| 2 .With simplified Leslie-stress tensor, the momentum and director equation (1) are replaced by respectively.In the first analysis of this system [33], the authors were able to prove the existence of weak solutions to (2).A rather general model including the full Leslie stress tensor is considered by [10], where again a Ginzburg-Landau penalization approach is introduced to replace the unit-norm constraint of the director.In this setting the authors prove that weak solutions exist and a blow-up criterion for local strong solutions.In [15] the existence of weak solutions is generalized to a larger class of free energy functions.An overview of the analytical results regarding the Ericksen-Leslie equations and its connection to other models for liquid crystals can be found in [14].In two spacial dimensions, for the limiting system of (2), where f ε (d d d) is replaced by −|∇d d d| 2 , it is known that a unique weak solution exists that is smooth except for finitely many points in time [34].
For a general model of the Ericksen-Leslie equations equipped with the naturally arising anisotropic Oseen-Frank energy, the concept of dissipative solutions is applied in [26], which are shown to be the local average of measure-valued solutions [27] and inherit their weak-strong uniqueness [28].In comparison to measure-valued solutions, dissipative solutions have the advantage that they have less degrees of freedom and can be approximated by numerical schemes.Nevertheless, they form only a subset of measure-valued solutions and are thus not as precise.In the work at hand, we introduce energy-variational solutions (cf.[30]), which have one degree of freedom more than dissipative solutions and can be argued to be as fine as measure-valued solutions, but as we will show, they can also be approximated by numerical schemes and have additional advantages.This solution concept is useful where one might either not be able to derive weak solutions for physically relevant models or where weak solutions admit unphysical behaviour, like unphysical non-uniqueness [22].The solution concept was first introduced in [29].We will prove that this concept is finer than the so-called dissipative solutions introduced by Lions [36] for the Euler equations and also applied to the Ericksen-Leslie equations by one of the authors [26].Additionally, in certain scenarios this concept is finer than the concept of measure-valued solutions.Energy-variational solutions do not only fulfill the standard weak-strong uniqueness property of generalized solution concepts (cf.[28]), but they also fulfill the semi-flow property such that prolongations and restrictions of solutions on larger, and smaller time intervals, respectively, are energy-variational solutions again.In [29] it was also argued that energy-variational solutions are amenable for different selection criteria.The set of energy-variational solutions can be seen as a convex, weakly * -closed superset of weak solutions but a subset of dissipative solutions and for a special choice of the regularity weight also a subset of measure-valued solutions.This may allow to introduce techniques from optimization theory in order to select the physically relevant solution maximizing the dissipation in every point-in-time [29].
In this work, we define the concept of energy-variational solutions for the Ericksen-Leslie equations in three spatial dimensions.This definition has some freedom, since it depends on the choice of a certain regularity weight K .For one choice of such a regularity weight, we prove the equivalence to a certain class of measure-valued solutions, and for another choice, we construct an energyvariational solution with the help of an implementable, structure-inheriting space-time discretization based on the finite element method.
Concerning the numerical approximation of the Ericksen-Leslie equations, a first study for the simplified system (2) equipped with the Ginzburg-Landau approximation is conducted in [37].They combined an implicit Euler scheme in time with Hermite type finite elements for the director and Q2-Q1-Taylor-Hood elements for the velocity and pressure.In their subsequent work [38], the authors replace the demanding Hermite finite elements for the director by piecewise quadratic functions.Even a relaxation from C 1 finite elements to C 0 finite elements is realized in [35].A different approach is proposed in [3], where the Ginzburg-Landau approximation of the unit-norm constraint is interpreted as a saddle-point structure.For the simplified and penalized Ericksen-Leslie system also decoupling techniques and mixed methods are examined in [18,9].In [5] two numerical schemes for the simplified model are proposed.The first one uses the Ginzburg-Landau approximation for the unit-norm constraint and the second proposed scheme does not depend on a regularization parameter and fulfills the unit-norm constraint in the limit.However both schemes do not fulfill the unit-norm constraint at the discrete level exactly.We therefore use an approach from [4] by implementing a midpoint rule at the finite element level which solves the sphere constraint exactly at every node of the mesh.But we have to refine this approach by introducing a special projection (see Remark 2.5), which will allow to identify the limit as an energy-variational solution.The proposed scheme is the first numerical scheme implementing the main properties of the continuous system including the algebraic norm restriction at every node of the mesh such that the approximate solutions converge to an energy-variational solution fulfilling the physically relevant semi-flow property.We think that the concept of energy-variational solutions is a strong tool for identifying limits of solutions to numerical schemes and can also serve as such in other related models.
After providing the considered system and an overview over the existing literature, we introduce the necessary notation in Section 2, we provide the definitions and the main results in Subsection 2.1.The weak-strong uniqueness proof and the relation to measure-valued solutions in considered in Subsection 2.2, whereas the preliminaries are provided in the following subsections on auxiliary lemmata 2.3, finite elements 2.4, and interpolation 2.5.The discrete system is introduced and analysed in Section 3, we provide the necessary a priori estimates 3.1, extract converging subsequences 3.2, and prove convergence to the director equations 3.3 and the energy-variational inequality 3.4.Finally, some computational studies are presented in Section 4 in order to show some evidence of the applicability of the proposed algorithm.
2 Preliminaries, main results, and continuous system We denote the space of smooth solenoidal functions with compact support by C ∞ c,σ (Ω; R 3 ).By L p σ (Ω), H 1 0,σ (Ω), and W 1,p 0,σ (Ω), we denote the closure of C ∞ c,σ (Ω; R 3 ) with respect to the norm of L p (Ω), H 1 (Ω), and W 1,p (Ω), respectively.By H 1 (Ω; S 2 ), we denote the functions d d d ∈ H 1 (Ω) such that |d d d| = 1 a.e. in Ω.The L 2 (Ω) inner product is thereby denoted by (., .).The Dual space of a Banach space V is denoted by V * , where the dual pairing is denoted as •, • .In order to define the cross product, we use the Levi-Cita symbol.The cross product of two vectors x, y ∈ R 3 is then defined as (x × y) i = ∑ j,k ε i, j,k x j y k , the cross product of a vector x ∈ R 3 with a matrix A ∈ R d×d as (x × A) i, j = ∑ l,m ε i,l,m x l A m, j .We further will make use of the matrix notation of the cross product for three dimensions, i.e. ([x] × ) ik = ∑ j ε i, j,k x j .We further introduce the discrete derivative d t as for a constant time-step size k > 0 and for a sequence in a normed space ( f j ) j=1,...,n , where we set d t f 0 = 0. Usually functions in the continuous setting are denoted by bold letters in contrast to the discrete functions.
By R d×d we denote d-dimensional quadratic matrices, by R d×d sym the symmetric subset, and by R d×d sym,+ the symmetric positive semi-definite matrices.The symmetric and skew-symmetric part of a matrix A ∈ R d×d are denoted by (A) sym and (A) skw , respectively.The positive and negative semidefinite part of a matrix A ∈ R d×d are denoted by (A) + and (A) − , respectively.We equip the last set with the usual spectral norm |A| 2 = max i∈{1,...,d} λ i , where λ i are the nonnegative eigenvalues of the matrix A ∈ R d×d sym,+ .The dual norm of the spectral norm with respect to the Frobenius product (A : The Radon measures taking values in R d×d sym are denoted by M (Ω; R d×d sym ), which may be interpreted as the dual space of the continuous functions,i.e., M (Ω; R d×d sym ) = (C (Ω; R d×d sym )) * .Note that an element µ ∈ M (Ω; M d×d sym,+ ) is a Radon measure taking values in the symmetric matrices such that for any ξ ξ ξ ∈ R d the measure ξ ξ ξ ⊗ ξ ξ ξ : µ is nonnegative.By I, we denote the identity matrix in R d×d .
For a given Banach space X, the space C w ([0, T ]; X) denotes the functions on [0, T ] taking values in X that are continuous with respect to the weak topology of X.The space L ∞ w * ([0, T ]; X * ) is the space of all function on [0, T ] taking values in X * that are Bochner measurable with respect to X * equipped with the weak-stark topology and essentially bounded.The total variation of a function where the supremum is taken over all finite partitions of the interval [0, T ].We denote the space of all bounded functions of bounded variations on [0, T ] by BV([0, T]).Note that the total variation of a monotone decreasing nonnegative function E only depends on the initial value, i.e.,

Main results
We define the total energy as The main new idea when defining energy-variational solutions is to introduce an auxiliary variable E as an upper bound of the total energy E (v v v,d d d) and add the difference of these two variables to a variational inequality of the energy-dissipation mechanism to close this formulation with respect to the appropriate weak topologies.The difference between the variable E and the energy E (v v v,d d d) can be interpreted as a measure of the difference between weak and strong convergence of the approximate solutions.The associated error term represents this difference in the limit of vanishing discretization parameters.
Definition 2.1 (Energy-variational solutions).We call for all s, t ∈ (0, T ) and for all ṽ a.e. in Ω × (0, T ).In this work, we choose , where c L is given in (22), with the associated Ericksen stresses being defined as Remark (Ericksen stress).The main obstacle when passing to the limit in an approximation of the Ericksen-Leslie equations is often the so-called Ericksen stress, which can be seen as a kind of Korteweg stress coupling the elastic deformations of the director field to the evolution of the velocity field.This term is the main reason that we have to consider a generalized solution framework different from usual weak solutions.Via Lemma 2.6 we observe that both formulations T T T E 1 and T T T E 2 are formally equivalent via the reformulation where the second term on the right-hand side vanish formally due to the fact that ṽ v v ∈ H 1 0,σ (Ω) and |d d d| = 1.Thus, in the case of a regular solution with d d d ∈ L 2 (0, T ; H 2 (Ω; S 2 )) and E = E (v v v,d d d) a.e. in (0, T ) both solution concepts coincide with usual weak solutions.Remark (Strong continuity of the initial value).From the monotony of t → E(t) and the weak continuity of t → (v v v(t),d d d(t)), we infer by the lower semi-continuity of the energy functional and the inequality E(t) ≥ E (v v v(t),d d d(t)) for all t ∈ [0, T ] that for any sequence t 0 it holds This implies that all inequalities are in-fact equations such that we infer from the uniform convexity of L 2 (Ω) that ) .Remark (Properties of generalized solutions).Energy-variational solutions fulfill several properties, which are desirable for generalized solution concepts.That generalized solutions coincide with classical solutions, if the latter exists, is the so-called weak-strong uniqueness property and covered by Theorem 2.5.
Furthermore, energy-variational solutions fulfill the so-called semi-flow property such that restrictions and concatinations of solutions are solutions again, which is an important property of a reasonable solution concept (cf.[30]).Remark (Different definition).We could also identify for a.e.t ∈ (0, T ) in the inequality (5).This can be done due to the strong convergence of v v v a.e. in Ω × (0, T ).We would also replace The disadvantage would be that the adapted inequality (5) is only fulfilled for a.e.s ∈ (0, T ) and all t ∈ (s, T ]. We are going to prove the theorem by the convergence of a fully discrete, implementable, unconditionally solvable scheme in the following sections.The proposed scheme in Section 3 is a fully implementable numerical scheme for the Ericksen-Leslie system that fulfills the norm restriction |d d d| = 1 for every approximate solution in every node of the mesh.We even infer that the norm of the approximate solutions converges to 1 with the order 1 (cf.Proposition 3.4).Note that the condition d d d Γ = tr(d d d 0 ) is a usual compatibility condition for parabolic problems.
Energy-variational solutions can be seen as a generalized solution concept, which generalizes the well-known weak solutions, but is finer than the usual concept of so-called measure-valued solutions [27].
)), and if ( 6) is fulfilled as well as the measure-valued formulation of the Navier-Stokes equation for all ṽ v v ∈ C ∞ c (Ω × (0, T )) and the energy-inequality for a.e.t ∈ (0, T ).d d d, E) be an energy-variational solution with Then there exists a measure-valued solution d d d, E) be an energy-variational solution in the sense of Definition 2.1 with such that K (ṽ v v) ∈ L 1 (0, T ).Then the relative energy inequality holds for a.e.t ∈ (0, T ), where ) is a classical solution to the Ericksen-Leslie system on (0, T * ), then both solutions coincide on this time-interval, i.e., e. on (0, T * ) .
Remark (Dissipative solutions).Theorem 2.5 does not only show the weak-strong uniqueness of energy-variational solutions, it also shows that every energy-variational solution is a dissipative solution.A dissipative solution is defined via the relative-energy inequality (10), i.e., (v v v,d d d) is called a dissipative solution, if it fulfills the relative energy-inequality 10 for all function (ṽ v v, d d d) fulfilling the regularity properties (9).The concept of energy-variational solutions is thus strictly finer than the concept of dissipative solutions.d d d, E) be an energy-variational solution according to Definition 2.1 with K (ṽ v v) = 2 ∇ṽ v v C (Ω) .Then, we may choose ṽ v v = 0 in order to infer

Proofs for the continuous system
Furthermore, we may choose ṽ v v = α ũ u u with α ≥ 0 and multiply the energy-variational inequality (5) by 1/α.In the limit α → ∞, we find . By the definition of the weak time derivative, we infer that First this is only well defined for ũ , we infer from the Hahn-Banach theorem [8, Thm.1.1] the existence of an element This allows to write the linear form for a.e.t ∈ (0, T ).This implies that the linear form l l l ∈ L ∞ w * (0, T ; (C 1 0 (Ω)) * ) can be identified via considering the gradient operator ) the range of the gradient, the set of all gradient fields equipped with the norm of L 1 (0, T ; C (Ω; R 3×3 )).Note that this mapping is injective due to the prescribed Dirichlet boundary conditions.We may set for all A A A ∈ C (Ω; R 3×3 ) a.e. in (0, T ), where the estimate follows from (12) and Hahn-Banach's theorem written in its point-wise a.e.-form.By the Riesz representation theorem, we know that there exists a measure which implies the formulation ( 7) by the definition of l l l in (12).From the estimate (13), we infer that the measure m m m vanishes for all continuous functions with values in the skew-symmetric matrices, such that m m m ∈ L ∞ w (0, T ; M (Ω; R 3×3 sym )).Additionally, (13) allows to insert From the estimate in ( 13), we additionally infer by choosing for a.e.t ∈ (0, T ), which in turn implies together with (11) the inequality (8).Note that we use the usual spectral-norm for symmetric matrices |A| 2 = max i∈{1,2,3} |λ i |, where λ i denote the eigenvalues of the matrix A ∈ R 3×3 sym .The dual-norm with respect to the Frobenius product : is given by For a positive definite matrix all eigenvalues are non-negative, such that we may write |A| 2 = ∑ 3 i=1 λ i = tr(A) = A : I.This implies that we may identify the norm |A| 2 = I : A.
Proof of Theorem 2.5.We observe that (6) tested with ∆ d d d and integrated-in-time gives Additionally, the equation (1a) for the strong solution tested with ṽ v v−v v v and (1c) for the strong solution tested with ∆(d We observe that Adding up ( 5), (14), and (15) implies We recall the identity, which is a special case of the identity (6.7) in [28], this identity helps to transform the second line on the right-hand side of ( 16) in the case of In case of T T T E = T T T E 1 we may observe from the uni-norm restriction and additional manipulations that The resulting terms can now be estimated in terms of the relative energy and the relative dissipation, i.e., the terms on the left-hand side.Similar arguments have already been applied in [28] and [25].
We remark, that in the term in the fifth line on the right-hand side of ( 16), we have to apply the product rule to infer In order to prove the weak-strong uniqueness result, the terms on the right-hand side of ( 16) have to estimated by terms on the left-hand side and in the end Gronwall's inequality is applied.Since most of the estimates are rather standard and done using Hölder's and Young's inequality, we refrain from displaying it in full detail here.The interested reader may consult the article [28] for more details, where these estimates were performed for a more involved system.
Estimating all terms appropriately, we find that there exists a constant C > 0 such that Gronwall's inequality implies the relative energy inequality (10).Choosing (ṽ v v, d d d) to be a solution, implies that the fourth and fifth line of (10) vanish.Since all terms on the left-hand side of (10) for a.e.t ∈ (0, T * ), which implies the assertion of Theorem 2.5.
Remark (Weak-strong uniqueness extended).It may seem strange that all calculations in the above proof are conducted on the time interval (0,t), when the definition of energy-variational solutions 2.1 is formulated on every interval (t, s) for all s < t ∈ (0, T ).Indeed, the calculations of the above proof can also be carried out for all intervals (s,t) such that we may infer inequality (17) with 0 replaced by s.Thus, the weak-strong uniqueness result can be sharpened: If there exists an s ∈ (0, T ) such that there exists a solution (ṽ v v, d d d) fulfilling the regularity assumptions (9) on (s, Remark.We note that the notation ∆d d d + |∇d d d| 2 d d d is the usual formulation in the Euler-Lagrange equation of harmonic maps into the sphere [19].The norm of this term appears as a dissipative term in the energy-dissipation mechanism of the system (see (5) below).Formally this implies that the solution to such a system is forced to approximately fulfill the Euler-Lagrange equation of the Dirichlet energy minimized over all functions with values in the sphere, in the large-time limit.

Auxiliary Lemmata
Proof.For all ϕ ϕ ϕ ∈ H 1 (Ω; S 2 ), we find by an integration-by-parts and ∇|d d d| 2 = 0 a.e. in Ω as well as a.e. on ∂ Ω that This implies the assertion.

Finite element spaces
For the rest of this paper we presume the following assumptions.
2. The subdivision T h of tetrahedra of our domain Ω ⊂ R 3 is quasi-uniform (in the sense of [7, 4.4.15]).
3. All T ∈ T h have at least one node that is not on the boundary ∂ Ω.
4. Every T ∈ T h is equipped with a Lagrangian finite element.
In order to apply well known results for the numerical approximation of incompressible flows, we choose (P2-P1) Taylor-Hood elements (cf.[21]) for the velocity.By X h and M h we denote the finite element spaces for the velocity and pressure, where P k (Ω) is the set of polynomials of degree k or less on the domain Ω.For better readability, we denote approximately divergence-free functions in X h as the space space V h , i.e.
It is a well-known result (cf.[49]) that this choice of spaces fulfills the inf-sup condition (sometimes also called LBB or Ladyzhenskaya-Babuska-Brezzi condition) under our given assumptions (1-4), that is sup for all q h ∈ M h and for some constant C > 0. The director and discrete Laplacian are approximated using piecewise linear functions, The homogeneous Dirichlet boundary conditions can be implied in the space by Our final mixed finite element space will be denoted by, Lemma 2.8 (Inverse estimate, cf.[7,Thm. 4.5.11]).Under the assumptions (1-4), let p, q ∈ [1, ∞] and 0 ≤ m ≤ l ≤ 1.Then there exists a constant C independent of h, such that v W l,p (Ω) ≤ Ch m−l+min{0,3/p−3/q} v W m,q (Ω) For the finite element spaces [Y h ] 3 and V h and a function f ∈ L 2 (Ω) the standard L 2 -projection, is denoted by →V h , respectively.Lemma 2.9 (Projection operator, cf.[50], [17,Lem. 1.131,Prop. 1.134]).Under the previously given assumptions (1 -4), the L 2 projection onto [Y h ] 3 is stable in L p (Ω) for all 1 ≤ p ≤ ∞ and H 1 (Ω), i.e. we have For the projection onto the finite element space for the velocity V h , Q h admits the error and stability estimates for all v v v ∈ H 1 0,σ (Ω) and ṽ v v ∈ H 2 (Ω) ∩ H 1 0,σ (Ω).This follows from [20,Lem. 4.3] and the regularity of the Stokes problem on convex polygons (see [13], see also [42,Proposition 2] or [40,Cor. 1.8]).

Interpolation and mass-lumping
We define the nodal interpolation operator for all functions f ∈ C(Ω), where N h is the set of all nodes in our mesh.For continuous functions y 1 , y 2 ∈ C(Ω; R 3 ) we introduce the lumped inner product (cf.[44,11]) and the according norm as For the interpolation operator we observe the following standard error estimate: Lemma 2.10.Under the assumptions (1-4), let 3/2 < p ≤ ∞.Then there exists a constant C > 0 independent of h and the function f , such that holds for all 0 ≤ k ≤ 2 and f ∈ W 2,p (Ω).
The above theorem is standard and can for example be found in [7,Thm. 4.4.20].There exists a generic constant c L > 0 independent of h > 0 such that for functions y 1 , y 2 ∈ [Z h ] 3 the estimates hold.
Proof.We proceed as in [11] and find for y 1 , The first inequality is a triangle inequality, the second one uses Hölder's inequality, the third bounds the interpolation error by Lemma 2.10.Since [Z h ] 3 consists of piecewise affine-linear functions the second derivative of a function in [Z h ] 3 vanishes.This together with Hölder's inequality, local and global inverse estimates (cf.2.8) allow to conclude.Note that the choice of exponents in the third step is necessary in order to make sure that we can indeed apply Lemma 2.10.
In case that one function is not a finite element functions, we have to estimate the difference with the interpolated function, i.e.,

|(I
and estimate the last factor on the right-hand side of ( 24) by We introduce another interpolation operator We observe the identity Lemma 2.12.Let f ∈ C (Ω).Then it holds Proof.Consider any y ∈ Ω.We obtain the pointwise estimate where the last term converges to 0 as h → 0, since g is continuous.Above, we used that ∑ z φ z (y) = 1 and that support of the shape functions φ z vanishes as h → 0. Since y was arbitrary the assertion follows.
Further, we introduce the discrete Laplacian in a slightly different way than usual, since we further equip it with mass lumping.Then the discrete Laplacian ∆ h : Similarly, we define an adapted projection via R h : This choice allows to infer with ( 22) that which implies together with (22 Remark (Adapted discretization of convection term).In comparison to the discretization in [4], we use the above projection R h in the definition of the convection term in comparison to the L 2 -projection P h .This change is needed in order to show the convergence to the finer concept of energy-variational solutions.In [4] only the convergence to dissipative solutions was shown, which is a bigger class of solutions.One key ingredient for this improvement are the above normequivalences.

Discrete system
In order to restrict ourselves to finite element spaces with a zero trace, we decompose the spatial approximation of our director into its interior and its non-homogeneous Dirichlet boundary condition, i.e.
with tr(d for all (a, c) ∈ U h .Hereby we discretize the Leslie stress tensor as where T j D collects the dissipative terms of the discrete Leslie stress tensor again, i.e.
Remark (Features of the discrete system).The mass lumping is applied on the director equation (32b) to guarantee the unit-norm restriction of d j in every node.In order to find a discrete energy law, the associated terms in the momentum equation (32a) and the discrete Laplacian are mass-lumped as well.Since the gradient of the director is piecewise constant and therefore not well-defined at the nodes of the mesh, we apply the mass-lumped L 2 projection R h in the convection term.Note that this is computationally not too expensive, since it only effects the gradient of the past iterate.Additionally, we note that solving the discrete system amounts in introducing two additional coupled equations, one to calculate the discrete Laplacian (28) and one to assure that the solution v j fulfills the discrete vanishing divergence condition (19).

A priori estimates
The proposed scheme is unconditionally solvable.Proposition 3.1 (Existence of discrete solutions).Let k, h > 0, j ∈ N and 1 ≤ j ≤ J = T /k .Then there exists a solution u j = (v j , d j • ) ∈ U h solving Scheme (32).
Proof.We show that the map Π : U h → U h implicitly defined by the above scheme (32), whose zero is the next iterate in terms of u = (v j , −∆ h d j−1/2 ), fulfills for all u U h ≥ R j−1 for some constant R j−1 > 0.Where the dependency on d j has to be interpreted via the discrete Greens function G h associated to the discrete Laplacian and given by Noting that this operator is a continuous bijection due to Lax-Milgram, we observe that the associated mapping Π corresponding to our discrete scheme is actually well defined and continuous.Then the existence of solutions to our discrete scheme (32) follows from a non-linear version of Brouwer's fix-point theorem (cf.[47]).In particular, we evaluate where we used the discrete identities Equation ( 35) is indeed positive, if holds, which only depends on u j−1 , k, h.This is implied by choosing an appropriate norm, since we have Note that the second term stems from a simple claculation obtained by applying the inverse estimate: . Proposition 3.2 (Properties of the discrete system).Let u j = (v j , d j • ) ∈ U h be a solution of the discrete scheme (32) for all 1 ≤ j ≤ n.Then the follwing energy equality holds, for all 1 ≤ n ≤ J.The algebraic restriction is fulfilled in every node where N h is the set of all nodes of the mesh.Additionally, it holds for some constant C > 0 that Proof.The discrete energy estimate (37) follows simply from iterating (35).In order to show (38), we adapt the approach in [4].Therefore, the director equation (32b) is tested with d j−1/2 (ẑ)φ ẑ, where ẑ ∈ N h is a node of the mesh.All terms except for the one including the time-derivative of the director vanish due to the mass lumping and the skew-symmetry.Then, we are left with which yields the result for all interior nodes ẑ.At the boundary the unit-norm constraint is fulfilled by the assumption on d d d 0 , our inhomogeneous Dirichlet boundary condition.Korn's inequality yields the estimate for the last addend in (39).Since d t d j ∈ Y h , the definition of the L 2 -projection yields for a smooth test function Using a duality argument, we are able to estimate the time derivative of the director in a L p -norm by The first term can easily be estimated using our discrete director equation (32b) for c ∈ Y h : Choosing c = P h (φ φ φ ) for φ φ φ ∈ H 1 (Ω) and using the embedding H 1 (Ω) → L 6 (Ω), the stability of the projection, Lemma 2.9, and the a priori estimates of Proposition 3.2 allow to obtain a bound for the right-hand side.The second term in (40) describes the error introduced by the mass-lumping.It must be estimated locally as in (24), i.e. by additionally applying the inverse estimate one obtains For the time derivative of the director d t d j in the L 2 -norm, we observe locally by choosing c = ∑ z z z∈N h ∩K d t d j (z z z)φ z denotes the in a local version of (41) that With the inverse estimate we can also bound the time derivative in the L 3 -norm by Putting the above together and reinserting it into (42) yields Multiplying by k and summing up leads to the estimate of d t d j , due to the a priori estimates (37).
For the velocity, we proceed as in [15]: We test the time derivative of the velocity with a smooth solenoidal test function φ φ φ , and apply the projection Q h such that d t v j ,φ φ φ = (d t v j , Q h φ φ φ ) and equation (32a), The convective terms can be estimated by a constant using (37) and For the Leslie stress tensor we infer the estimate Note that we hereby reformulated the term including the skew-symmetric part of the velocity's gradient on the algebraic level, i.e.
The Ericksen stress tensor can be estimated in a similar fashion Summing over j, using the embedding H 2 (Ω) → L ∞ (Ω) and the stability of the projection onto V h due to (20) yields the assertion.
) ∈ U h be a solution to scheme (32).Then the discrete energyvariational inequality with the energy E j and the regularity weight K given by holds for all ṽ (37), we find the variational-energy inequality By defining the variable E j and adding

Converging subsequences
Let {u j } j=0,...,J be a sequence of measurable functions in space.We define their constant and linear interpolates in time as follows, This yields the standard relation between the continuous and discrete time derivative, Analogously we define the interpolates for a smooth test func- The a priori estimates (37) allow us to extract converging subsequences which we do not relabel and (v v v,d d d, E) fulfilling ( 4) such that where the last convergence follows from Helly's selection theorem (cf.[8,Ex. 8.3]).We can refine the convergence for the director by using the following compact embeddings, referred to as Aubin-Lions-Simon Lemma [48], which yield strong convergence for the director A standard interpolation inequality (cf.[6, p. 192 ff.]) and the uniform bound in L ∞ (0, T ; L ∞ (Ω)) yields strong convergence, For the term d d d × ∆d d d, note that due to the above convergences and ( 23), we have ).Now note that the first term vanishes as h → 0 and the last term converges to has to be understood in the distributional sense.The fact that all interpolates in (44) converge against the same limit can be confirmed as in [43].Consider examplary for the velocity, which vanishes for k → 0, since the term on the right-hand side is bounded by our energy estimate (37).In order to deduce that the velocity is solenoidal, we use standard arguments (e.g. as in [12]).We test the divergence of the velocity with a test function φ φ φ ∈ C ∞ c (Ω × (0, T ]) with Ω φ φ φ dx = 0 and use our discrete divergence-zero condition to obtain where Π h is the L 2 -projection onto M h .The first term vanishes due to the weak convergence of the velocity in (44).The second summand vanishes for h → 0 simply due to the uniform boundedness v h k and a standard approximation property of the space M h (see e.g.[7, Lemma 12.4.3]),by The convergence for our non-homogeneous Dirichlet boundary conditions of the director, follows from the trace theorem (cf.[2, Thm.6.3.10]) and the interpolation error estimate (2.10), for every t ∈ [0, T ].Analogously, we observe for the initial condition of the director For the initial condition of the velocity, we can derive an analogous estimate -however in the L 2 norm instead -by making use of ( 20), since we approximated the initial condition with the respective projection operator, v 0 = Q h (v 0 ).

Unit-norm restriction and convergence for the director equations
From (38), we can infer that the unit-norm restriction |d(t, x)| = 1 holds at every node z ∈ N h .In the limit as h → 0, this restriction will be fulfilled almost everywhere on Ω, since I h (|d j | 2 − 1) = 0 holds.By subsequently applying the interpolation error bound (cf.Lemma 2.10) and the inverse estimate (cf.Lemma 2.8), one obtains where the right-hand side is bounded uniformly due to our a priori energy estimate (37).This gives us the result of Proposition 3.4.
Proposition 3.4.Let d k h be the approximate solution in the sense of the discrete scheme (32) that converges to a limit as in (44).Then the norm of d k h converges to 1 in a linear order, i.e.
In order to prove the convergence of the approximate solution to (32b) to a solution of (6), we consider equation (32b) tested with the projection P h c c c of a smooth test function c c c ∈ C ∞ c (Ω × (0, T )) and observe, for instance by Lemma 2.11 that This term vanishes as h → 0 since the last factor is bounded for smooth c c c due to 2.9.For the convection term, we infer due to the definition 29 For the last term, we infer such that we find The first inequality is due to Lemma 2.10, in the second inequality, we calculate the derivative, where we used that the second derivative of P1 finite elements vanishes on every tetrahedron.In order to infer the last inequality, we use Hölder's inequality and inverse estimates.
In order to pass from the mass-lumped inner product to the standard L 2 (Ω) inner product, other occurring terms have to be estimated in the same fashion, but locally.Regarding the second term for example, we observe , where we used that the second derivatives of the Lagrangian finite element functions vanish on the set K. Using Hölder's inequality and inverse estimates for all appearing terms, we infer .
The right-hand side is bounded due to the H 1 -stability of the L 2 projection (cf.Lemma 2.9).The estimates for the difference of the mass-lumped and L 2 (Ω) inner product for the other terms follow analogously.Thus, it remains to pass to the limit in (32b) with L 2 (Ω)inner products instead of the lumped-inner products.Note that by standard approximation arguments Lemma 2.9.The weak convergences of ( 44) and ( 45) allow to pass to the limit in all occurring terms, and we infer it also holds a.e. in Ω × (0, T ), which is nothing else than (6).

Convergence to the energy-variational formulation
Multiplying the discrete energy-variational formulation (43) by φ φ φ j for a φ φ φ ∈ C ∞ c (0, T ) applying an discrete integration as well as a discrete integration-by-parts, we find Since φ φ φ ≥ 0 for all t ∈ [0, T ], the second term can be estimated from below by zero.We observe that the terms T D k h ; (∇v k h ) sym are non-negative and quadratic in ∇v k h such that they are convex and weakly-lower continuous.Furthermore, we observe that , and from from inequality (30 h such that we may find By writing a a a 2 L ∞ (Ω) I − a a a ⊗ a a a = a a a for a a a ∈ L ∞ (Ω), we may use (22) to infer .
The first equality holds due to the definition of the interpolation operator in (27).The second equality is a transformation and the third is the definition of (29).In the fourth equation, a zero is added and the last inequality follows from Hölder's inequality.Furthermore, we estimate the second term on the right-hand side of the previous inequality by .
The first inequality follows from the local estimate of the interpolation error (cf.Lemma 2.10) and the quasi-uniformity of the mesh.In order to infer the second inequality, we calculate the derivatives via the product rule and observe that the second derivative vanishes for all function in [Z h ] 3 .In the third estimate, we use Hölder's inequality and local inverse estimates.We note that the interpolation I h converges in C (Ω) (cf.Lemma 2.12).Moreover, from (20) and Gagliardo-Nirenberg's inequality [41] we infer ṽ for all ṽ v v ∈ H 2 (Ω) ∩ H 1 0,σ (Ω).The weak convergence (44) and the strong convergence (47) implies the weak convergence In a similar fashion, we may infer the convergences The pointwise convergence of E k h → E allows to pass to the limit in the occuring terms, even thought the prefactors are possibly negative.In all other terms, the quantities in (50) (∇v )) only occur linearly and can thus, be handled in a standard fashion.The switch from masslumping to the L 2 inner product and the convergence of the projection of the gradient can be dealt with as in the previous section.

Computational studies
An efficient implementation of our discrete scheme (32) faces two challenges: The non-linearity in Scheme 32 does not allow to rely on well-known and effective solvers for the resulting linear systems.Secondly, all equations are coupled which leads to a high-dimensional mixed problem formulation.Therefore, by fully linearizing and decoupling all equations we make use of an iterative fixed-point algorithm which was presented in [45, ch. 5] and follows the ideas of [43, Algorithm A 1 ].This approach can be described by the following scheme.Table 1: Parameter choices for the experiments Accordingly the fourth and last term in the director equation (32b) of our scheme are also neglected.Both schemes can be solved by the fixed-point iteration described above.The tolerance of the iterative fixpoint solver was set to θ = 10 −6 .The python implementation of the fixed-point solver relies on the API of the finite element package FEniCS (cf.[1,39]).The implementation can be found in [46].
If not mentioned otherwise, we employ homogeneous Dirichlet boundary (no-slip) conditions for the velocity and non-homogeneous Dirichlet boundary conditions for the director as in (1).The parameter choices for the following experiments can be found in table 1.Our choice of parameters is rather exemplary and analogous to the experiments considered in [4,5,37,9].In particular regarding the choice of µ, other choices are most definitely conceivable.Note that the dissipative character of the system is fulfilled, i.e., µ 4 > 0, µ 5 + µ 6 − λ 2 > 0, µ 1 + λ 2 > 0. Further, ν = µ 4 is fulfilled which shall allow a comparison of the models.Since the convergence of our fixed-point solver relies on a contraction property, a smaller choice of our temporal discretization parameter k might sometimes still lead to faster computation time, if this leads to less total iterations of the fixed-point solver itself.

A smooth example in two spatial dimensions
We consider the domain Ω = (−1, 1) 2 equipped with the initial conditions d 0 (x) = sin (2π (cos(x 1 ) − sin(x 2 ))) cos (2π (cos(x 1 ) − sin(x 2 ))) for x ∈ Ω, Although the above analysis is only done for three spatial dimensions, the results can be transferred to the two-dimensional case as well.In this case, the cross product has to be understood in the sense x .In this setting, at least for the simplified model, a smooth solution exists except for finitely many points in time (cf.[32,Thm. 1.3]).The results for our simplified model (see (52), Fig. 2) deviate from the results of previous authors (cf.[5]) due to our Dirichlet boundary conditions.The director exhibits a long-term self-alignment which causes a fluid flow.In opposite to e.g.[5], this alignment of the director is not uniform.The same holds for our full Ericksen-Leslie model (see Fig. 3) except for the fact that induced velocity field exhibits less symmetry due to the anisotropic dynamics of our Leslie-stress tensor.
It can be observed that for both models the defects slowly annihilate by rotation of the director (see Fig. 4 and Fig. 5).This causes rough swirls to develop in the velocity field.The evolution of the director qualitatively agrees with the results in [5,4,37,9] 1 .For the full Ericksen-Leslie model, we observe again anisotropic dynamics transferred to the velocity field as well.

Annihilation of two defects in a rotating flow
The setting of the experiment equals the setting of the one before except for the choice of the initial condition and boundary condition of the velocity.Instead we choose v 0 = 10(−x 2 , x 1 , 0) T for all x ∈ Ω, v = 10(−x 2 , x 1 , 0) T for all x ∈ ∂ Ω.
Note that the inhomogeneous Dirichlet boundary conditions for the velocity are a deviation from our so far considered setting such that we can only expect the energy law (37) to hold with a modified right-hand side accounting for the boundary conditions (see Fig. 6).In Fig. 6 we can observe that in the simplified model the defects are swirled anticlockwise around the center by the velocity field before annihilation.This is in qualitative agreement with the results obtained by [4,37] and also holds for the more general model (see Fig. 7).two dimensions as done in [5] because the defects would otherwise be conserved.

∂
t d d d + (v v v • ∇)d d d − (∇v v v) skw d d d + (I −d d d ⊗d d d) (λ (∇v v v) sym d d d − ∆d d d) = 0, (1c) |d d d| = 1, (1d)where we employ the initial conditionsv v v(0) = v v v 0 , d d d(0) = d d d 0 with |d d d 0 | = 1 in Ω and boundary conditions v v v = 0, d d d = d d d Γ with |d d d Γ | = 1 on ∂ Ω.The Leslie stress tensor T L is defined byT T T L := T T T D + λ [d d d ⊗ [d d d] T x [d d d] x • ∆d d d] sym + [d d d ⊗ ∆d d d] skw ,where T D collects the dissipative terms of the Leslie stress tensor, i.e.T T T D :=(µ 1 + λ 2 )(d d d • (∇v v v) sym d d d)(d d d ⊗d d d) + µ 4 (∇v v v) sym + (µ 5 + µ 6 − λ 2 )(d d d ⊗ (∇v v v) sym d d d) sym and E ≥ E (v v v,d d d) on [0, T ] as well as |d d d| = 1 a.e. in Ω × (0, T ).The term d d d × ∆d d d has to be understood as the weak divergence of d d d × ∇d d d.The solution fulfills the energy-variational inequality The initial values are fulfilled in a weak sense v v v(0) = v v v 0 , d d d(0) = d d d 0 and d d d fulfills the inhomogeneuous Dirichlet boundary conditions in the sense of the trace tr(d d d(t)) = tr(d d d Γ ) for a.e.t ∈ (0, T ).Additionally, it holds

Lemma 2 . 6 .
Let d d d ∈ L ∞ (Ω) ∩ H 1 (Ω) with d d d × ∆d d d ∈ L 2 (Ω) and |d d d(x,t)| = 1 a.e. in Ω × (0, T ) such that for the trace | tr(d d d)| = 1 a.e. on ∂ Ω.Then we can Identify the terms −d d d × (d d d × ∆d d d) = [d d d] T x [d d d] x ∆d d d = (I −d d d ⊗d d d)∆d d d = ∆d d d + |∇d d d| 2 d d d in L 2 (Ω), which implies Ω |d d d × ∆d d d| 2 dx = Ω |∆d d d + |∇d d d| 2 d d d| 2 dx .

∂
t d d d •c c c − (∇v v v) skw d d d •c c c + (d d d × (∇d d dv v v + λ (∇v v v) sym d d d − ∆d d d)) • (d d d ×c c c) dx dt = 0 for all c c c ∈ L 2 (0, T ; L 3 (Ω)) by density arguments.The regularity of d d d suffices to infer 2d d d •(∇d d dv v v) = (v v v• ∇)|d d d| 2 such that (d d d × ∇d d dv v v) • (d d d ×c c c) = (v v v • ∇)d d d •c c c. Since the relation holds for all c c c ∈ L 2 (0, T ; L 3 (Ω))

Ω h k µ 1 µ 4 ν µ 5 , µ 6
of the matrix I − d d d ⊗ d d d (see (6)), which in three dimensions fulfills I − d d d ⊗ d d d = [d d d] T x [d d d]

Figure 2 :
Figure 2: Experiment 4.1, simplified model (52): Evolution of the director at time t = 0, 2.0 (from left to right), evolution of the energy (bottom left), velocity field at times t = 0.2 (bottom right).The colour marks the magnitude of the vectors.

Figure 3 :
Figure 3: Experiment 4.1, model (32): Evolution of the director at time t = 0, 2.0 (from left to right), evolution of the energy (bottom left), velocity field at times t = 0.2 (bottom right).The colour marks the magnitude of the vectors.

Figure 4 :
Figure 4: Experiment 4.2, simplified model (52): Evolution of the director in the plane x 3 = 0 at time t = 0, 0.03, 0.05, 0.1 (from left to right, from top to bottom), evolution of the energy (bottom left), velocity field in the plane x 3 = 0 at time t = 0.05 (bottom right).

Figure 5 :
Figure 5: Experiment 4.2, model (32): Evolution of the director in the plane x 3 = 0 at time t = 0, 0.03, 0.05, 0.1 (from left to right, from top to bottom), evolution of the energy (bottom left), velocity field in the plane x 3 = 0 at time t = 0.05 (bottom right).The colour marks the (absolute) value of the z-component.

Figure 6 :
Figure 6: Annihilation in a rotating flow, simplified model (52): Evolution of the director in the plane x 3 = 0 at time t = 0.03, 0.05, 0.15, 0.5, evolution of the energy, velocity field in the plane x 3 = 0 at time t = 0.15 (from top left to bottom right).The colour marks the (absolute) value of the z-component.

Figure 7 :
Figure 7: Annihilation in a rotating flow, model (32): Evolution of the director in the plane x 3 = 0 at time t = 0.03, 0.05, 0.15, 0.5, evolution of the energy, velocity field in the plane x 3 = 0 at time t = 0.15 (from top left to bottom right).The colour marks the (absolute) value of the z-component.