Full paper - Institute for Computational Mathematics - Kent State ...
Structures and Structural Phase Transitions in Confined Liquid Crystal Systems Eugene C. Gartland, Jr.∗ Department of Mathematics and Computer Science Kent State University P.O. Box 5190 Kent, Ohio 44242-0001 USA
Abstract Discussed are two numerical packages for computing equilibrium structures of liquid crystals. The first uses Frank elasticity in a rectangular geometry in two spatial dimensions. The second uses a tensor order parameter and Landau-deGennes free energy in the geometry of a finite circular cylinder. The Frank code is being used to study the phenomenon of the “stripe phase” in thin nematic films in the bendFr´eedericksz geometry. The cylinder code is being used to investigate configurations in liquid crystals confined in microscopic cylindrical pores subject to homeotropic anchoring conditions. Preliminary results of these projects are presented.
1
Introduction
The Liquid Crystal Institute at Kent State University is the locus of the Science and Technology Center for “Advanced Liquid Crystalline Optical Materials” (ALCOM), which is funded by the National Science Foundation and the State of Ohio. This Center is principally concerned with the development (analysis, characterization, and synthesis) of new liquid crystalline materials for applications in various kinds of display devices. The effort is broad and inter-disciplinary; it involves researchers in physics, chemistry, mathematics, ∗
This research was partially supported by the National Science Foundation Science and Technology Center on “Advanced Liquid Crystalline Optical Materials” (ALCOM) under grant DMR 89-20147 and by NSF grant DMS-9310733.
1
and computer science, working on problems from basic theory to practical application, in programs involving theory, experiment, and computation. In recent years, we have been involved in various activities of numerical computation in support of this effort. Most of this has been concerned with exploring the static-equilibrium structures of the molecular orientational properties of liquid crystal materials in confinement. Of particular interest are the geometries of droplets, cylindrical pores, and films. We give here a preliminary report on two projects in this area, both of which are currently in progress. The first project uses Frank elasticity to study the development of so-called “stripe” or “ripple” phases in nematic thin films in the presence of a magnetic field at temperatures just above the nematic to smectic-A transition point. The divergence of the bend elastic constant in such a parameter range favors the development of a translational-symmetry-breaking configuration with periodicity in the plane of the film, perpendicular to the applied field. These textures have been observed experimentally, and there exists some theory, upon which our computational results build. The second project uses the tensor order parameter and Landau-deGennes free energy to study co-existing structural phases of materials confined in microscopic cylindrical pores subject to homeotropic anchoring conditions. In addition to the familiar “escaped” radial director field, two other configurations are found. One preserves the rotational symmetry of the problem; while the other one breaks it. Both involve line defects: one of unit strength running down the axis of the cylinder, the other possessing a pair of half-strength disclinations. The results are consistent with experimental evidence and other numerical simulations.
2
Frank minimization
We have developed a general numerical code to compute equilibrium director configurations in a rectangular geometry in two spatial dimensions by minimizing a discretized Frank free energy. We have adapted the code to numerically study “stripe” phases in thin nematic films. We report on aspects of the code and on the results obtained so far on the study of the stripe phases.
2.1
Model and code
The free energy utilized is b = F (n)
1 Z dZ b b dxdy , f (n) 2 c a
2
with the density given by the standard Oseen-Z¨ ocher-Frank elastic model with magnetic field: b := K1 (∇ · n) b 2 + K2 (n b · (∇ × n)) b 2 f (n) b × (∇ × n)| b 2 − χa (H · n) b 2 . + K3 |n
(1)
At present, boundary conditions are allowed to be Dirichlet or free (on any of the four sides) or periodic in the x direction and/or periodic in the y direction. We are in the process of adding a general surface free-energy density and including the K24 “saddle-splay” elastic constant term (as a surface term). In the future, we plan to modify the code to include a scalar order parameter, s(x, y), to allow for a locally variable degree of order (along the lines of the treatment by Ericksen in [10]). The problem is discretized by introducing a rectangular grid, which is not required to be uniform in either direction. The discrete director is represented at each mesh point by a unit vector. The free-energy functional is approximated by a composite mid-point quadrature rule, using cell-centered finite differences and averages. Local minimizers are found by using a projectedconjugate-gradient line search procedure similar to that employed in [7]. The minimization algorithm proceeds in the following way. At each step of the iteration, the gradient of the discrete free energy is computed and projected onto the tangent space of the constraint set at that configuration—this process is local in nature, involving projecting threevectors onto tangent planes to unit spheres. A search direction is built up from these using the Fletcher-Reeves conjugate gradient recurrence relations (see, for example, [11]). A one-dimensional search is then conducted in this direction but always remaining on the constraint set. The approach differs from other published approaches (e.g., [3, 6, 7, 15]) in this regard. Those schemes all perform a true “line” search, allowing the discrete director field to move off of the constraint set, violating its pointwise unit-length constraint, and then re-normalizing all the vectors once the line search is terminated. In all of the situations in which these other codes have been employed, this re-normalization step was provably free-energy reducing [3]. For our application below (and for others we envision), this is not the case. The structures and configurations we seek are also rather fragile, and so we feel that this greater fidelity to the constraints should enhance the robustness of the code. Some additional features have been added to the code to improve its utility and performance. These include a primitive multi-level capability, whereby successive refinements of an initial mesh can be generated by cell-edge bisections. Computed approximate solution fields are interpolated to finer meshes via a bi-cubic spline interpolation routine. The discrete free-energy functional 3
admits an asymptotic error expansion in even powers of the cell dimensions, and this makes extrapolation a very effective procedure to achieve highly accurate free-energy approximations (which are needed for producing good phase diagrams in the ill-conditioned situations, such as near Fr´eedericksz transitions). A simple, single-parameter continuation procedure is also implemented. It uses an application-specific routine to update the problem parameters in terms of changes in a single path parameter. This allows one to follow an arbitrary trajectory in the parameter space. Initial guesses for the unknown discrete director field are generated from a secant-line approximation through the two most recently computed solution fields. The approach is rather simplistic, but it is helpful in collecting data for phase diagrams and the like. The code is written in Fortran double precision and developed for a workstation environment. In its current incarnation (for a Hewlett-Packard workstation), it utilizes HP Matrix-Vector Library routines as well as routines from the NAG Software Library. Matlab is used for the graphics and visualizations. The code was benchmarked against a uniform Fr´eedericksz distortion computed independently by iteratively solving the implicit analytical solution in terms of elliptic integrals. For a realistic computation of 64 × 64 mesh cells, there are on the order of 4,000 nodes (and constraints) and 12,000 unknowns. With no preconditioning, the code can take 3–5 minutes to converge to the level of the discretization error from a moderately good initial guess. The accuracy of the approximate free energy is of the order of four decimal digits (in absolute error terms); the accuracy of the discrete director field is of the order of two digits. In continuation mode, with high-accuracy initial guesses available, convergence is much faster. We have built in a very simple preconditioner that utilizes rank-one perturbations of the identity. This can be useful, for instance, when collecting data near a point where one has some analytical (or intuitive) information about a mode that is losing stability (e.g., near a Fr´eedericksz transition). This can help dramatically in some cases, but it is not general, and it is of limited applicability. There are several positive features that we can point to in such a code. It is relatively easy to construct: it just requires a discrete free energy and its gradient. No second-derivative (Hessian) information is used, and so no largesparse matrix linear algebra and data structures are required. When dealing with periodic (or doubly periodic) conditions, no “phase condition” is required: the gradients (and projected gradients and conjugate projected gradients) are all guaranteed to be “orthogonal to the periodicity.” Negative features include the relatively low order of the discretization
4
method. Though helped by extrapolation, it makes it very difficult to compute fragile structures on coarse meshes. The averaging used in the cell-centered difference formulas also admits some troublesome non-physical solutions—these typically have a high-frequency micro-structure at the level of the mesh. Finally we mention that the lack of second-derivative information makes it impossible to compute (or path follow) unstable stationary points; this limits the usefulness of the code for full bifurcation studies.
2.2
Application: stripes in nematic thin films
An outstanding problem in the physics of nematic films is the stripe phase that appears in the bend-Fr´eedericksz geometry at temperatures slightly above that of the nematic to smectic-A phase transition in a sufficiently strong magnetic field. Here we have a material confined to a cell 0 < x < d of assumed infinite extent in the y and z directions, satisfying homeotropic anchoring conditions b = x) b at the plate boundaries x = 0 and x = d. A magnetic field is applied (n in the plane of the film, in the z direction: H = H zb . The problem geometry is depicted in Figure 1. z
y H
l periodic
x d
Figure 1: Geometry and boundary conditions for the bend-Fr´eedericksz transition and competing stripe phase solution. In [2], a linear stability analysis was used with a two-parameter ansatz to show that at magnetic fields slightly above the critical Fr´eedericksz threshold, a phase with periodicity in the y direction has a lower free energy than the uniform Fr´eedericksz distortion. We have used the code discussed in §2.1 to numerically study this problem. We model the system as follows. Take as free energy 1 Z lZ d b l) = b dxdy , F (n, f (n) 2l 0 0
5
b given as before by (1). Thus F is a free-energy density per unit area, with f (n) the free-energy density per unit length in z of one period in y (normalized by b and l, subject the length of the period). We seek to minimize F over both n b =x b on x = 0 and x = d, and periodic to homeotropic anchoring conditions n conditions on y = 0 and y = l. We then compare the free-energy density of this single periodic cell with that of the Fr´eedericksz configuration, which is uniform in the y direction. We non-dimensionalize as in [2]: b l) = F (n,
Z
1/2Z 1
0
b l) dxdy , f (n,
0
where
2
b l) := M ∇ · n b f (n,
2
b · ∇×n b + MN n
2
b × ∇×n b − h2 π 2 n2z . + n
Here x := x/d, y := y/2l, l := 2l/d (twice the period, in units of the cell thickness, or the ratio of the y scaling to the x scaling), and !
∇ :=
∂ 1 ∂ , . ∂x l ∂y
F is in units of K3 /d. The parameters M, N, and h are defined by K1 M := , K3
K2 N := , K1
Hd h := π
s
χa , K3
this last representing the magnetic field in units of the critical field (i.e., h = 1.0 corresponds to the critical Fr´eedericksz threshold). We have fixed the computational domain of the problem on (0, 1) × (0, 1/2) because previous experiments and analyses have indicated that the period of the stripe phases is in the range of .3 to .7 times the cell thickness. We are interested in the range where M is small: the bend elastic constant (K3 ) diverges relative to the splay constant (K1 ) as TN A is approached from the right. The twist constant (K2 ) also diverges, but at a much more moderate rate, and so N is treated as O(1). b l) (on a fixed computational domain) Our problem then is to minimize F (n, b and l, subject to the appropriate boundary and pewith respect to both n riodic conditions (and the pointwise unit-length constraint). The free-energy functional F has the structure 1 1 b l) = F 0 (n) b + F 1 (n) b b 2. F (n, + F 2 (n) l l 6
We perform the numerical minimization by a “decoupling iteration,” first minb (using the projected-conjugate-gradient line search imizing with respect to n discussed above), and then minimizing with respect to l via the exact relation lν+1 = −2
h b hν+1 F2 n h
b hν+1 F1 n
,
properly safeguarded—repeating until convergence. This approach can be thought of as a type of fractional-step method. It can also be viewed as a form of nonlinear block Gauss-Seidel. It has the advantage that it is relatively easy to implement on top of a basic existing Frank-minimization code, such as we have already developed. The main results we have obtained thus far include the verification of the first-order nature of the phase transition (from the uniform Fr´eedericksz distortion to the periodic stripe phase); there is a definite overlap in the regions of meta-stability, which can be relatively large in certain ranges of parameters. We have also developed a refinement of the phase diagram (Figure 2) from [2]. It retains the qualitative features of the one in that study but indicates the global stability of the stripe phase at higher elastic-constant ratios M. Stripe Regime Thresholds 1.1 1.09
original new
Scaled Applied Field ( h )
1.08 1.07 N = .75
1.06
N = .50
1.05 1.04
N = .25
1.03 1.02 1.01 1 0
0.02
0.04 0.06 0.08 0.1 0.12 Elastic Constants Ratio ( M = K1 / K3 )
0.14
0.16
Figure 2: Stripe threshold field, original vs new, N = .25, .50, .75. For the indicated values of N, the stripe phase has lower free energy than the uniform Fr´eedericksz distortion for values of M to the left of the threshold (crossover) lines. The most significant result we have established so far is the fact that these stripe phases are actually stable below the critical Fr´eedericksz threshold. Figure 3 shows representative structures for stripe phases at three different field 7
M = 0.126, N = 0.25, h = 1.03
M = 0.14, N = 0.25, h = 1.07
M = 0.09, N = 0.25, h = 0.99
0.8 0.7 0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0
−0.1
−0.1
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.2
M = 0.14, N = 0.25, h = 1.07
0.4
0.6
0.8
1
0
0.2
0.1
0.1
0.1
0
0
0
−0.1
−0.1
−0.1
−0.2
−0.2 0.2
0.4
0.6
0.8
0
0.2
0.4
0.6
0.8
1
0.2
0.2
0.2
0.1
0.1
0.1
0
0
0
−0.1
−0.1
−0.1
−0.2 0
0.2
0.4
0.6
0.8
0
0.2
M = 0.14, N = 0.25, h = 1.07
0.4
0.6
0.8
1
1
0.5
0.15
0.2
0.25 y
0.3
0.35
0.4
0.45
0.5
0.05
0.1
0.15
0.2
0.25 y
0.3
0.35
0.4
0.45
0 0
0.5
0
0
0
−0.5
−0.5
−0.5
20
25
30
1
0.4
0.6
0.8
1
0.05
0.1
0.15
0.2
0.25 y
0.3
0.35
0.4
0.45
0.5
A
1 0.5
A
1 0.5
A
1
15 frequency
0.8
1.5
0.5
10
0.6
Discrete Fourier Cosine Coefficients
1.5
5
0.2
Discrete Fourier Cosine Coefficients
1.5
0
0
0.4
0.5
0 0
Discrete Fourier Cosine Coefficients
−1
0.2
theta
1 theta
1 theta
1.5
0.1
1
M = 0.09, N = 0.25, h = 0.99
1.5
0.05
0
M = 0.126, N = 0.25, h = 1.03
1.5
0 0
0.8
−0.2
1
0.5
0.6
−0.2
1
−0.2
0.4
M = 0.09, N = 0.25, h = 0.99
0.2
0
0.2
M = 0.126, N = 0.25, h = 1.03
0.2
−1
0
5
10
15 frequency
20
25
30
−1
0
5
10
15 frequency
20
25
30
Figure 3: Stripe phases (down-sampled from 64 × 63 discretization), N = .25, (M, h) = (.140, 1.07), (.126, 1.03), (.090, 0.99), (a) top view, (b) side view vs (c) uniform Fr´eedericksz distortion, (d) vertical angle in cell center (as a function of y, vs θm for the uniform Fr´eedericksz distortion) and (e) its discrete Fourier cosine coefficients. strengths: 7% and 3% above critical (h = 1.07, 1.03) and 1% below (h = 0.99). At high fields (where there is little “room to maneuver”), the stripe phase appears as a modest sinusoidal perturbation on top of the uniform Fr´eedericksz distortion. At lower field strengths (in particular, below the critical threshold), the stripe phase uses relatively cheap splay and twist deformations to saturate with the applied field to a high degree in the core of the cell. Negative results include the fact that we are unable to obtain any of these stripe configurations on meshes coarser than 32 × 32. Thus we have not been able to take full advantage of multi-level computing and extrapolation. This would appear to be related to the relative fragility of these structural phases and to the low order of accuracy of the discretization model. Our path8
following steps also were required to be relatively short, as the iterations would readily become attracted to spurious, non-physical solutions.
3
Landau-deGennes minimization
We have developed a numerical code to compute equilibrium configurations of the tensor order parameter in the geometry of a circular cylinder by minimizing a discretized Landau-deGennes free energy [8, 9]. We are using this code numerically to study structural phases in microscopic cylindrical pores. We report on aspects of the code and on the results we have obtained so far on the study of the phases in a cylinder.
3.1
Model and code
We model our system in terms of the tensor order parameter, Q. Thus Q(x) is a symmetric, traceless tensor field defined (at each point x ∈ Ω) by Qαβ :=
1 2
< 3lα lβ − δαβ > ,
α, β = 1, 2, 3 ,
where < · > denotes an average with respect to the orientational probability distribution function at the point, and lα are the components of a unit vector in the direction of the long axis of a molecule. The Q tensor contains information about the preferred direction(s) of orientation as well as the degree of order. This information is most readily interpreted in terms of the eigenvalues and eigenvectors of Q : (λi , ui ), i = 1, 2, 3. Ordering these real eigenvalues from right to left, they satisfy − 12 ≤ λ3 ≤ λ2 ≤ λ1 ≤ 1 . The associated eigenvectors can be taken to be orthonormal and give extremals of the mean-square projection of the axis of a molecule: n
o
2λi + 1 . 3 The method we currently use to visualize the Q tensor (see Figure 4) is to draw an ellipsoid, with axes in the directions of the eigenvectors ui , scaled relatively by the mean-square-projection factors (2λi + 1)/3 (the values of which are in [0, 1]). Special (limiting, extreme) cases are max < (u · l)2 > : u ⊥ {u1 , . . . , ui−1 } , |u| = 1 = < (ui · l)2 > =
perfectly ordered with a positive order parameter: λ3 = λ2 = − 12 ,
λ1 = 1 ,
equivalently a3 = a2 = 0, a1 = 1, in which case the ellipsoid degenerates to a line segment (in the u1 direction); 9
perfectly ordered with a negative order parameter: λ3 = − 12 ,
λ2 = λ1 = 14 ,
or a3 = 0, a2 = a1 = 1/2, in which case the ellipsoid degenerates to a disk (perpendicular to the u3 direction); isotropic: λ3 = λ2 = λ1 = 0 , or a3 = a2 = a1 = 1/3, in which case the ellipsoid is a sphere. The material is said to be “biaxial” (at a point) if λ3 < λ2 < λ1 . It is “uniaxial” if λ3 = λ2 < λ1 (positive order parameter) or λ3 < λ2 = λ1 (negative order parameter). We utilize a Landau-deGennes free energy of the following form: F (Q) =
Z Ω
{ felastic (∂Q) + fbulk (Q) + ffield (Q) } +
Z ∂Ω
fsurf (Q) ,
where the densities are given (in compact tensor notation) by L1 Qαβ,γ Qαβ,γ + 12 L2 Qαβ,β Qαγ,γ + 12 L3 Qαβ,γ Qαγ,β A tr(Q2 ) − 13 B tr(Q3 ) + 14 C tr(Q2 )2 + 15 D tr(Q2 )tr(Q3 ) + 16 E tr(Q2 )3 + 16 E 0 tr(Q3 )2 ffield (Q) := − tr(F Q)
felastic (∂Q) := fbulk (Q) :=
fsurf (Q) :=
1 2 1 2
1 2
(2)
W tr (Q − Q0 )2 .
Here L1 , L2 , . . . , are material-dependent constants; Qαβ,γ := ∂Qαβ /∂xγ ; tr(·) denotes the ordinary trace of a matrix; and summation over repeated indices is implied. The elastic free-energy density includes all three invariants. The L2 and L3 terms are coupled by a surface relation [16], and so in the case of strong anchoring (Dirichlet boundary conditions), stationary points of the functional only depend on the sum L2 +L3 —this is all that appears in the Euler-Lagrange partial differential equations. Under more general boundary conditions, these separate terms will have a role. There is a definite “elastic constant degeneracy” here: when one makes a comparison with the Frank elastic energy—and there are various ways to do this—one comes to a situation with only two independent elastic coefficients (see, for example, [16]). The bulk free-energy density includes terms in the Landau expansion up to sixth order to allow for a stable biaxial phase in the bulk (see [13]). The
10
field term involves a general linear coupling between the Q tensor and a given tensor-field F . For a typical coupling to a magnetic field, say, we would have 1 χ H Q H 2 a α αβ β
= tr(F Q) ,
with Fαβ := 12 χa Hα Hβ .
The surface free energy acts as a finite-element penalty term coercing conformance between Q and the prescribed tensor Q0 on the boundary. The surface anchoring strength, W , is allowed to be piecewise constant. A strong anchoring (Dirichlet boundary) condition is imposed (at the discretization level) by taking appropriately large values for W . And so the problem is to minimize F (Q) over admissible tensor fields. While the free-energy functional is somewhat large and complicated, its overall structure is relatively clean. It can be written in the form F (Q) = 12 A(Q, Q) + G(Q) , where the quadratic form A(·, ·) contains only the elastic part, and the functional G contains all the remaining terms. The quadratic form A(·, ·) will be H01 -elliptic—the space referred to here is a natural tensorial analogue of the Sobolev space H01 (Ω) (see [8, 9])—provided the elastic constants satisfy certain conditions. These can be written 0 < L1 ,
−L1 < L3 < 2L1 ,
− 35 L1 −
1 L 10 3
< L2
and are essentially established in [16] using “spherical tensors.” The current version of our numerical package is for the geometry of a finite circular cylinder. We plan to modify it to handle more general cross sections, but we will continue to make essential use of the cylindrical nature of the geometry. The code allows for strong or weak anchoring conditions on the lateral surface of the cylinder, and for strong or weak anchoring or periodic conditions on the end caps. It utilizes a piecewise-linear finite-element discretization. Local minimizers are sought via a Newton iterative method. For each “outer” iteration, the Newton step is computed by a preconditioned conjugate-gradient method. The Newton method is safeguarded (globalized) by using an augmented Hessian and back-tracking. This approach was deliberately chosen to enable the extension to computing unstable stationary points and doing full bifurcation and path following; although the code at present is only capable of minimization (computing meta-stable stationary points). More detailed discussion can be found in [8] or [9]. We are continuing to try to enhance the utility and effectiveness of the code, while at the same time using it to perform numerical simulations on systems of physical interest and to compare with theory and experiments. We 11
plan to add a more general surface free-energy term. We are continuing to experiment with different preconditioners to enhance the convergence of the inner conjugate-gradient iterations. We are moving toward full bifurcation capabilities. We are looking for ways to reduce the dimensionality of some of the calculations by considering nonlinear “Krylov Subspace” methods (see, for example, [4]). A major barrier to the usefulness of such a code for practical purposes is the degeneracy of the elastic constants; we are evaluating ways to overcome this by adding additional terms to the elastic free-energy density.
3.2
Application: 2-D phases in a cylinder
In recent years, a lot of attention has been devoted to studying liquid crystal structures in sub-micron-size cylindrical cavities (see, for example, [1, 5, 14, 17, 18] and references contained therein). We are in the process of using the Landau-deGennes minimization code discussed in §3.1 to numerically study these systems. As a first effort, and to reduce the number of free parameters, we have taken L2 = L3 = 0 (which corresponds to an “equal elastic constants” approximation) and D = E = E 0 = 0 in (2), and have non-dimensionalized as follows: F(Q) =
Z n ω
1 2
Qαβ,γ Qαβ,γ + 12 (T − 1)tr(Q2 ) − 13 B tr(Q3 ) 2 2
1 4
+ C tr(Q )
o
Z
1 2
+ W ∂ω
tr (Q − Q0 )2 . q
Here we have used A = A0 (T − 1) and have taken X := L1 /A0 as our length scale: xα := xα /X. The region ω is a disk of radius R. The free-energy density per unit length in the x3 direction is expressed in units of L1 ( F := F /L1 ), and the remaining coefficients are defined by B :=
B , A0
C=
C , A0
W W := √ . L1 A0
We have considered the case B = C = W = 1.0. For these values of the parameters, the critical temperatures in the bulk are given in Table 1. For the boundary conditions, we have utilized a form of homeotropic anchoring condition. We use a uniaxial tensor of the form q
Q0 = 12 S (3 rb ⊗ rb − I) ,
S=
1+
1 − 24(T − 1) . 6
Here rb is the unit vector in the radial direction, and the latter equation above expresses the bulk (S, T ) relationship for this model. 12
Super Cooling: Nematic-Isotropic Transition: Super-Heating:
T = 1.0 . T = 28/27 = 1.037 . T = 25/24 = 1.042
Table 1: Bulk critical temperatures for the case B = C = 1.0. We find three distinct structural phases co-existing for various values of R and T . These are depicted in Figure 4. In addition to the “escaped radial” configuration of Cladis and Kl´eman [5], we find two structures with line defects. These are familiarly referred to as the “planar polar”—referred to as ˇ “planar polar with line defects (PPLD)” by Zumer [14]—and “planar radial” configurations. The escaped-radial configuration is essentially uniaxial throughout: the eigenvalues λ3 and λ2 agree at the level of the discretization error and are essentially indistinguishable from each other in the plots. The degree of order varies little from the bulk value imposed at the boundary; although there is some relaxation of order indicated in the core. The defect structures exhibit a type of “local biaxial transition” in the vicinity of the disclinations (as observed by Schopohl and Sluckin in [20]). We see this “eigenvalue exchange” mechanism commonly in various systems (see [12, 19]). We have constructed a crude phase diagram (Figure 5), which indicates regions (as a function of R and T ) where these various phases are stable or meta-stable. It agrees qualitatively with other evidence, in that it indicates the preference for the escaped solution at sufficiently large radii, the existence of the planar radial configuration only for very small radii, and the stability of the rotational-symmetry-breaking planar polar configuration for intermediate values of the radius of the disk. The configuration referred to as “isotropic” (ISO) is an essentially isotropic phase in which the Q tensor rapidly decays to zero away from the boundary, where an ordered phase condition is imposed. There is another structure that we encounter in parameter ranges with weaker surface anchoring strengths; it is the so-called “axial” phase [17], in which the director breaks free of the surface anchoring influence and forms a uniform uniaxial phase, oriented in the direction of the axis of the cylinder. Much work remains to be done on this study. We are in the process of performing a more thorough exploration of the dependence of these structures on these various parameters. We also intend to more carefully examine the natures of the transitions between the various phases and to establish some connection with data for some common liquid-crystal materials—with respect
13
Axial View of Planar Polar Field ( T = 0.98 )
20
20
20
10
10
10
0
0
0
y
30
−10
−10
−10
−20
−20
−20
−30
−30 −30
−20
−10
0 x
10
20
30
−30 −30
−20
−10
0 x
10
20
30
−30
0
0
0
z
10
−10
−10 −20
−10
0 x
10
20
30
−20
−10
0 x
10
20
30
−30
0.3
0.3
0.2
0.2
0.2
Eigenvalues
0.3
Eigenvalues
0.4
0.1
0
0
−0.1
−0.1
−0.1
−0.2 −30
−0.2 −30
10
20
30
20
30
−20
−10
0 x
10
−10
0 x
10
20
30
0.1
0
0 x
10
Planar Radial Field ( T = 0.98 )
0.4
−10
−20
Planar Polar Field ( T = 0.98 )
0.4
−20
0 x
−10 −30
Escaped Radial Field ( T = 0.98 )
0.1
−10
Side View of Planar Radial Field
10
−30
−20
Side View of Planar Polar Field
10
z
z
Side View of Escaped Radial Field
Eigenvalues
Axial View of Planar Radial Field ( T = 0.98 )
30
y
y
Axial View of Escaped Radial Field ( T = 0.98 ) 30
20
30
−0.2 −30
−20
−10
0 x
10
20
30
Figure 4: Disk phases (down-sampled, R = 30, T = 0.98): “escaped radial,” “planar polar,” “planar radial,” (a) top view, (b) side view, and (c) eigenvalues of Q tensor along x axis. to material parameters and experimental results on such systems. The code is performing well, and we anticipate submitting these results soon.
Acknowledgment We are grateful to Rein Zeller (Hillsdale College) for coding several modules that were used in a preliminary version of the Frank minimization code. The work on the stripe phases has been in collaboration with David Allender (Department of Physics and Liquid Crystal Institute, Kent State University) and Dick Hornreich (Weizmann Institute). Much of the work on the Landau-deGennes minimization code and accompanying application has been initiated in the dissertation of Tim Davis (Department of Mathematics, Case Western Reserve University) [8]. This work is in collaboration with him and with Peter Palffy-Muhoray (Liquid Crystal 14
Phase Diagram ( Bbar = Cbar = Wbar = 1.0 ) 100 ER ( PP ) 90 80 ISO 70 Radius (Rbar)
PP ( ER ) 60 50 PP ( ER, PR ) 40 30 20 PP ( PR ) 10 0 0.9
PR 0.92
0.94
0.96 0.98 1 Reduced Temperature (Tbar)
1.02
1.04
Figure 5: Phase diagram for disk phases (L2 = L3 = 0, B = C = W = 1): “escaped radial” (ER), “planar polar” (PP), “planar radial” (PR), and “isotropic” (ISO). Radius (R) vs temperature (T ), in reduced units. Metastable phases in parentheses. Institute and Department of Physics, Kent State University). We are grateful to Sami Mkaddem for assistance with the numerical results and graphics reported in §3.2.
References [1] D. W. Allender, G. P. Crawford, and J. W. Doane. Determination of the liquid-crystal surface elastic constant k24 . Physical Review Letters, 67(11):1442–1445, 1991. [2] D. W. Allender, R. M. Hornreich, and D. L. Johnson. Theory of the stripe phase in bend-Fr´eedericksz-geometry nematic films. Physical Review Letters, 59(23):2654–2657, December 1987. [3] F. Alouges. An energy-decreasing algorithm for harmonic map. In J.-M. Coron, J.-M. Ghidaglia, and F. H´elein, editors, Nematics: Mathematical and Physical Aspects, volume 332 of NATO ASI Series. Series C: Mathematical and Physical Sciences, pages 1–13, Dordrecht / Boston / London, 1991. NATO, Kluwer Academic Publishers. [4] P. N. Brown and Y. Saad. Projection methods for solving nonlinear systems of equations. In J.-M. Coron, J.-M. Ghidaglia, and F. H´elein, edi15
tors, Nematics: Mathematical and Physical Aspects, volume 332 of NATO ASI Series. Series C: Mathematical and Physical Sciences, pages 341–355, Dordrecht / Boston / London, 1991. NATO, Kluwer Academic Publishers. [5] P. E. Cladis and M. Kl´eman. Non-singular disclinations of strength s = +1 in nematics. Journal de Physique, 33(5–6):591–598, May–June 1972. [6] R. Cohen, R. Hardt, D. Kinderlehrer, S.-Y. Lin, and M. Luskin. Minimum energy configurations for liquid crystals: Computational results. In J. L. Ericksen and D. Kinderlehrer, editors, Theory and Applications of Liquid Crystals, volume 5 of IMA Volumes in Mathematics and its Applications, pages 99–122, New York, 1987. Institute for Mathematics and Its Applications, Springer-Verlag. [7] R. Cohen, S.-Y. Lin, and M. Luskin. Relaxation and gradient methods for molecular orientation in liquid crystals. Computer Physics Communications, 53:455–465, 1989. [8] T. A. Davis. Finite-Element Analysis of the Landau-deGennes Minimization Problem for Liquid Crystals in Confinement. PhD thesis, Kent State University, August 1994. [9] T. A. Davis and E. C. Gartland, Jr. Finite-element analysis of the landaudeGennes minimization problem for liquid crystals. In preparation. [10] J. L. Ericksen. Liquid crystals with variable degree of orientation. Archive for Rational Mechanics and Analysis, 113:97–120, 1991. [11] R. Fletcher. Practical Methods of Optimization. Wiley, New York, 2nd edition, 1987. [12] E. C. Gartland, Jr., P. Palffy-Muhoray, and R. S. Varga. Numerical minimization of the Landau-de Gennes free energy: Defects in cylindrical capillaries. Molecular Crystals and Liquid Crystals, 199:429–452, 1991. Proceedings of the 13th International Liquid Crystal Conference Vancouver, July, 1990. [13] E. F. Gramsbergen, L. Longa, and W. H. de Jeu. Landau theory of the nematic-isotropic phase transition. Physics Reports (Review Section of Physics Letters), 135(4):195–257, 1986. ˇ [14] S. Kralj and S. Zumer. Saddle-splay elasticity of nematic structures confined to a cylindrical capillary. Physical Review E, 51(1):366–379, January 1995.
16
[15] S.-Y. Lin and M. Luskin. Relaxation methods for liquid crystal problems. SIAM Journal on Numerical Analysis, 26(6):1310–1324, December 1989. [16] L. Longa, D. Monselesan, and H.-R. Trebin. An extension of the LandauGinzburg-de Gennes theory for liquid crystals. Liquid Crystals, 2(6):769– 796, 1987. [17] R.-M. Marroum, G. S. Iannacchione, D. Finotello, and M. A. Lee. Numerical study of cylindrically confined nematic liquid crystals. Physical Review E, 51(4):R2743–R2746, April 1995. [18] R. B. Meyer. On the existence of even indexed disclinations in nematic liquid crystals. Philosophical Magazine, 27:405–424, 1973. [19] P. Palffy-Muhoray, E. C. Gartland, Jr., and J. R. Kelly. A new configurational transition in inhomogeneous nematics. Liquid Crystals, 16(4):713– 718, 1994. [20] N. Schopohl and T. J. Sluckin. Defect core structure in nematic liquid crystals. Physical Review Letters, 59(22):2582–2584, 1987.
17
Abstract Discussed are two numerical packages for computing equilibrium structures of liquid crystals. The first uses Frank elasticity in a rectangular geometry in two spatial dimensions. The second uses a tensor order parameter and Landau-deGennes free energy in the geometry of a finite circular cylinder. The Frank code is being used to study the phenomenon of the “stripe phase” in thin nematic films in the bendFr´eedericksz geometry. The cylinder code is being used to investigate configurations in liquid crystals confined in microscopic cylindrical pores subject to homeotropic anchoring conditions. Preliminary results of these projects are presented.
1
Introduction
The Liquid Crystal Institute at Kent State University is the locus of the Science and Technology Center for “Advanced Liquid Crystalline Optical Materials” (ALCOM), which is funded by the National Science Foundation and the State of Ohio. This Center is principally concerned with the development (analysis, characterization, and synthesis) of new liquid crystalline materials for applications in various kinds of display devices. The effort is broad and inter-disciplinary; it involves researchers in physics, chemistry, mathematics, ∗
This research was partially supported by the National Science Foundation Science and Technology Center on “Advanced Liquid Crystalline Optical Materials” (ALCOM) under grant DMR 89-20147 and by NSF grant DMS-9310733.
1
and computer science, working on problems from basic theory to practical application, in programs involving theory, experiment, and computation. In recent years, we have been involved in various activities of numerical computation in support of this effort. Most of this has been concerned with exploring the static-equilibrium structures of the molecular orientational properties of liquid crystal materials in confinement. Of particular interest are the geometries of droplets, cylindrical pores, and films. We give here a preliminary report on two projects in this area, both of which are currently in progress. The first project uses Frank elasticity to study the development of so-called “stripe” or “ripple” phases in nematic thin films in the presence of a magnetic field at temperatures just above the nematic to smectic-A transition point. The divergence of the bend elastic constant in such a parameter range favors the development of a translational-symmetry-breaking configuration with periodicity in the plane of the film, perpendicular to the applied field. These textures have been observed experimentally, and there exists some theory, upon which our computational results build. The second project uses the tensor order parameter and Landau-deGennes free energy to study co-existing structural phases of materials confined in microscopic cylindrical pores subject to homeotropic anchoring conditions. In addition to the familiar “escaped” radial director field, two other configurations are found. One preserves the rotational symmetry of the problem; while the other one breaks it. Both involve line defects: one of unit strength running down the axis of the cylinder, the other possessing a pair of half-strength disclinations. The results are consistent with experimental evidence and other numerical simulations.
2
Frank minimization
We have developed a general numerical code to compute equilibrium director configurations in a rectangular geometry in two spatial dimensions by minimizing a discretized Frank free energy. We have adapted the code to numerically study “stripe” phases in thin nematic films. We report on aspects of the code and on the results obtained so far on the study of the stripe phases.
2.1
Model and code
The free energy utilized is b = F (n)
1 Z dZ b b dxdy , f (n) 2 c a
2
with the density given by the standard Oseen-Z¨ ocher-Frank elastic model with magnetic field: b := K1 (∇ · n) b 2 + K2 (n b · (∇ × n)) b 2 f (n) b × (∇ × n)| b 2 − χa (H · n) b 2 . + K3 |n
(1)
At present, boundary conditions are allowed to be Dirichlet or free (on any of the four sides) or periodic in the x direction and/or periodic in the y direction. We are in the process of adding a general surface free-energy density and including the K24 “saddle-splay” elastic constant term (as a surface term). In the future, we plan to modify the code to include a scalar order parameter, s(x, y), to allow for a locally variable degree of order (along the lines of the treatment by Ericksen in [10]). The problem is discretized by introducing a rectangular grid, which is not required to be uniform in either direction. The discrete director is represented at each mesh point by a unit vector. The free-energy functional is approximated by a composite mid-point quadrature rule, using cell-centered finite differences and averages. Local minimizers are found by using a projectedconjugate-gradient line search procedure similar to that employed in [7]. The minimization algorithm proceeds in the following way. At each step of the iteration, the gradient of the discrete free energy is computed and projected onto the tangent space of the constraint set at that configuration—this process is local in nature, involving projecting threevectors onto tangent planes to unit spheres. A search direction is built up from these using the Fletcher-Reeves conjugate gradient recurrence relations (see, for example, [11]). A one-dimensional search is then conducted in this direction but always remaining on the constraint set. The approach differs from other published approaches (e.g., [3, 6, 7, 15]) in this regard. Those schemes all perform a true “line” search, allowing the discrete director field to move off of the constraint set, violating its pointwise unit-length constraint, and then re-normalizing all the vectors once the line search is terminated. In all of the situations in which these other codes have been employed, this re-normalization step was provably free-energy reducing [3]. For our application below (and for others we envision), this is not the case. The structures and configurations we seek are also rather fragile, and so we feel that this greater fidelity to the constraints should enhance the robustness of the code. Some additional features have been added to the code to improve its utility and performance. These include a primitive multi-level capability, whereby successive refinements of an initial mesh can be generated by cell-edge bisections. Computed approximate solution fields are interpolated to finer meshes via a bi-cubic spline interpolation routine. The discrete free-energy functional 3
admits an asymptotic error expansion in even powers of the cell dimensions, and this makes extrapolation a very effective procedure to achieve highly accurate free-energy approximations (which are needed for producing good phase diagrams in the ill-conditioned situations, such as near Fr´eedericksz transitions). A simple, single-parameter continuation procedure is also implemented. It uses an application-specific routine to update the problem parameters in terms of changes in a single path parameter. This allows one to follow an arbitrary trajectory in the parameter space. Initial guesses for the unknown discrete director field are generated from a secant-line approximation through the two most recently computed solution fields. The approach is rather simplistic, but it is helpful in collecting data for phase diagrams and the like. The code is written in Fortran double precision and developed for a workstation environment. In its current incarnation (for a Hewlett-Packard workstation), it utilizes HP Matrix-Vector Library routines as well as routines from the NAG Software Library. Matlab is used for the graphics and visualizations. The code was benchmarked against a uniform Fr´eedericksz distortion computed independently by iteratively solving the implicit analytical solution in terms of elliptic integrals. For a realistic computation of 64 × 64 mesh cells, there are on the order of 4,000 nodes (and constraints) and 12,000 unknowns. With no preconditioning, the code can take 3–5 minutes to converge to the level of the discretization error from a moderately good initial guess. The accuracy of the approximate free energy is of the order of four decimal digits (in absolute error terms); the accuracy of the discrete director field is of the order of two digits. In continuation mode, with high-accuracy initial guesses available, convergence is much faster. We have built in a very simple preconditioner that utilizes rank-one perturbations of the identity. This can be useful, for instance, when collecting data near a point where one has some analytical (or intuitive) information about a mode that is losing stability (e.g., near a Fr´eedericksz transition). This can help dramatically in some cases, but it is not general, and it is of limited applicability. There are several positive features that we can point to in such a code. It is relatively easy to construct: it just requires a discrete free energy and its gradient. No second-derivative (Hessian) information is used, and so no largesparse matrix linear algebra and data structures are required. When dealing with periodic (or doubly periodic) conditions, no “phase condition” is required: the gradients (and projected gradients and conjugate projected gradients) are all guaranteed to be “orthogonal to the periodicity.” Negative features include the relatively low order of the discretization
4
method. Though helped by extrapolation, it makes it very difficult to compute fragile structures on coarse meshes. The averaging used in the cell-centered difference formulas also admits some troublesome non-physical solutions—these typically have a high-frequency micro-structure at the level of the mesh. Finally we mention that the lack of second-derivative information makes it impossible to compute (or path follow) unstable stationary points; this limits the usefulness of the code for full bifurcation studies.
2.2
Application: stripes in nematic thin films
An outstanding problem in the physics of nematic films is the stripe phase that appears in the bend-Fr´eedericksz geometry at temperatures slightly above that of the nematic to smectic-A phase transition in a sufficiently strong magnetic field. Here we have a material confined to a cell 0 < x < d of assumed infinite extent in the y and z directions, satisfying homeotropic anchoring conditions b = x) b at the plate boundaries x = 0 and x = d. A magnetic field is applied (n in the plane of the film, in the z direction: H = H zb . The problem geometry is depicted in Figure 1. z
y H
l periodic
x d
Figure 1: Geometry and boundary conditions for the bend-Fr´eedericksz transition and competing stripe phase solution. In [2], a linear stability analysis was used with a two-parameter ansatz to show that at magnetic fields slightly above the critical Fr´eedericksz threshold, a phase with periodicity in the y direction has a lower free energy than the uniform Fr´eedericksz distortion. We have used the code discussed in §2.1 to numerically study this problem. We model the system as follows. Take as free energy 1 Z lZ d b l) = b dxdy , F (n, f (n) 2l 0 0
5
b given as before by (1). Thus F is a free-energy density per unit area, with f (n) the free-energy density per unit length in z of one period in y (normalized by b and l, subject the length of the period). We seek to minimize F over both n b =x b on x = 0 and x = d, and periodic to homeotropic anchoring conditions n conditions on y = 0 and y = l. We then compare the free-energy density of this single periodic cell with that of the Fr´eedericksz configuration, which is uniform in the y direction. We non-dimensionalize as in [2]: b l) = F (n,
Z
1/2Z 1
0
b l) dxdy , f (n,
0
where
2
b l) := M ∇ · n b f (n,
2
b · ∇×n b + MN n
2
b × ∇×n b − h2 π 2 n2z . + n
Here x := x/d, y := y/2l, l := 2l/d (twice the period, in units of the cell thickness, or the ratio of the y scaling to the x scaling), and !
∇ :=
∂ 1 ∂ , . ∂x l ∂y
F is in units of K3 /d. The parameters M, N, and h are defined by K1 M := , K3
K2 N := , K1
Hd h := π
s
χa , K3
this last representing the magnetic field in units of the critical field (i.e., h = 1.0 corresponds to the critical Fr´eedericksz threshold). We have fixed the computational domain of the problem on (0, 1) × (0, 1/2) because previous experiments and analyses have indicated that the period of the stripe phases is in the range of .3 to .7 times the cell thickness. We are interested in the range where M is small: the bend elastic constant (K3 ) diverges relative to the splay constant (K1 ) as TN A is approached from the right. The twist constant (K2 ) also diverges, but at a much more moderate rate, and so N is treated as O(1). b l) (on a fixed computational domain) Our problem then is to minimize F (n, b and l, subject to the appropriate boundary and pewith respect to both n riodic conditions (and the pointwise unit-length constraint). The free-energy functional F has the structure 1 1 b l) = F 0 (n) b + F 1 (n) b b 2. F (n, + F 2 (n) l l 6
We perform the numerical minimization by a “decoupling iteration,” first minb (using the projected-conjugate-gradient line search imizing with respect to n discussed above), and then minimizing with respect to l via the exact relation lν+1 = −2
h b hν+1 F2 n h
b hν+1 F1 n
,
properly safeguarded—repeating until convergence. This approach can be thought of as a type of fractional-step method. It can also be viewed as a form of nonlinear block Gauss-Seidel. It has the advantage that it is relatively easy to implement on top of a basic existing Frank-minimization code, such as we have already developed. The main results we have obtained thus far include the verification of the first-order nature of the phase transition (from the uniform Fr´eedericksz distortion to the periodic stripe phase); there is a definite overlap in the regions of meta-stability, which can be relatively large in certain ranges of parameters. We have also developed a refinement of the phase diagram (Figure 2) from [2]. It retains the qualitative features of the one in that study but indicates the global stability of the stripe phase at higher elastic-constant ratios M. Stripe Regime Thresholds 1.1 1.09
original new
Scaled Applied Field ( h )
1.08 1.07 N = .75
1.06
N = .50
1.05 1.04
N = .25
1.03 1.02 1.01 1 0
0.02
0.04 0.06 0.08 0.1 0.12 Elastic Constants Ratio ( M = K1 / K3 )
0.14
0.16
Figure 2: Stripe threshold field, original vs new, N = .25, .50, .75. For the indicated values of N, the stripe phase has lower free energy than the uniform Fr´eedericksz distortion for values of M to the left of the threshold (crossover) lines. The most significant result we have established so far is the fact that these stripe phases are actually stable below the critical Fr´eedericksz threshold. Figure 3 shows representative structures for stripe phases at three different field 7
M = 0.126, N = 0.25, h = 1.03
M = 0.14, N = 0.25, h = 1.07
M = 0.09, N = 0.25, h = 0.99
0.8 0.7 0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0
−0.1
−0.1
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.2
M = 0.14, N = 0.25, h = 1.07
0.4
0.6
0.8
1
0
0.2
0.1
0.1
0.1
0
0
0
−0.1
−0.1
−0.1
−0.2
−0.2 0.2
0.4
0.6
0.8
0
0.2
0.4
0.6
0.8
1
0.2
0.2
0.2
0.1
0.1
0.1
0
0
0
−0.1
−0.1
−0.1
−0.2 0
0.2
0.4
0.6
0.8
0
0.2
M = 0.14, N = 0.25, h = 1.07
0.4
0.6
0.8
1
1
0.5
0.15
0.2
0.25 y
0.3
0.35
0.4
0.45
0.5
0.05
0.1
0.15
0.2
0.25 y
0.3
0.35
0.4
0.45
0 0
0.5
0
0
0
−0.5
−0.5
−0.5
20
25
30
1
0.4
0.6
0.8
1
0.05
0.1
0.15
0.2
0.25 y
0.3
0.35
0.4
0.45
0.5
A
1 0.5
A
1 0.5
A
1
15 frequency
0.8
1.5
0.5
10
0.6
Discrete Fourier Cosine Coefficients
1.5
5
0.2
Discrete Fourier Cosine Coefficients
1.5
0
0
0.4
0.5
0 0
Discrete Fourier Cosine Coefficients
−1
0.2
theta
1 theta
1 theta
1.5
0.1
1
M = 0.09, N = 0.25, h = 0.99
1.5
0.05
0
M = 0.126, N = 0.25, h = 1.03
1.5
0 0
0.8
−0.2
1
0.5
0.6
−0.2
1
−0.2
0.4
M = 0.09, N = 0.25, h = 0.99
0.2
0
0.2
M = 0.126, N = 0.25, h = 1.03
0.2
−1
0
5
10
15 frequency
20
25
30
−1
0
5
10
15 frequency
20
25
30
Figure 3: Stripe phases (down-sampled from 64 × 63 discretization), N = .25, (M, h) = (.140, 1.07), (.126, 1.03), (.090, 0.99), (a) top view, (b) side view vs (c) uniform Fr´eedericksz distortion, (d) vertical angle in cell center (as a function of y, vs θm for the uniform Fr´eedericksz distortion) and (e) its discrete Fourier cosine coefficients. strengths: 7% and 3% above critical (h = 1.07, 1.03) and 1% below (h = 0.99). At high fields (where there is little “room to maneuver”), the stripe phase appears as a modest sinusoidal perturbation on top of the uniform Fr´eedericksz distortion. At lower field strengths (in particular, below the critical threshold), the stripe phase uses relatively cheap splay and twist deformations to saturate with the applied field to a high degree in the core of the cell. Negative results include the fact that we are unable to obtain any of these stripe configurations on meshes coarser than 32 × 32. Thus we have not been able to take full advantage of multi-level computing and extrapolation. This would appear to be related to the relative fragility of these structural phases and to the low order of accuracy of the discretization model. Our path8
following steps also were required to be relatively short, as the iterations would readily become attracted to spurious, non-physical solutions.
3
Landau-deGennes minimization
We have developed a numerical code to compute equilibrium configurations of the tensor order parameter in the geometry of a circular cylinder by minimizing a discretized Landau-deGennes free energy [8, 9]. We are using this code numerically to study structural phases in microscopic cylindrical pores. We report on aspects of the code and on the results we have obtained so far on the study of the phases in a cylinder.
3.1
Model and code
We model our system in terms of the tensor order parameter, Q. Thus Q(x) is a symmetric, traceless tensor field defined (at each point x ∈ Ω) by Qαβ :=
1 2
< 3lα lβ − δαβ > ,
α, β = 1, 2, 3 ,
where < · > denotes an average with respect to the orientational probability distribution function at the point, and lα are the components of a unit vector in the direction of the long axis of a molecule. The Q tensor contains information about the preferred direction(s) of orientation as well as the degree of order. This information is most readily interpreted in terms of the eigenvalues and eigenvectors of Q : (λi , ui ), i = 1, 2, 3. Ordering these real eigenvalues from right to left, they satisfy − 12 ≤ λ3 ≤ λ2 ≤ λ1 ≤ 1 . The associated eigenvectors can be taken to be orthonormal and give extremals of the mean-square projection of the axis of a molecule: n
o
2λi + 1 . 3 The method we currently use to visualize the Q tensor (see Figure 4) is to draw an ellipsoid, with axes in the directions of the eigenvectors ui , scaled relatively by the mean-square-projection factors (2λi + 1)/3 (the values of which are in [0, 1]). Special (limiting, extreme) cases are max < (u · l)2 > : u ⊥ {u1 , . . . , ui−1 } , |u| = 1 = < (ui · l)2 > =
perfectly ordered with a positive order parameter: λ3 = λ2 = − 12 ,
λ1 = 1 ,
equivalently a3 = a2 = 0, a1 = 1, in which case the ellipsoid degenerates to a line segment (in the u1 direction); 9
perfectly ordered with a negative order parameter: λ3 = − 12 ,
λ2 = λ1 = 14 ,
or a3 = 0, a2 = a1 = 1/2, in which case the ellipsoid degenerates to a disk (perpendicular to the u3 direction); isotropic: λ3 = λ2 = λ1 = 0 , or a3 = a2 = a1 = 1/3, in which case the ellipsoid is a sphere. The material is said to be “biaxial” (at a point) if λ3 < λ2 < λ1 . It is “uniaxial” if λ3 = λ2 < λ1 (positive order parameter) or λ3 < λ2 = λ1 (negative order parameter). We utilize a Landau-deGennes free energy of the following form: F (Q) =
Z Ω
{ felastic (∂Q) + fbulk (Q) + ffield (Q) } +
Z ∂Ω
fsurf (Q) ,
where the densities are given (in compact tensor notation) by L1 Qαβ,γ Qαβ,γ + 12 L2 Qαβ,β Qαγ,γ + 12 L3 Qαβ,γ Qαγ,β A tr(Q2 ) − 13 B tr(Q3 ) + 14 C tr(Q2 )2 + 15 D tr(Q2 )tr(Q3 ) + 16 E tr(Q2 )3 + 16 E 0 tr(Q3 )2 ffield (Q) := − tr(F Q)
felastic (∂Q) := fbulk (Q) :=
fsurf (Q) :=
1 2 1 2
1 2
(2)
W tr (Q − Q0 )2 .
Here L1 , L2 , . . . , are material-dependent constants; Qαβ,γ := ∂Qαβ /∂xγ ; tr(·) denotes the ordinary trace of a matrix; and summation over repeated indices is implied. The elastic free-energy density includes all three invariants. The L2 and L3 terms are coupled by a surface relation [16], and so in the case of strong anchoring (Dirichlet boundary conditions), stationary points of the functional only depend on the sum L2 +L3 —this is all that appears in the Euler-Lagrange partial differential equations. Under more general boundary conditions, these separate terms will have a role. There is a definite “elastic constant degeneracy” here: when one makes a comparison with the Frank elastic energy—and there are various ways to do this—one comes to a situation with only two independent elastic coefficients (see, for example, [16]). The bulk free-energy density includes terms in the Landau expansion up to sixth order to allow for a stable biaxial phase in the bulk (see [13]). The
10
field term involves a general linear coupling between the Q tensor and a given tensor-field F . For a typical coupling to a magnetic field, say, we would have 1 χ H Q H 2 a α αβ β
= tr(F Q) ,
with Fαβ := 12 χa Hα Hβ .
The surface free energy acts as a finite-element penalty term coercing conformance between Q and the prescribed tensor Q0 on the boundary. The surface anchoring strength, W , is allowed to be piecewise constant. A strong anchoring (Dirichlet boundary) condition is imposed (at the discretization level) by taking appropriately large values for W . And so the problem is to minimize F (Q) over admissible tensor fields. While the free-energy functional is somewhat large and complicated, its overall structure is relatively clean. It can be written in the form F (Q) = 12 A(Q, Q) + G(Q) , where the quadratic form A(·, ·) contains only the elastic part, and the functional G contains all the remaining terms. The quadratic form A(·, ·) will be H01 -elliptic—the space referred to here is a natural tensorial analogue of the Sobolev space H01 (Ω) (see [8, 9])—provided the elastic constants satisfy certain conditions. These can be written 0 < L1 ,
−L1 < L3 < 2L1 ,
− 35 L1 −
1 L 10 3
< L2
and are essentially established in [16] using “spherical tensors.” The current version of our numerical package is for the geometry of a finite circular cylinder. We plan to modify it to handle more general cross sections, but we will continue to make essential use of the cylindrical nature of the geometry. The code allows for strong or weak anchoring conditions on the lateral surface of the cylinder, and for strong or weak anchoring or periodic conditions on the end caps. It utilizes a piecewise-linear finite-element discretization. Local minimizers are sought via a Newton iterative method. For each “outer” iteration, the Newton step is computed by a preconditioned conjugate-gradient method. The Newton method is safeguarded (globalized) by using an augmented Hessian and back-tracking. This approach was deliberately chosen to enable the extension to computing unstable stationary points and doing full bifurcation and path following; although the code at present is only capable of minimization (computing meta-stable stationary points). More detailed discussion can be found in [8] or [9]. We are continuing to try to enhance the utility and effectiveness of the code, while at the same time using it to perform numerical simulations on systems of physical interest and to compare with theory and experiments. We 11
plan to add a more general surface free-energy term. We are continuing to experiment with different preconditioners to enhance the convergence of the inner conjugate-gradient iterations. We are moving toward full bifurcation capabilities. We are looking for ways to reduce the dimensionality of some of the calculations by considering nonlinear “Krylov Subspace” methods (see, for example, [4]). A major barrier to the usefulness of such a code for practical purposes is the degeneracy of the elastic constants; we are evaluating ways to overcome this by adding additional terms to the elastic free-energy density.
3.2
Application: 2-D phases in a cylinder
In recent years, a lot of attention has been devoted to studying liquid crystal structures in sub-micron-size cylindrical cavities (see, for example, [1, 5, 14, 17, 18] and references contained therein). We are in the process of using the Landau-deGennes minimization code discussed in §3.1 to numerically study these systems. As a first effort, and to reduce the number of free parameters, we have taken L2 = L3 = 0 (which corresponds to an “equal elastic constants” approximation) and D = E = E 0 = 0 in (2), and have non-dimensionalized as follows: F(Q) =
Z n ω
1 2
Qαβ,γ Qαβ,γ + 12 (T − 1)tr(Q2 ) − 13 B tr(Q3 ) 2 2
1 4
+ C tr(Q )
o
Z
1 2
+ W ∂ω
tr (Q − Q0 )2 . q
Here we have used A = A0 (T − 1) and have taken X := L1 /A0 as our length scale: xα := xα /X. The region ω is a disk of radius R. The free-energy density per unit length in the x3 direction is expressed in units of L1 ( F := F /L1 ), and the remaining coefficients are defined by B :=
B , A0
C=
C , A0
W W := √ . L1 A0
We have considered the case B = C = W = 1.0. For these values of the parameters, the critical temperatures in the bulk are given in Table 1. For the boundary conditions, we have utilized a form of homeotropic anchoring condition. We use a uniaxial tensor of the form q
Q0 = 12 S (3 rb ⊗ rb − I) ,
S=
1+
1 − 24(T − 1) . 6
Here rb is the unit vector in the radial direction, and the latter equation above expresses the bulk (S, T ) relationship for this model. 12
Super Cooling: Nematic-Isotropic Transition: Super-Heating:
T = 1.0 . T = 28/27 = 1.037 . T = 25/24 = 1.042
Table 1: Bulk critical temperatures for the case B = C = 1.0. We find three distinct structural phases co-existing for various values of R and T . These are depicted in Figure 4. In addition to the “escaped radial” configuration of Cladis and Kl´eman [5], we find two structures with line defects. These are familiarly referred to as the “planar polar”—referred to as ˇ “planar polar with line defects (PPLD)” by Zumer [14]—and “planar radial” configurations. The escaped-radial configuration is essentially uniaxial throughout: the eigenvalues λ3 and λ2 agree at the level of the discretization error and are essentially indistinguishable from each other in the plots. The degree of order varies little from the bulk value imposed at the boundary; although there is some relaxation of order indicated in the core. The defect structures exhibit a type of “local biaxial transition” in the vicinity of the disclinations (as observed by Schopohl and Sluckin in [20]). We see this “eigenvalue exchange” mechanism commonly in various systems (see [12, 19]). We have constructed a crude phase diagram (Figure 5), which indicates regions (as a function of R and T ) where these various phases are stable or meta-stable. It agrees qualitatively with other evidence, in that it indicates the preference for the escaped solution at sufficiently large radii, the existence of the planar radial configuration only for very small radii, and the stability of the rotational-symmetry-breaking planar polar configuration for intermediate values of the radius of the disk. The configuration referred to as “isotropic” (ISO) is an essentially isotropic phase in which the Q tensor rapidly decays to zero away from the boundary, where an ordered phase condition is imposed. There is another structure that we encounter in parameter ranges with weaker surface anchoring strengths; it is the so-called “axial” phase [17], in which the director breaks free of the surface anchoring influence and forms a uniform uniaxial phase, oriented in the direction of the axis of the cylinder. Much work remains to be done on this study. We are in the process of performing a more thorough exploration of the dependence of these structures on these various parameters. We also intend to more carefully examine the natures of the transitions between the various phases and to establish some connection with data for some common liquid-crystal materials—with respect
13
Axial View of Planar Polar Field ( T = 0.98 )
20
20
20
10
10
10
0
0
0
y
30
−10
−10
−10
−20
−20
−20
−30
−30 −30
−20
−10
0 x
10
20
30
−30 −30
−20
−10
0 x
10
20
30
−30
0
0
0
z
10
−10
−10 −20
−10
0 x
10
20
30
−20
−10
0 x
10
20
30
−30
0.3
0.3
0.2
0.2
0.2
Eigenvalues
0.3
Eigenvalues
0.4
0.1
0
0
−0.1
−0.1
−0.1
−0.2 −30
−0.2 −30
10
20
30
20
30
−20
−10
0 x
10
−10
0 x
10
20
30
0.1
0
0 x
10
Planar Radial Field ( T = 0.98 )
0.4
−10
−20
Planar Polar Field ( T = 0.98 )
0.4
−20
0 x
−10 −30
Escaped Radial Field ( T = 0.98 )
0.1
−10
Side View of Planar Radial Field
10
−30
−20
Side View of Planar Polar Field
10
z
z
Side View of Escaped Radial Field
Eigenvalues
Axial View of Planar Radial Field ( T = 0.98 )
30
y
y
Axial View of Escaped Radial Field ( T = 0.98 ) 30
20
30
−0.2 −30
−20
−10
0 x
10
20
30
Figure 4: Disk phases (down-sampled, R = 30, T = 0.98): “escaped radial,” “planar polar,” “planar radial,” (a) top view, (b) side view, and (c) eigenvalues of Q tensor along x axis. to material parameters and experimental results on such systems. The code is performing well, and we anticipate submitting these results soon.
Acknowledgment We are grateful to Rein Zeller (Hillsdale College) for coding several modules that were used in a preliminary version of the Frank minimization code. The work on the stripe phases has been in collaboration with David Allender (Department of Physics and Liquid Crystal Institute, Kent State University) and Dick Hornreich (Weizmann Institute). Much of the work on the Landau-deGennes minimization code and accompanying application has been initiated in the dissertation of Tim Davis (Department of Mathematics, Case Western Reserve University) [8]. This work is in collaboration with him and with Peter Palffy-Muhoray (Liquid Crystal 14
Phase Diagram ( Bbar = Cbar = Wbar = 1.0 ) 100 ER ( PP ) 90 80 ISO 70 Radius (Rbar)
PP ( ER ) 60 50 PP ( ER, PR ) 40 30 20 PP ( PR ) 10 0 0.9
PR 0.92
0.94
0.96 0.98 1 Reduced Temperature (Tbar)
1.02
1.04
Figure 5: Phase diagram for disk phases (L2 = L3 = 0, B = C = W = 1): “escaped radial” (ER), “planar polar” (PP), “planar radial” (PR), and “isotropic” (ISO). Radius (R) vs temperature (T ), in reduced units. Metastable phases in parentheses. Institute and Department of Physics, Kent State University). We are grateful to Sami Mkaddem for assistance with the numerical results and graphics reported in §3.2.
References [1] D. W. Allender, G. P. Crawford, and J. W. Doane. Determination of the liquid-crystal surface elastic constant k24 . Physical Review Letters, 67(11):1442–1445, 1991. [2] D. W. Allender, R. M. Hornreich, and D. L. Johnson. Theory of the stripe phase in bend-Fr´eedericksz-geometry nematic films. Physical Review Letters, 59(23):2654–2657, December 1987. [3] F. Alouges. An energy-decreasing algorithm for harmonic map. In J.-M. Coron, J.-M. Ghidaglia, and F. H´elein, editors, Nematics: Mathematical and Physical Aspects, volume 332 of NATO ASI Series. Series C: Mathematical and Physical Sciences, pages 1–13, Dordrecht / Boston / London, 1991. NATO, Kluwer Academic Publishers. [4] P. N. Brown and Y. Saad. Projection methods for solving nonlinear systems of equations. In J.-M. Coron, J.-M. Ghidaglia, and F. H´elein, edi15
tors, Nematics: Mathematical and Physical Aspects, volume 332 of NATO ASI Series. Series C: Mathematical and Physical Sciences, pages 341–355, Dordrecht / Boston / London, 1991. NATO, Kluwer Academic Publishers. [5] P. E. Cladis and M. Kl´eman. Non-singular disclinations of strength s = +1 in nematics. Journal de Physique, 33(5–6):591–598, May–June 1972. [6] R. Cohen, R. Hardt, D. Kinderlehrer, S.-Y. Lin, and M. Luskin. Minimum energy configurations for liquid crystals: Computational results. In J. L. Ericksen and D. Kinderlehrer, editors, Theory and Applications of Liquid Crystals, volume 5 of IMA Volumes in Mathematics and its Applications, pages 99–122, New York, 1987. Institute for Mathematics and Its Applications, Springer-Verlag. [7] R. Cohen, S.-Y. Lin, and M. Luskin. Relaxation and gradient methods for molecular orientation in liquid crystals. Computer Physics Communications, 53:455–465, 1989. [8] T. A. Davis. Finite-Element Analysis of the Landau-deGennes Minimization Problem for Liquid Crystals in Confinement. PhD thesis, Kent State University, August 1994. [9] T. A. Davis and E. C. Gartland, Jr. Finite-element analysis of the landaudeGennes minimization problem for liquid crystals. In preparation. [10] J. L. Ericksen. Liquid crystals with variable degree of orientation. Archive for Rational Mechanics and Analysis, 113:97–120, 1991. [11] R. Fletcher. Practical Methods of Optimization. Wiley, New York, 2nd edition, 1987. [12] E. C. Gartland, Jr., P. Palffy-Muhoray, and R. S. Varga. Numerical minimization of the Landau-de Gennes free energy: Defects in cylindrical capillaries. Molecular Crystals and Liquid Crystals, 199:429–452, 1991. Proceedings of the 13th International Liquid Crystal Conference Vancouver, July, 1990. [13] E. F. Gramsbergen, L. Longa, and W. H. de Jeu. Landau theory of the nematic-isotropic phase transition. Physics Reports (Review Section of Physics Letters), 135(4):195–257, 1986. ˇ [14] S. Kralj and S. Zumer. Saddle-splay elasticity of nematic structures confined to a cylindrical capillary. Physical Review E, 51(1):366–379, January 1995.
16
[15] S.-Y. Lin and M. Luskin. Relaxation methods for liquid crystal problems. SIAM Journal on Numerical Analysis, 26(6):1310–1324, December 1989. [16] L. Longa, D. Monselesan, and H.-R. Trebin. An extension of the LandauGinzburg-de Gennes theory for liquid crystals. Liquid Crystals, 2(6):769– 796, 1987. [17] R.-M. Marroum, G. S. Iannacchione, D. Finotello, and M. A. Lee. Numerical study of cylindrically confined nematic liquid crystals. Physical Review E, 51(4):R2743–R2746, April 1995. [18] R. B. Meyer. On the existence of even indexed disclinations in nematic liquid crystals. Philosophical Magazine, 27:405–424, 1973. [19] P. Palffy-Muhoray, E. C. Gartland, Jr., and J. R. Kelly. A new configurational transition in inhomogeneous nematics. Liquid Crystals, 16(4):713– 718, 1994. [20] N. Schopohl and T. J. Sluckin. Defect core structure in nematic liquid crystals. Physical Review Letters, 59(22):2582–2584, 1987.
17
