The Ericksen Model of Liquid Crystals with Colloidal and Electric Effects

The Ericksen Model of Liquid Crystals with Colloidal and Electric Effects


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 quasi-gradient 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 non-trivial 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.

liquid crystals, finite element method, gamma-convergence, gradient flow, line defect, plane defect, Saturn ring defect
65N30, 49M25, 35J70

1 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 work-horse 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 self-assembly (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 one-constant 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 so-called degree-of-orientation, which minimize the energy


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 non-trivial defects (depending on boundary conditions) (63); (10); (64); (39); (38); (65). The presence of in (1) leads to an Euler-Lagrange 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 quasi-gradient 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 non-obtuse 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 (non-obtuse mesh) and a non-conforming 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 non-obtuse 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 quasi-gradient 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 vector-valued function with unit length (see Figure 1 for a description of the meaning of ). The degree-of-orientation 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).

Figure 1: Macroscopic order parameter: the director variable . Nematic liquid crystal molecules have an elongated rod-like shape (see elongated ellipsoids), which gives the material its anisotropic nature. The value of (a unit vector), at the point , represents a probabilistic average over a local ensemble of liquid crystal molecules “near” (32).
Figure 2: Macroscopic order parameter: the degree-of-orientation variable . It is a probabilistic average over the angle between and an individual liquid crystal molecule; the average is taken over a local ensemble (32). The case represents the state of perfect alignment in which all molecules (in the ensemble) are parallel to . Likewise, represents the state of microscopic order in which all molecules (in the ensemble) are perpendicular to . When , the molecules do not lie along any preferred direction which represents the state of an isotropic, uniformly random, distribution of molecules (in the ensemble). The state is called a defect in the liquid crystal material.

2.1 Ericksen’s one constant energy

The equilibrium state of the liquid crystal material is described by the pair minimizing a bulk-energy functional (1) which we split as


where . The double well potential is a function defined on that satisfies (1); (35); (39)

  1. ,

  2. for some ,

  3. .

It was shown in (35); (39) that introducing an auxiliary variable allows the energy to be rewritten as


which follows from the orthogonal decomposition (and is due to the constraint ). Hence, (35); (39) define the admissible class of solutions (minimizers) as


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


for some given functions that satisfy the following in a neighborhood of : and , for some . Note that if we further assume


then the function is in a neighborhood of and satisfies on .

When the degree of orientation is a non-zero constant, the energy in (2) effectively reduces to the Oseen-Frank 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


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


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 non-obtuse mesh (one where all interior angles are bounded by ) is automatically weakly-acute.

The method uses the following finite element spaces:


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:


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:


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


where we define


Therefore, we define the discrete energy to be


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:


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:


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


Note that both , are in . We now state a discrete version of (3).

Lemma 2 (discrete energy inequality).

Let the mesh satisfy (8). If , then, for any , the discrete energy (14) satisfies


as well as


In fact, the following inequalities are valid (61)




Note that because for . We refer to (61) for further details.

3.2 Minimization scheme

We summarize the discrete quasi-gradient 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


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:


Next, the first order variation of in the direction reads


The first order variation of in the direction consists of two terms


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


which yields the inequality


for any and in (61). Note that

Discrete quasi-gradient 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 quasi-gradient 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 quasi-gradient flow scheme has a monotone energy decreasing property, provided the mesh is weakly acute (8) (71); (72).

Theorem 3 (monotonicity (61)).

Let satisfy (8). The iterate of the Algorithm (discrete quasi-gradient flow) of Section (3.2.3) exists and satisfies

Equality holds if and only if (equilibrium state).

We extend this result in Theorem 8 to include additional energy terms; see (101).

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 quasi-gradient flow algorithm in 3-D using the tangent space. For all 3-D simulations, we used the algebraic multi-grid solver (AGMG) (79); (80); (81); (82) to solve the linear systems in parts (a) and (c) of the quasi-gradient flow algorithm. In 2-D, 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


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 non-italicized capital letters, i.e. , , such that


where and


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:


where is the coefficient vector corresponding to , and is the symmetric matrix defined by

Lemma 4.

Let be an arbitrary matrix, and , be arbitrary column vectors. Then


where is a column vector and is a diagonal matrix formed from .

Continuing with (31), and using (33), we get


Next, we write out from (24) in a different form:


where we defined


with the coefficient vector corresponding to .

Let us now focus on :


where is the symmetric matrix defined by


Using symmetry gives


where denotes the first components of , etc., and


Applying the same argument to , we get


Therefore, recalling (35), we obtain


where is block diagonal (with identical blocks), where each block equals defined in (40); similarly, is block diagonal (with identical blocks), where each block equals defined in (38).

4.4 Discrete quasi-gradient 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 :


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 semi-definite, 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.

When at a sufficient number of nodes, the matrix will be singular. In this case, one could use a conjugate gradient method (83); note that the right-hand-side of (43) is guaranteed to be in the column space of the system matrix.

In 3-D, we solve (43) using AGMG (79); (80);