The Ericksen Model of Liquid Crystals with Colloidal and Electric Effects
Abstract
We present a robust discretization of the Ericksen model of liquid crystals with variable degree of orientation coupled with colloidal effects and electric fields. The total energy consists of the Ericksen energy, a weak anchoring (or penalized Dirichlet) energy to model colloids, and an electrical energy for a given electric field. We describe our special discretization of the total energy along with a method to compute minimizers via a discrete quasigradient flow algorithm which has a strictly monotone energy decreasing property. Numerical experiments are given in two and three dimensions to illustrate that the method is able to capture nontrivial defect patterns, such as the Saturn ring defect. We conclude with a rigorous proof of the convergence of our discrete energy to the continuous energy.
keywords:
liquid crystals, finite element method, gammaconvergence, gradient flow, line defect, plane defect, Saturn ring defectMsc:
65N30, 49M25, 35J701 Introduction
This paper presents a method for solving the Ericksen model of liquid crystals (1); (2), with additional effects due to colloidal domains and electric fields. Liquid crystals are a workhorse technology enabling electronic displays (3); (4); (5), for instance. Moreover, they have a host of potential applications in material science (6); (7); (8); (9); (10); (11); (12); (13); (14); (15); (16); (17); (18); (19); (20); (21). One avenue is to use external fields (e.g. electric fields) and colloidal dispersions to build new materials through directed selfassembly (22); (7); (12); (23); (24); (25); (26); (27); (28); (29); (16); (19); (30); (31).
A significant amount of mathematical analysis has been done on liquid crystals (32); (33); (34); (35); (36); (37); (38); (39); (40); (41); (42); (43); (44). Moreover, a host of numerical methods have been developed for statics and dynamics (45); (46); (47); (48); (49); (50); (51). In particular, the methods in (52); (43); (53); (54); (55) are for harmonic mappings and liquid crystals with fixed degree of orientation, i.e. a unit vector field (called the director field) represents the orientation of liquid crystal molecules. See (56); (57); (58); (59); (60) for methods that couple liquid crystals to Stokes flow. We also refer to the survey paper (45) for more numerical methods.
The method we present (61); (62) is for the oneconstant model of liquid crystals with variable degree of orientation (1); (2); (32) (Ericksen’s model). The state of the liquid crystal is described by a director field and a scalar function , the socalled degreeoforientation, which minimize the energy
(1) 
Hereafter, is a material constant, is a bounded Lipschitz domain in with , and is a double well potential (defined below).
Minimizers of the Ericksen model may exhibit nontrivial defects (depending on boundary conditions) (63); (10); (64); (39); (38); (65). The presence of in (1) leads to an EulerLagrange equation for that is degenerate. This allows for line and plane defects (singularities of ) in dimension when vanishes; these types of defects are important for applications, especially defects that lie on three dimensional space curves (7); (66). Regularity properties of minimizers, and the size of defects, were studied in (39). This leads to the study of dynamics (33) and corresponding numerics (46), which are relevant to our paper. But in both cases they regularize the model to avoid the degeneracy associated with the order parameter vanishing.
Our finite element method (FEM) does not require any regularization. We discussed the mathematical foundation of our method in (61); (62): we proved stability and convergence properties via convergence (67) (as the mesh size goes to zero) and developed a quasigradient flow method to solve the discrete problem. Our discretization of the energy (1), defined in (16), requires that the mesh be weakly acute (or the stronger condition of having nonobtuse angles). This discretization preserves the underlying structure and robustly handles the unit length constraint on and the degeneracy present when vanishes. Our previous paper (61) showed a variety of simulations of minimizers with interesting defect structures.
The present paper demonstrates the ability of the Ericksen model, and of our method, to capture defect structures induced by colloidal inclusions (i.e. holes in the domain) and effects due to electric fields. We are able to recreate the famous Saturn ring defect (68), which occurs around colloidal particles in different situations, by using both a conforming (nonobtuse mesh) and a nonconforming cube mesh with an immersed boundary approach to model the colloid. In addition, we include electric field effects by incorporating an electric energy term into the total energy (1), and demonstrate the classic Freedericksz transition (32); (2); (41); (69); (70). We also investigate the coupling of colloidal and electric effects.
The paper is organized as follows. In Section 2, we recall the Ericksen model for liquid crystals with variable degree of orientation, and describe our discretization of the continuous energy. Section 3 recalls properties of the discretization, and our initial minimization scheme. Section 4 describes the details for properly implementing our method. Section 5.1 illustrates our method in the presence of a colloidal inclusion with a conforming nonobtuse mesh. Section 5.2 shows an alternative way to model colloids by an immersed boundary approach (along with supporting simulations). In Section 6, we show how to include electric field effects in the model and describe a modified minimization procedure to compute minimizers. Section 7 presents the monotone energy decreasing property of the quasigradient flow algorithm to compute discrete minimizers. Section 8 provides a summary of the convergence for our discrete energy. We close in Section 9 with some discussion.
2 Ericksen’s model
Let the director field be a vectorvalued function with unit length (see Figure 1 for a description of the meaning of ). The degreeoforientation is a real valued function (see Figure 2 for a description of the meaning of ). The variable , by itself, cannot properly describe a “loss of order” in the liquid crystal material because it has unit length. The variable provides a way to characterize the local order (see Figure 2).
2.1 Ericksen’s one constant energy
The equilibrium state of the liquid crystal material is described by the pair minimizing a bulkenergy functional (1) which we split as
(2) 
where . The double well potential is a function defined on that satisfies (1); (35); (39)

,

for some ,

.
It was shown in (35); (39) that introducing an auxiliary variable allows the energy to be rewritten as
(3) 
which follows from the orthogonal decomposition (and is due to the constraint ). Hence, (35); (39) define the admissible class of solutions (minimizers) as
(4) 
We may also enforce boundary conditions on , possibly on different parts of the boundary. Let be open subsets of where we set Dirichlet boundary conditions for . Then we have the following restricted admissible class
(5) 
for some given functions that satisfy the following in a neighborhood of : and , for some . Note that if we further assume
(6) 
then the function is in a neighborhood of and satisfies on .
When the degree of orientation is a nonzero constant, the energy in (2) effectively reduces to the OseenFrank energy . The purpose of the degree of orientation is to relax the energy of defects. In fact, discontinuities in (i.e. defects) may still occur in the singular set
(7) 
with finite energy: . The existence of minimizers in the admissible class, subject to Dirichlet boundary conditions, was shown in (35); (39). Minimizers with defects are constructed explicitly in (32) or discovered numerically in (61).
The parameter in (2) plays a major role in the occurrence of defects. Assuming the boundary condition for is a positive constant well away from zero, if is large, then dominates the energy and stays close to a positive constant within the domain . Thus, defects are less likely to occur. If is small (say ), then dominates the energy, and may vanish in regions of and induce a defect. This was confirmed by our numerical experiments in (61); (62). The physically relevant case is the more difficult case with regard to proving convergence (see (61)) because the energy is no longer convex.
2.2 Discretization of the energy
Let be a conforming simplicial triangulation of . The set of nodes (vertices) of is denoted and has cardinality . We demand that be weakly acute, namely
(8) 
where is the standard “hat” basis function associated with node . We indicate with the patch of a node (i.e. the “star” of elements in that contain the vertex ). Of course, (8) imposes a severe geometric restriction on (71); (72) (especially in three dimensions). We recall the following characterization of (8) for .
Lemma 1 (weak acuteness in two dimensions).
For any pair of triangles , in in two space dimensions that share a common edge , let be the angle in opposite to (for ). Then (8) holds if and only if for every edge .
Generalizations of Lemma 1 to three dimensions, involving interior dihedral angles of tetrahedra, can be found in (73); (74). Note that a nonobtuse mesh (one where all interior angles are bounded by ) is automatically weaklyacute.
The method uses the following finite element spaces:
(9) 
where imposes the unit length constraint at the vertices of the mesh.
Let denote the piecewise linear Lagrange interpolation operator on mesh with values in either or . We have the following discrete version of the admissible class:
(10) 
Next, we let and be the discrete Dirichlet data, and introduce the discrete spaces that include (Dirichlet) boundary conditions
as well as the discrete admissible class with boundary conditions:
(11) 
In view of (6), we can also impose the Dirichlet condition on .
Our discrete version of is “derived” by invoking basic properties of the stiffness matrix entries . First note for all , and for piecewise linear we have
Thus, using and the symmetry , we get
(12)  
where we define
(13) 
Therefore, we define the discrete energy to be
(14) 
where the second term is a first order approximation of , which is novel in the finite element literature (61). The double well energy is discretized in the usual way:
(15) 
The specific form of (14) lies on the fact that it makes the nodal values of and readily accessible for analysis. The identity in (3) is obtained (at the continuous level) by taking advantage of the unit length constraint . However, at the discrete level, we only impose the unit length constraint at the nodes of the mesh, and we cannot hope for much more because is a piecewise polynomial. Hence, in order to obtain a similar identity to (3) (see Lemma 2 below), we need access to nodal values. In (61), we show that the discrete energy (14) allows us to handle the degenerate coefficient without regularization.
The discrete formulation is as follows. Find such that the following energy is minimized:
(16) 
3 Review of the method
3.1 Energy inequality
Our discrete energy (14) satisfies a discrete version of (3) (35); (39), which is a key component of our analysis in (61). To see this, we introduce and two discrete versions of the vector field
(17) 
Note that both , are in . We now state a discrete version of (3).
Lemma 2 (discrete energy inequality).
3.2 Minimization scheme
We summarize the discrete quasigradient flow scheme in (61) which we use to compute discrete minimizers.
Boundary conditions
In the continuous setting, Dirichlet boundary conditions are enforced in the space. Let
(22) 
where is either , and , . We assume and (6) to be valid on . The trace is thus well defined on .
The superscript will stand for an iteration counter. Therefore, and indicate iterates satisfying Dirichlet boundary conditions where and . We will further simplify the notation in some places upon writing:
First order variation
We start with the energy . Due to the unit length constraint at the nodes in (see (9)), we introduce the space of discrete tangential variations:
(23) 
Next, the first order variation of in the direction reads
(24) 
The first order variation of in the direction consists of two terms
(25) 
We next consider the energy . In order to guarantee a monotonically energy decreasing scheme, we employ the convex splitting technique in (75); (76); (77), i.e. we split the double well potential into a convex and concave part. Let and be both convex for all so that , and set
(26) 
which yields the inequality
(27) 
for any and in (61). Note that
Discrete quasigradient flow algorithm
Our scheme for minimizing , defined in (16), is an alternating direction method, which minimizes with respect to and evolves separately in the steepest descent direction during each iteration. Therefore, this algorithm is not a standard gradient flow but rather a quasigradient flow.
Algorithm: Given in , iterate Steps (a)(c) for .
Step (a): Minimization. Find such that minimizes the energy for all in , i.e. satisfies
Step (b): Projection. Normalize at all nodes .
Step (c): Gradient flow. Using , find in such that
In the numerical experiments, we impose Dirichlet boundary conditions for both and on subsets of the boundary. Note that the scheme has no restriction on the time step thanks to the implicit Euler method in Step (c).
The quasigradient flow scheme has a monotone energy decreasing property, provided the mesh is weakly acute (8) (71); (72).
Theorem 3 (monotonicity (61)).
4 Implementation
We implemented our method using the MATLAB/C++ finite element toolbox FELICITY (78). In this section, we give details on forming the ensuing discrete systems, and how to solve part (a) of the quasigradient flow algorithm in 3D using the tangent space. For all 3D simulations, we used the algebraic multigrid solver (AGMG) (79); (80); (81); (82) to solve the linear systems in parts (a) and (c) of the quasigradient flow algorithm. In 2D, we simply used the “backslash” command in MATLAB.
4.1 Finite element matrices
Implementing the algorithm requires construction of the discrete energy, as well as its variational derivative. This requires the symmetric mass and stiffness finite element matrices: , , where
(28) 
and is the set of basis functions of the space .
4.2 Finite element functions and coefficient vectors
The function is represented by a linear combination of . If the dimension of is , then the vector field has components, where each component is written as a linear combination of . The nodal values of and , at node , are denoted by and .
The corresponding coefficient vectors (arrays) are denoted with nonitalicized capital letters, i.e. , , such that
(29) 
where and
(30) 
for , where are the canonical basis vectors of . In other words, we store the coefficients of so that the components are first, followed by the components, and so on. Therefore,
4.3 Discrete variations
Let us write from (25) in a different form:
(31) 
where is the coefficient vector corresponding to , and is the symmetric matrix defined by
(32) 
Lemma 4.
Let be an arbitrary matrix, and , be arbitrary column vectors. Then
(33) 
where is a column vector and is a diagonal matrix formed from .
Next, we write out from (24) in a different form:
(35) 
where we defined
(36) 
with the coefficient vector corresponding to .
Let us now focus on :
(37) 
where is the symmetric matrix defined by
(38) 
Using symmetry gives
(39) 
where denotes the first components of , etc., and
(40) 
4.4 Discrete quasigradient flow
Given , we have the corresponding coefficient vectors . We now rewrite the Algorithm in Section 3.2.3 in terms of the matrices and vectors introduced earlier.
Step (a): By (42), we solve the following linear system in the tangent space (see Section 4.5) to obtain :
(43) 
where in is the coefficient vector corresponding to in . Note that (43) must be modified to enforce Dirichlet boundary conditions (if necessary).
Remark 5 (solving a degenerate system).
The system matrix in (43) is symmetric positive semidefinite, which is easily verified from the properties of . Moreover, it is positive definite if everywhere. Hence, the system can be solved by any method for symmetric positive definite matrices.