A theoretically informed coarse-grained Monte Carlo method is proposed for studying liquid crystals. The free energy functional of the system is described in the framework of the Landau-de Gennes formalism. The alignment field and its gradients are approximated by finite differences, and the free energy is minimized through a stochastic sampling technique. The validity of the proposed method is established by comparing the results of the proposed approach to those of traditional free energy minimization techniques. Its usefulness is illustrated in the context of three systems, namely, a nematic liquid crystal confined in a slit channel, a nematic liquid crystal droplet, and a chiral liquid crystal in the bulk. It is found that for systems that exhibit multiple metastable morphologies, the proposed Monte Carlo method is generally able to identify lower free energy states that are often missed by traditional approaches. Importantly, the Monte Carlo method identifies such states from random initial configurations, thereby obviating the need for educated initial guesses that can be difficult to formulate.

Liquid crystals (LCs) are finding a number of emerging applications in areas ranging from biomaterials1 to sensors2 and optical devices.3–5 For many of these applications, it is essential to generate detailed theoretical and computational predictions of the system’s structure and properties as a function of material characteristics, including any interfaces that may be used to control or manipulate the material’s behavior.6–8 

Liquid crystals exhibit intriguing structural features at multiple lengths scales. Defects, which represent an integral part of the structure, have characteristic dimensions in the range of 10 nm. Elastic interactions and interfacial effects can be felt over micrometer length scales. And typical device dimensions can range from micrometers to millimeters. One of the most widely used models of liquid crystal structure over such a wide range of length scales is that proposed by Landau and de Gennes in which the system is described at the level of a local tensorial order parameter. The model consists of a free energy functional which is minimized by resorting to well-known minimization techniques, e.g., by direct solution of the corresponding Euler-Lagrange equations,9 or through dynamic relaxation in terms of a time-dependent Ginzburg-Landau (GL) formalism.10–13 An important challenge in this type of work is that of finding relevant global free energy minima.

Many-particle molecular dynamics or Monte Carlo (MC) simulations of coarse-grained molecular models can also be used to predict liquid crystal structure. One can rely on advanced sampling techniques to identify equilibrium morphologies but, at present, many-body simulations are highly computationally demanding for systems exceeding characteristic dimensions on the order of a hundred nanometers.7,8,14

An alternative approach consists of adopting a continuum representation of the free energy of the liquid crystal, e.g., in terms of the Landau-de Gennes model, but to resort to a Monte Carlo method to identify the equilibrium (free-energy minima) structure of the system. Such a strategy has been particularly fruitful in the context of ordered polymeric materials, where it is often referred to as theoretically informed coarse grained simulation,15–19 and in this work, we develop it in the context of liquid crystalline systems.

To do so, we build on pioneering work by Gruhn and Hess20 and Ruhwandl and Terentjev,21 who were among the first to study LC behavior by stochastic sampling of a vectorial director field. These authors proposed stochastic algorithms to sample the configuration space of free energy models expressed in terms of a local alignment vector for uniaxial systems. More specifically, by relying on Metropolis sampling, they obtained the equilibrium director field n(x) of nematic LCs subject to an external field. Ruhwandl and Terentjev21 implemented a similar method to study the nematic director field around a spherical colloid with homeotropic anchoring. In addition to an elastic free energy, these authors included a surface anchoring term in the Hamiltonian originally proposed by Rapini and Papoular.22 A related approach was later employed by Saielli et al. to describe the three Fréedericksz transitions corresponding to different elastic distortions.23 Luckhurst added a chiral term to the Gruhn-Hess model and examined the director distribution of chiral LCs confined between slabs, including the effects of thermal fluctuations.24 More recently, Dong-Lai25 used the Gruhn-Hess model and reported on the uniaxial-biaxial (UB) transition of LCs confined between plates. While their results are generally in good agreement with previous works,26–29 they also found that the characteristic director profile for the biaxial structure is not symmetric when the splay and bend elastic constants are different (i.e., k11k33). For spherical geometries (e.g., droplets), Goyal and Denn presented a particularly important study in which they examined the effects of homeotropic anchoring30 on morphology; to do so, these authors extended the Gruhn-Hess formalism to include a saddle-splay term (k24), which is essential for descriptions of the transition from a radial to a uniform structure under spherical confinement. A common feature of the coarse-grained studies mentioned above is that the underlying free energy models were based on a vectorial representation of the nematic field.

In this paper, we propose an alternative theoretically informed MC sampling method for liquid crystals. The method differs from previous MC approaches in that the sampling is performed in terms of an orthonormal basis set representation of the tensor order parameter, thereby enabling efficient sampling in systems that exhibit biaxiality, including chiral materials, or systems with defects or suspended particles. In our approach, the free energy is written in terms of a Landau-de Gennes free energy model,31 where the local LC state is represented by its alignment tensor field. Such a tensorial description enables a detailed representation of defect cores. By expressing the alignment tensor in terms of an orthonormal basis set, it becomes possible to sample the alignment tensor parameters independently without relying on cumbersome constraints. Therefore, scalar order parameters, uniaxiality S and biaxiality η, as well as their corresponding director vectors, n and n′, are sampled uniformly.

The validity of the proposed method is addressed for three examples: (1) a LC confined in a slit14 (a so-called hybrid nematic cell (HNC)), (2) a LC droplet,4 and (3) a periodic chiral LC. For confined systems, the Landau-de Gennes free energy31 is supplemented by a Rapini-Papoular surface contribution in the case of homeotropic anchoring22 and by a Fournier-Galatola 4th order expression for the case of degenerate planar anchoring.32 The free energy functional adopted in Monte Carlo simulations is also minimized through a GL isothermal relaxation, thereby providing a reference to examine the performance of the proposed MC approach. We find that for systems that exhibit a rugged free energy landscape, the GL method often gets trapped in local free energy minima. In contrast, the MC sampling method is able to find the global free energy minimum.

We adopt a representation in which the free energy of the LC is a function of the local alignment tensor Q(x). It is defined in terms of the tensor of second moments of the local molecular orientational distribution function, ψ(u, x, t),

Q ( x ) = ( uu δ 3 ) ψ ( u , x , t ) d u ,
(1)

where u(x) is the ensemble average of molecular orientations and δ is the 3 × 3 identity matrix. At each configuration of the alignment tensor field Q(x), a probability is assigned proportional to the Boltzmann factor,

exp β F ( Q ( x ) ) ,
(2)

where β 1 = k B T ˆ (kB is Boltzmann constant and T ˆ plays the role of a temperature). The coarse-grained free energy functional is a sum of three terms: a short-range Landau polynomial expansion in Q(x)—to describe the isotropic-nematic (IN) transition—a long-range free energy density—for elastic distortions—and a surface free energy contribution,

F ( Q ) = d 3 x f L ( Q ) + f E ( Q ) + d 2 x f S ( Q ) .
(3)

The Landau free energy density is given by a phenomenological approximation obtained by truncating the expansion with respect to the invariants of Q, i.e.,

f L = A 2 1 U 3 tr Q 2 A U 3 tr Q 3 + A U 4 tr Q 2 2 ,
(4)

where parameter A sets an energy density scale for the model. A dimensionless parameter U determines the IN transition10,33,34 and the magnitude of the average scalar order parameter. For Eq. (4), we adopt Doi’s notation35,36 and consider only a unique value for the polynomial coefficients, i.e., A.

The elastic contribution is given in terms of gradients of the alignment tensor according to

f E = 1 2 L 1 Q i j x k Q i j x k + 1 2 L 5 ϵ i k l Q i j Q l j x k ,
(5)

where L1 and L5 are elastic constants and ϵijk is the Levi-Civita pseudotensor. For a uniaxial system, the elastic constants in Eq. (5) are related to the Frank-Oseen37–41 elastic constants42 by the following expressions:

L 1 = k 33 k 11 + 3 k 22 6 S 2 ,
(6)
L 5 = 2 q 0 k 22 S 2 ,
(7)

where S is the scalar order parameter and k11, k22, and k33 are non-vanishing elastic moduli corresponding to independent splay, twist, and bend modes.43 Here, L5 quantifies the chirality of the LC and the pitch pch = 2π/q0. Without loss of generality, in this work, we adopt a single elastic constant approximation, i.e., k11 = k22 = k33. We satisfy the Ericksen inequalities44 and the stability ranges derived by Gartland and Mkaddem45 and Longa et al.46 

For homeotropic anchoring, we adopt a Rapini-Papoular22,47–50 expression of the form

f S , R = 1 2 W R Q Q 2 ,
(8)

where WR represents the strength of the homeotropic anchoring and Q is the perpendicular tensor orientation preferred at the surface, defined as follows:

Q = S ν ν 1 3 δ ,
(9)

where ν(x) is the normal unit vector at any point x along the surface. To describe degenerate planar anchoring, we adopt a 4th order Fournier-Galatola expression32 

f S , F = 1 2 W F Q ¯ Q ¯ 2 + 1 4 W F Q ¯ : Q ¯ S 2 2 ,
(10)

where Q ¯ = Q + S δ / 3 , Q ¯ = p Q ¯ p is the tensor projection on the surface, and p = δνν.

The LC physical properties will set two types of characteristic length scales: the nematic coherence length ξN and the surface extrapolation length ξS; they are defined as follows:

ξ N = L 1 / A , ξ S = L 1 / W R  or  F .
(11)

The last term is also known as the Kleman-de Gennes extrapolation distance.36 Parameters in the simulations are non-dimensionalized in terms of the nematic coherence length ξN and material parameter A. The free energy functional of the system is expressed in dimensionless form as

F ˆ ( Q ) = F ( Q ) A ξ N 3 .
(12)

For simplicity, all functions and variables are assumed to be dimensionless throughout this manuscript.

A Metropolis MC scheme is used to generate a Markov chain of states. In particular, trial random updates of the alignment tensor field are proposed and accepted according to Metropolis criteria. A Markov chain of configurations is constructed by proposing transitions between an old configuration “o” and a new one “n” with probability Pprop(on). Trial transitions are accepted with probability,

P acc ( o n ) = min [ 1 , exp ( β Δ F ) ] ,
(13)

where ΔF = F(Q(n)) − F(Q(o)) is the free energy difference between the new and old configurations. Trial displacements of Q’s are proposed by recognizing that the tensor order parameter can be expressed in terms of its proper values, eigenvalues, and eigenvectors, according to43,51

Q = S n n δ 3 + η n n ( n × n ) ( n × n ) ,
(14)

where S(x) is the scalar order parameter, related to the maximum eigenvalue 2S/3, and η(x) is the biaxiality, related to the other two eigenvalues ± η S 3 . The scalar order parameter S and biaxiality η are bounded by

S [ 1 2 , 1 ]

and

η [ 1 3 ( 1 S ) , + 1 3 ( 1 S ) ] .
(15)

The eigenvectors, n and n′, corresponding to the maximum and second largest eigenvalues, define an orthonormal basis n , n , n × n for LC orientation.

According to Eq. (14), Q = Q(S, η, n, n′); therefore, a correct MC approach must sample each variable randomly and without a bias, within the constrains provided by the LC physics. Random sampling of S, η, n, and n′ involves cumbersome schemes. Recall that the alignment tensor is symmetric and traceless, and out of its nine cartesian components, only five are actual degrees of freedom. Therefore, it can be expressed in terms of an orthonormal tensor basis, as originally discussed by Hess and co-workers52 and, more recently, by Bhattacharjee et al.,11,53

Q x , t = ν 5 a ν x , t T ν .
(16)

The orthonormal basis is defined by five tensors,

T 1 = 3 2 z z ST = 3 2 δ 3 i δ 3 j δ i j / 3 , T 2 = 2 x y ST = 2 δ 1 i δ 2 j + δ 2 i δ 1 j / 2 , T 3 = 2 x z ST = 2 δ 1 i δ 3 j + δ 3 i δ 1 j / 2 , T 4 = 1 2 x x y y = 1 2 δ 1 i δ 1 j δ 2 i δ 2 j , T 5 = 2 y z ST = 2 δ 2 i δ 3 j + δ 3 i δ 2 j / 2 ,
(17)

where x, y, and z are the canonical ℜ3 basis, A ST is the symmetric-traceless projection operator, and δij is the Kronecker delta. Because the T m basis is orthonormal,

tr ( T m T n ) = T i j m T i j n = δ m n ,
(18)

thereby ensuring that the five scalar components aν of the alignment tensor represent simple projections,

a ν = tr ( Q T ν ) .
(19)

The a set provides a unique and independent way to specify a configuration of the alignment tensor field. For instance, the trace operations in the Landau free energy density in Eq. (4) are simple polynomials in a,

tr Q 2 = a 1 2 + a 2 2 + a 3 2 + a 4 2 + a 5 2 ,
(20)
tr Q 3 = a 1 3 6 3 2 a 1 a 2 2 a 4 2 + a 3 2 + a 5 2 2 + 3 2 a 4 a 3 2 a 5 2 2 + a 2 a 3 a 5 ,
(21)

and, by the Cayley-Hamilton theorem, we can compute directly higher powers of Q,

Q 3 = tr Q 2 2 Q + tr Q 3 3 δ .
(22)

In addition, since the characteristic polynomial of Q is

p ( λ ) = λ 3 tr Q 2 2 λ tr Q 3 3 ,
(23)

eigenvalues (and eigenvectors) can be obtained quickly using Vieta’s formulae for the roots of p(λ). Contributions to the free energy functional can be expressed directly in terms of a.

A Monte Carlo minimization of the free energy starts from a random distribution of the alignment tensor field in the bulk (and surfaces, as appropriate). The field is approximated numerically on a discrete mesh, defining a set of N collocation points. The Hamiltonian, F = F(ξN, ξS, a, ∇a), is obtained by numerical integration. A collocation point at x is selected randomly, in the bulk or surface, and a MC trial move is attempted by randomly displacing a(x) according to

a μ , n ( x ) = a μ , o ( x ) + δ μ ̄ ( ζ 0 . 5 ) ,
(24)

where component μ is selected at random from 1 to 5, ζ is a random number distributed uniformly between [0, 1], and δ μ ̄ is the maximum allowed displacement of aμ. We defined the MC move according to Eq. (24) to maximize the MC performance. First, for every MC move, the collocation point, the amount of perturbation, and the orthonormal basis are selected randomly; a single basis μ is modified in each MC step. In addition, the maximum displacement δ μ ̄ is adjusted periodically during the minimization to achieve an average 30% acceptance rate. During the MC move, the physical bounds for the tensor Q, Eq. (15), are always enforced. Note that only the local free energy difference at collocation point x needs be calculated, resulting in an efficient O(1) method. To minimize the free energy, we follow an annealing process that is controlled through the effective temperature T ˆ .54 The simulation starts from a high non-dimensional effective temperature T = A ξ N 3 / β in the Metropolis criterion. The system is isotropic at that temperature. The effective temperature is then progressively lowered, following an exponential decrease, until the LC system reaches equilibrium at a low temperature. For all minimizations reported in this work, 10 < β < 1012 over the total number of MC steps. A step is defined by the average number of accepted moves per collocation node, and we use between 200 000 and 500 000 accepted MC moves per collocation point.

In this work, the numerical approximation for F = F(ξN, ξS, a, ∇a) is evaluated within a finite-difference (FD) scheme on a uniform three-dimensional mesh. A weighted residual scheme is applied to the free energy functional with Dirac delta distributions as weighting functions. Additional details on the free energy density and free energy functional approximations are provided in  Appendix B. As mentioned earlier, a traditional Ginzburg-Landau relaxation approach was also used to minimize the free energy.  Appendix C delineates the method; details can be also found in the literature.12,13,48 For simple systems, the GL relaxation approach is approximately four times faster than the proposed MC algorithm. The main advantage of MC, however, is that it is able to find the global minimum of energy in cases where GL relaxation techniques get trapped in local minima. Also note that no attempts have been made in this work to accelerate the performance of the MC technique by resorting to multiple-spin, collective trial moves, or other advanced sampling techniques. It is conceivable that, with such improvements, the performance of MC methods could even surpass that of traditional techniques.

The parameters adopted in all nematic simulations are chosen to be representative of 4′-pentyl-4-cyanobiphenyl: A ≈ 1 × 105 J/m3 and k11 = k22 = k33 ≈ 6 pN, resulting in ξN = 7.15 nm. For chiral systems, we use ξ N = 40 nm and L1 = 25 pN. In the Landau-de Gennes framework, parameter U sets a value for the nematic scalar order parameter,

S bulk = 1 4 + 3 4 1 8 3 U .
(25)

Bulk systems are expected to be in the isotropic phase for values of U < 2.8 and in the nematic phase for U ≥ 2.8. We use a value of U = 5 for non-chiral LCs, which corresponds to a bulk scalar order parameter S ≈ 0.76. The value of U for the chiral LC is modified according to the phase of interest.55 

The walls induce opposite LC orientations, homeotropic at the top and degenerative planar at the bottom, thereby forming a HNC. This particular configuration forces the LC to undergo a mixed splay-bend distortion.56 The homeotropic anchoring follows a Rapini-Papoular 2nd order potential, while the degenerative planar is captured by the 4th order Fournier-Galatola potential. The anchoring strengths in these potentials define a set of extrapolation lengths that we kept equal by enforcing WR = WF = W. The results presented here encompass strengths from weak (W = 10−7 J/m2) to strong (W = 10−2 J/m2) anchoring.

We start the analysis by fixing the HNC thickness to d = 20ξn and by performing two different types of relaxations that are referred to as R-1 and R-2. For R-1, the top wall imposes infinite homeotropic anchoring, while the bottom wall has finite degenerative planar anchoring. Conversely, for R-2, the bottom wall imposes infinite planar anchoring and the top wall exhibits finite homeotropic anchoring.57 The minimum free energy for different values of anchoring strength is shown in Fig. 1 for both cases. As expected, we see a transition from zero elastic distortions, at weak anchoring, to large elastic distortions, at strong anchoring. In the figure, contours of the orientation angle between the bottom surface and the vector order parameter n are included. Overall, the results of the proposed MC method and those of simple GL minimization are in good agreement. Note that the energies predicted by MC simulations are at times slightly lower than those calculated by standard GL relaxation, suggesting that the MC approach is in some cases more effective in identifying the true global minimum.

FIG. 1.

Minimum free energy functional (Hamiltonian) as a function of anchoring strength ξS for a liquid crystal confined in a slit pore of thickness d = 20ξN. R-1 corresponds to results for infinite homeotropic anchoring at the top surface and finite planar anchoring at the bottom surface. R-2 corresponds to a slit with finite homeotropic anchoring at the top and infinite planar anchoring at the bottom. (a) Total free energy functional and (b) elastic and surface contributions. Continuous lines show results for Ginzburg-Landau minimizations, and symbols are from Monte-Carlo simulations. The color contours represent the orientation angle between the bottom surface and the vector order parameter n, where red means parallel to the bottom wall, and blue means perpendicular to the bottom wall.

FIG. 1.

Minimum free energy functional (Hamiltonian) as a function of anchoring strength ξS for a liquid crystal confined in a slit pore of thickness d = 20ξN. R-1 corresponds to results for infinite homeotropic anchoring at the top surface and finite planar anchoring at the bottom surface. R-2 corresponds to a slit with finite homeotropic anchoring at the top and infinite planar anchoring at the bottom. (a) Total free energy functional and (b) elastic and surface contributions. Continuous lines show results for Ginzburg-Landau minimizations, and symbols are from Monte-Carlo simulations. The color contours represent the orientation angle between the bottom surface and the vector order parameter n, where red means parallel to the bottom wall, and blue means perpendicular to the bottom wall.

Close modal

We next consider a situation in which we have a UB transition. Strong anchoring conditions are imposed at both walls and the slit thickness is changed. It has been predicted that a critical thickness exists at which such a transition occurs.26,29,57 A high mesh resolution is employed in this case, namely, 0.25ξN. Figure 2 summarizes our results for d = 2.5ξN, 3.0ξN, and 4.0ξN. The figure includes the scalar order parameter S, the biaxiality η, and the vector order parameter orientation (n) with respect to the bottom wall φ along the confinement direction z/d. The homeotropic wall is at z/d = 0.5, whereas the planar wall is at z/d = 0.5. For d = 4.0ξN, the value of S is always higher than η and φ varies linearly from top to bottom; there is no UB transition. At d = 3.0ξN, the system shows an increment in the biaxiality, and S and η are almost equal at the mid-plane. The orientation angle φ is now a non-linear function of z/d. Finally, for d = 2.5ξN, η is higher than S at the mid-plane. This implies that there is an exchange in the magnitude of the eigenvalues that correspond to S and η and that a UB transition has taken place. These results indicate that d≊3ξN is the critical thickness value for the UB transition. Figure 2 includes the final LC phases, at each value of the confinement. The colors represent the value of the orientational angle φ: red corresponds to orientations along the bottom wall, and blue corresponds to orientations perpendicular to the bottom wall.

FIG. 2.

Scalar order parameter S, biaxiality η, and orientational angle φ of the vector order parameter as a function of confinement (along z axis) for different slit thicknesses. Continuous lines show results from Ginzburg-Landau minimization, and symbols show results from Monte-Carlo simulations. The insets show the final LC phases, for each value of the confinement; the colors represent the value of the orientational angle φ: red is for orientations along the bottom wall, and blue is for orientations perpendicular to the walls.

FIG. 2.

Scalar order parameter S, biaxiality η, and orientational angle φ of the vector order parameter as a function of confinement (along z axis) for different slit thicknesses. Continuous lines show results from Ginzburg-Landau minimization, and symbols show results from Monte-Carlo simulations. The insets show the final LC phases, for each value of the confinement; the colors represent the value of the orientational angle φ: red is for orientations along the bottom wall, and blue is for orientations perpendicular to the walls.

Close modal

Liquid crystal droplets exhibit a wide variety of mesophases and therefore provide a stringent test on which to test models and simulation algorithms.12,13,58–63 The model parameters considered here for droplets are the same as those used for the HNC system, including the length scale for confinement, namely, r = 20ξN, where r is the droplet radius.

Figure 3 shows the free energy of LC droplets as a function of the anchoring strength for homeotropic and planar degenerate droplets. For strong homeotropic anchoring, a radial hedgehog (R) texture is observed with a saturn ring close to the core. This phase is characterized by elastic distortions in the bulk and negligible penalizations at the surface. For weak homeotropic anchoring, the radial phases are unstable and the LC adopts a uniform (U) morphology. In contrast, for degenerate planar anchoring, in the strong anchoring regime, one finds a bipolar (BP) morphology, with boojum defects at each pole. For weak anchoring, the uniform (U) phase is also the most stable configuration. The figure includes the final phases at strong and weak anchoring, where red iso-surfaces represent low values of the scalar order parameter S, and where the blue streamlines follow the vector order parameter n. The final states predicted by GL and MC relaxations are equivalent and, as before, the magnitude of the energy values differ by less than 1%, with the MC method yielding slightly lower numbers.

FIG. 3.

Minimum free energy as a function of anchoring strength ξS for liquid crystal nano-droplet. The droplet size is r = 20ξN. Results are shown for two different anchoring conditions, namely, homeotropic and degenerate planar. Continuous lines show results from a Ginzburg-Landau minimization, and symbols show results from Monte-Carlo simulations. The final phases for strong and weak anchoring are included in the figure. Red iso-surfaces represent low values of the scalar order parameter S, and the blue lines depict the vector order parameter n.

FIG. 3.

Minimum free energy as a function of anchoring strength ξS for liquid crystal nano-droplet. The droplet size is r = 20ξN. Results are shown for two different anchoring conditions, namely, homeotropic and degenerate planar. Continuous lines show results from a Ginzburg-Landau minimization, and symbols show results from Monte-Carlo simulations. The final phases for strong and weak anchoring are included in the figure. Red iso-surfaces represent low values of the scalar order parameter S, and the blue lines depict the vector order parameter n.

Close modal

Note that the slight differences between MC and GL can be explained by the fact that the MC method directly calculates the free energy by changing a at a collocation point; these changes are independent and local. The highest gradients that need to be resolved in MC are therefore given directly by the free energy densities of the Hamiltonian. On are sampled uniformly hand, GL minimizes through the Volterra derivatives of the free energy functional. Thus, any local change in the collocation point induces an averaging global effect. In addition, mathematically, there is an increase in the order of differentiation in the GL. For instance, the L1 elastic term results in a laplacian after the Volterra derivative. This fundamental difference between GL and MC will dictate the efficiency of the method when exploring the depths of free energy minima, according to the dominant long-range elastic term.

In what follows, we address a slightly different issue, namely, that of reaching a global minimum regardless of the starting point (or initial configuration). We do so in the context of blue phases for chiral liquid crystals, where traditional minimization techniques are known to get easily trapped in local energy minima. The chiral LC considered here is defined by ξ N = 40 nm and L1 = 25 pN; we explore a periodic cell of size L.

We first consider conditions that correspond to a blue phase I (BP-I); we use U = 3 and a cell size of L = 57.3ξN for a pitch of pch = 56.6ξN (N = 2L/pch = 2.026). Figure 4 shows multiple phases during the relaxation process. The MC relaxation starts from a completely random initial state. The fluctuating temperature is decreased exponentially as the minimization progresses. The final blue phase I has a O 8 symmetry.55 

FIG. 4.

Intermediate and final blue phase I of a periodic chiral liquid crystal with ξ N = 40 nm , U = 3, L = 57.3ξN, and pch = 56.6ξN. Iso-surfaces represent values of S = 0.6.

FIG. 4.

Intermediate and final blue phase I of a periodic chiral liquid crystal with ξ N = 40 nm , U = 3, L = 57.3ξN, and pch = 56.6ξN. Iso-surfaces represent values of S = 0.6.

Close modal

Figure 5 shows the evolution of the average scalar order parameter 〈S〉 and the local director 〈n〉 during a MC simulation and a GL relaxation. The average local director is defined as the dot product between the director at a point x and the average director surrounding it. In the figure, 〈n〉 contours of the initial and final phases with the corresponding S iso-surfaces are shown. For the 〈n〉 contours, red indicates an orientation perpendicular to the average director, and blue indicates a direction parallel to the average director. As the MC simulations approach the equilibrium morphology, local directors (n) fluctuate about a mean value that converges rapidly (in approximately 105 MC steps) to the final director configuration, shown in the inset of the figure (and the corresponding final value of 〈n〉). On the other hand, in the MC algorithm, the eigenvalues (S and η) undergo large changes in the early stages of the calculation and approach equilibrium after about 4 × 105 steps. In contrast, for the GL minimization, the eigenvalues (S and η) approach their equilibrium values early in the minimization (after approximately 3 × 105 time steps), whereas the local directors (and the corresponding morphology, given by n and shown in the inset) continue to change after much longer time scales (4 × 105 time steps). These differences help explain why GL gets more easily trapped in metastable states.

FIG. 5.

Average scalar order parameter 〈S〉 and local director 〈n〉 during MC and GL relaxations. Initial and final phases are included for each minimization. Note that the final state for the GL minimization does not correspond to equilibrium.

FIG. 5.

Average scalar order parameter 〈S〉 and local director 〈n〉 during MC and GL relaxations. Initial and final phases are included for each minimization. Note that the final state for the GL minimization does not correspond to equilibrium.

Close modal

We next proceed to consider different cell sizes, namely, L = 15.81ξN and L = 31.62ξN, and to change the chirality. The final morphologies predicted by GL minimization depend entirely on the initial configuration. Past literature studies of chiral liquid crystals using traditional techniques attempt to circumvent this problem by proposing multiple initial states, and then assume that the configuration having the lowest free energy represents the “correct” final phase.

The GL simulations considered here start from five different initial states, namely, random, cholesteric, and systems with O5, O2, O 8 , and O 8 + symmetries.55 To generate a cholesteric initial state, the z-component of eigenvector n is defined by

n z = cos ( 2 π z ) , sin ( 2 π z ) , 0 .
(26)

For O2 symmetry, the tensor order parameter is defined by

Q x x = A cos ( 2 q 0 z ) cos ( 2 q 0 y ) , Q x y = A sin ( 2 q 0 z ) ,
(27)

whereas for O5 by

Q x x = A { 2 cos ( q y ) cos ( q z ) cos ( q x ) cos ( q z ) cos ( q x ) cos ( q y ) } , Q x y = A { 2 cos ( q y ) sin ( q z ) 2 cos ( q x ) sin ( q z ) sin ( q x ) sin ( q y ) } ,
(28)

and for O8 by

Q x x = A { 2 cos ( q y ) sin ( q z ) + sin ( q x ) cos ( q z ) + cos ( q x ) sin ( q y ) } , Q x y = A { 2 cos ( q y ) cos ( q z ) 2 sin ( q x ) sin ( q z ) sin ( q x ) sin ( q y ) } ,
(29)

where q x = 2 q 0 x , q y = 2 q 0 y , and q z = 2 q 0 z . Constant A defines the amplitude and A > 0 for O2, O5, and O 8 + , whereas A < 0 for O 8 .

Figure 6 shows the free energy difference between GL and MC simulations for different chiralities. The GL calculations start from the initial conditions mentioned above, whereas the MC simulations always start from a random state. As one can appreciate in the figure, all GL minimizations get trapped; in some cases, the local minima are several kBT units above the relevant global minimum. As a general trend, we find that the deviations between the final phase free energies increase as the chirality increases. Here, we note that the simulated cell sizes were purposely chosen in a way that frustrates the system and prevents any known phase to be form. Two final MC states at N = 1.5 for L = 15.81ξN and N = 2.5 with L = 31.62ξN are included, where iso-surfaces represent values of S = 0.6. GL relaxations that started from the final MC state are also shown in the figure (red MC-GL). As discussed before, the GL remains the same state, but the free energy increases due to the smoothing of an averaged minimization. This issue may be diminished by increasing the numerical resolution in a way that smaller scales in gradients are captured. The problem here is of course that FD, being a collocation method, relies on a large number of nodes to resolve a fine gradient, thereby increasing computational demands.

FIG. 6.

Free energy difference for a chiral LC in the bulk with periods (a) L = 15.81ξN and (b) L = 31.62ξN. The chirality is increased according to the number of turns in the cell, N = 2L/pch. Monte Carlo relaxations, which are able to find the minimum value of the free energy, are used for the reference value in all cases. The final MC states are shown for N = 1.5 with L = 15.81ξN cell and for N = 2.5 with L = 31.62ξN cell. Iso-surfaces represent values of S = 0.6.

FIG. 6.

Free energy difference for a chiral LC in the bulk with periods (a) L = 15.81ξN and (b) L = 31.62ξN. The chirality is increased according to the number of turns in the cell, N = 2L/pch. Monte Carlo relaxations, which are able to find the minimum value of the free energy, are used for the reference value in all cases. The final MC states are shown for N = 1.5 with L = 15.81ξN cell and for N = 2.5 with L = 31.62ξN cell. Iso-surfaces represent values of S = 0.6.

Close modal

Figure 7 shows representative final states, with S = 0.6 iso-surfaces, from GL and MC minimizations. These phases are at N = 2.5 for L = 31.62ξN. All final states, predicted by GL, are different than the MC global minimum. This phase is an abnormal combination of a chiral and a blue phase due to the lack of symmetry. Therefore, any blue phase predicted by GL represents, in reality, local minima. For the GL that started from a random configuration, no rotation was found to match the final MC state, and we are confident that the phases are different.

FIG. 7.

Final states from GL and MC relaxations of a chiral LC with N = 2.5 for L = 31.62ξN. Ginzburg-Landau minimizations start from different initial configurations and become trapped in local minima. The MC simulations are able to find the global minimum. Iso-surfaces represent values of S = 0.6.

FIG. 7.

Final states from GL and MC relaxations of a chiral LC with N = 2.5 for L = 31.62ξN. Ginzburg-Landau minimizations start from different initial configurations and become trapped in local minima. The MC simulations are able to find the global minimum. Iso-surfaces represent values of S = 0.6.

Close modal

A Monte Carlo method has been proposed to identify the liquid crystal morphologies that correspond to the minimum free energy of widely used Landau-de Gennes continuum models. This new implementation represents a special case of a wider class of methods referred to as theoretically informed MC methods, which have been widely used in polymer physics. Our results are that, for simple nematic materials, the proposed approach is able to identify the same morphologies that are generated by traditional minimization techniques. For more complex materials, such as the blue phases of chiral liquid crystals, the proposed approach avoids metastable states and, for all cases considered in this work, is able to find the morphologies that correspond to the global free energy minimum. In this respect, the proposed MC method is superior to traditional GL minimization techniques, which get easily trapped in metastable morphologies. The proposed technique offers the additional advantage of being easy to implement, easily parallelizable, and able to handle composite systems, including suspensions of particles in liquid crystals.

The authors acknowledge support from the Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering Division, Biomaterials Program under Grant No. DE-SC0004025. Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program of the Argonne Leadership Computing Facility at Argonne National Laboratory. Additional development work was performed using the Argonne Laboratory Resource Computing Center (LCRC) and The University of Chicago Midway cluster. J.C.A.-P. is thankful to CONACYT for Postdoctoral Fellowships Nos. 186166 and 203840. J.P.H.-O. is grateful to funding provided by the Universidad Nacional de Colombia Ph.D. grant and COLCIENCIAS under the Contract No. 110-165-843-748, “Patrimonio Autónomo Fondo Nacional de Financiamiento para la Ciencia, Tecnología y la Innovación Francisco José de Caldas.”

For biaxial systems, there is no direct transformation (or at least linear) from the Frank-Oseen representation in terms of ∇n and ∇n′ to the ∇Q representation. However, the free energy density in Eq. (5) may be used for biaxial systems as long as it is understood that the elastic constants Li are not exactly as defined by Eq. (7); they should be a function of S, η, and k s. We use expressions in Eq. (7) for convenience and to simplify comparison to available experimental data. A complete description, however, would require the L constants to be measured directly, experimentally, or by molecular dynamics.64 

For an i-th collocation point in the bulk, the free energy density contribution fbulk is expressed in terms of the tensor order parameter and its gradients as follows:

f bulk ( Q i , Q i , Q j ) = f L ( Q i ) + f E ( Q i ) + j f E ( Q j ) ,
(B1)

where fL(Qi) and fE(∂Qi) are the Landau and elastic energy density contribution at the collocation point, respectively, and fE(∂Qi) are the contributions from the elastic density from the j neighbor to the collocation point i. For a typical second order FDM scheme, the neighbor elastic contribution includes six surrounding points. Similarly, for collocation points at the surface, the surface energy fsurface is computed as follows:

f surface ( Q i ) = f S ( Q i ) + j f E ( Q j ) ,
(B2)

where the surface term fS is calculated locally, whereas the elastic contribution must include the effects on the neighboring bulk point j.

It is important to point out that second order central finite-differences are not desired because, as the tensor is changed at the collocation point, the contributions at the neighboring nodes cancel each other. Therefore, higher order non-central schemes are used to approximate the gradients. For instance, as the point approaches the surface, a 3rd order scheme is used, resulting in Ref. 65,

( Q i ) 1 6 h ( 2 Q i + 1 + 3 Q i 6 Q i 1 + Q i 2 ) ,
(B3)

where h the distance between points in the mesh.

Numerical integration for the free energy density and the surface free energy is carried out directly accounting to the uniform mesh used in the FDM scheme, i.e.,

d 3 x f ( Q , Q ) i , j , k f ( Q i , j , k , Q i , j , k ) Δ V
(B4)

and

d 2 x f ( Q , Q ) i , j f ( Q i , j , Q i , j ) Δ A ,
(B5)

where ΔV = V/bulk nodes and ΔA = A/surface nodes.

For the particular case of a spherical droplet, the FD mesh is obtained by intersecting a uniform mesh on a cubic box with a sphere of radius R ± δR; points with a radius less than RδR are designated as internal points. The mesh resolution was selected on the basis of the smallest length scale to be resolved. For example, a 1 μm cube with a ξN = 10 nm nematic coherence length required at least 201 nodes in each direction, resulting in a 2 × 106 node mesh.

We decided to compare the equilibrium states obtained by the new MC method with those from the Ginzburg-Landau minimization. In the latter method, in order to minimize free energy functional (3), a projector operator ΠQ is defined according to the natural construction of the tensor Q, which is symmetric and traceless,

Π Q B = 1 2 B + B T 1 3 tr B δ ,
(C1)

for any tensor B x 3 × 3 , for any bulk point xV, and surface point xσS. The free energy must satisfy the following Euler-Lagrange equations:

Π Q δ F δ Q = 0 ,
(C2)

for xV and

Π Q δ F δ Q ν = 0 ,
(C3)

as boundary conditions (i.e., for xσS), where ν is the unit outward normal to the surface. The Volterra derivatives are defined by Ref. 9,

δ F δ Q = F Q x F Q .
(C4)

The solution to these equations is found by allowing the tensor order parameter Q to evolve towards equilibrium according to a Ginzburg-Landau relaxation equation of the form

Q t = 1 γ Π Q δ F δ Q ,
(C5)

with boundary conditions given in (C3), and where γ is a rotational viscosity (or diffusion) coefficient.48,66,67 Here, the effects of fluid flow are neglected; they are not expected to influence the equilibrium configurations.

The Ginzburg-Landau relaxation is carried out numerically using a semi-implicit Euler scheme for the time derivative, whereas for the spatial derivatives, a finite difference approach is employed.65,68 For the GL implementation, we may express the Volterra derivatives in terms of the independent components of the tensor order parameter Qij or following the orthonormal basis (the a form, Eq. (16)). Both formulations are equivalent and lead to the same free energy minima.

1.
N. A.
Lockwood
,
J. C.
Mohr
,
L.
Ji
,
C. J.
Murphy
,
S.
Palecek
,
J.
de Pablo
, and
N.
Abbott
,
Adv. Funct. Mater.
16
,
618
(
2006
).
2.
I.-H.
Lin
,
D. S.
Miller
,
P. J.
Bertics
,
C. J.
Murphy
,
J. J. D.
Pablo
, and
N.
Abbott
,
Science
332
,
1297
(
2011
).
3.
S.
Sivakumar
,
J. K.
Gupta
,
N.
Abbott
, and
F.
Caruso
,
Chem. Mater.
20
,
2063
(
2008
).
4.
J.
Gupta
,
S.
Sivakumar
,
F.
Caruso
, and
N.
Abbott
,
Angew. Chem., Int. Ed.
48
,
1652
(
2009
).
5.
U.
Manna
,
Y. M.
Zayas-Gonzalez
,
R. J.
Carlton
,
F.
Caruso
,
N. L.
Abbott
, and
D. M.
Lynn
,
Angew. Chem., Int. Ed. Engl.
52
,
14011
(
2013
).
6.
T.
Wang
,
J.
Zhuang
,
J.
Lynch
,
O.
Chen
,
Z.
Wang
,
X.
Wang
,
D.
LaMontagne
,
H.
Wu
,
Z.
Wang
, and
Y. C.
Cao
,
Science
338
,
358
(
2012
).
7.
J.
Moreno-Razo
,
E.
Sambriski
,
G.
Koenig
,
E.
Díaz-Herrera
,
N.
Abbott
, and
J.
de Pablo
,
Soft Matter
7
,
6828
(
2011
).
8.
J. A.
Moreno-Razo
,
E. J.
Sambriski
,
N. L.
Abbott
,
J. P.
Hernandez-Ortiz
, and
J. J.
de Pablo
,
Nature
485
,
86
(
2012
).
9.
E.
Kreyszig
,
Advanced Engineering Mathematics
, 8th ed. (
John Wiley & Sons
,
New York
,
1999
).
10.
S.
Grollau
,
O.
Guzman
,
N.
Abbott
, and
J. J.
de Pablo
,
J. Chem. Phys.
122
,
024703
(
2005
).
11.
A. K.
Bhattacharjee
,
G. I.
Menon
, and
R.
Adhikari
,
Phys. Rev. E
78
,
026707
(
2008
).
12.
V.
Tomar
,
T. F.
Roberts
,
N. L.
Abbott
,
J. P.
Hernandez-Ortiz
, and
J. J.
de Pablo
,
Langmuir
28
,
6124
(
2012
).
13.
V.
Tomar
,
S. I.
Hernandez
,
N. L.
Abbott
,
J. P.
Hernandez-Ortiz
, and
J. J.
de Pablo
,
Soft Matter
8
,
8679
(
2012
).
14.
M.
Ricci
,
M.
Mazzeo
,
R.
Berardi
,
P.
Pasini
, and
C.
Zannoni
,
Faraday Discuss.
144
,
171
(
2010
).
15.
F. A.
Detcheverry
,
D. Q.
Pike
,
P. F.
Nealey
,
M.
Müller
, and
J. J.
de Pablo
,
Phys. Rev. Lett.
102
,
197801
(
2009
).
16.
B. L.
Peters
,
A.
Ramírez-Hernández
,
D. Q.
Pike
,
M.
Müller
, and
J. J.
de Pablo
,
Macromolecules
45
,
8109
(
2012
).
17.
F. A.
Detcheverry
,
D. Q.
Pike
,
P. F.
Nealey
,
M.
Müller
, and
J. J.
de Pablo
,
Faraday Discuss.
144
,
111
(
2010
).
18.
J.
Hobdell
and
A.
Windle
,
Liq. Cryst.
23
,
157
(
1997
).
19.
P.
Gemünden
,
C.
Poelking
, and
K.
Kremer
,
Macromolecules
46
,
5762
(
2013
).
20.
T.
Gruhn
and
S.
Hess
,
Z. Naturforsch., A
51
,
693
(
1996
).
21.
R. W.
Ruhwandl
and
E. M.
Terentjev
,
Phys. Rev. E
56
,
5561
(
1997
).
22.
M.
Papoular
and
A.
Rapini
,
Solid State Commun.
7
,
1639
(
1969
).
23.
P. J.
Le Masurier
,
G. R.
Luckhurst
, and
G.
Saielli
,
Liq. Cryst.
28
,
769
(
2001
).
24.
G. R.
Luckhurst
and
G.
Saielli
,
Mol. Cryst. Liq. Cryst.
395
,
183
(
2003
).
25.
Z.
Zhi-Dong
,
C.
Chun-Rui
, and
M.
Dong-Lai
,
Chin. Phys. B
18
,
1560
(
2009
).
26.
C.
Chiccoli
,
P.
Pasini
,
A.
Šarlah
,
C.
Zannoni
, and
S.
Žumer
,
Phys. Rev. E
67
,
050703
(
2003
).
27.
H. G.
Galabova
,
N.
Kothekar
, and
D. W.
Allender
,
Liq. Cryst.
23
,
803
(
1997
).
28.
P.
Palffy-muhoray
,
E. C.
Gartland
, and
J. R.
Kelly
,
Liq. Cryst.
16
,
713
(
1994
).
29.
A.
Šarlah
and
S.
Žumer
,
Phys. Rev. E
60
,
1821
(
1999
).
30.
R. K.
Goyal
and
M. M.
Denn
,
Phys. Rev. E
75
,
021704
(
2007
).
31.
L.
Landau
and
E.
Lifshitz
,
Statistical Physics
, 3rd ed. (
Pergamon-Press
,
Oxford
,
1980
).
32.
J.
Fournier
and
P.
Galatola
,
Europhys. Lett.
72
,
403
(
2005
).
33.
O.
Guzmán
,
N.
Abbott
, and
J.
de Pablo
,
J. Polym. Sci., Part B
43
,
1033
(
2005
).
34.
O.
Guzman
,
N.
Abbott
, and
J. J.
de Pablo
,
J. Chem. Phys.
122
,
184711
(
2005
).
35.
A. N.
Beris
and
B. J.
Edwards
,
Thermodynamics of Flowing Systems
(
Oxford University Press
,
Oxford
,
1994
).
36.
P.
de Gennes
and
J.
Prost
,
The Physics of Liquid Crystals
(
Oxford University Press
,
Oxford
,
1995
).
37.
C.
Oseen
,
Trans. Faraday Soc.
29
,
883
(
1933
).
38.
C.
Oseen
,
Fortschr. Chem., Phys. Phys. Chem.
20
,
1
(
1929
).
39.
C.
Oseen
,
Ark. Mat., Astron. Fys. A
19A
,
16
(
1925
).
40.
Bernal
,
Lawrence
,
Oseen
,
H.
Zocher
,
Kast
,
Ornstein
,
Feachem
,
N. K.
Adam
,
Desch
,
Lowry
,
Rawlins
,
Herrmann
,
Foex
,
Malkin
,
Porter
,
Stewart
,
G. W.
Stewart
,
A. M.
Taylor
,
Ostwald
,
Carpenter
, and
G. S.
Hartley
,
Trans. Faraday Soc.
29
,
1060
(
1933
).
41.
F.
Frank
,
Discuss. Faraday Soc.
25
,
19
(
1958
).
42.
H.
Mori
,
E.
Gartland
,
J.
Kelly
, and
P.
Bos
,
Jpn. J. Appl. Phys.
38
,
135
(
1999
).
43.
T. C.
Lubensky
,
Phys. Rev. A
2
,
2497
(
1970
).
44.
J.
Ericksen
,
Phys. Fluids
9
,
1205
(
1966
).
45.
E. C.
Gartland
and
S.
Mkaddem
,
Phys. Rev. E
59
,
563
(
1999
).
46.
L.
Longa
,
D.
Monselesan
, and
H.-R.
Trebin
,
Liq. Cryst.
2
,
769
(
1987
).
47.
M.
Nobili
and
G.
Durand
,
Phys. Rev. A
46
,
R6174
(
1992
).
48.
M.
Ravnik
,
G. P.
Alexander
,
J. M.
Yeomans
, and
S.
Žumer
,
Faraday Discuss.
144
,
159
(
2010
).
49.
M.
Skarabot
,
M.
Ravnik
,
S.
Žumer
,
U.
Tkalec
,
I.
Poberaj
,
D.
Babic
,
N.
Osterman
, and
I.
Musevic
,
Phys. Rev. E
76
,
051406
(
2007
).
50.
F. R.
Hung
,
Phys. Rev. E
79
,
021705
(
2009
).
51.
A.
Saupe
,
J. Chem. Phys.
75
,
5118
(
1981
).
52.
P.
Kaiser
,
W.
Wiese
, and
S.
Hess
,
J. Non-Equilib. Thermodyn.
17
,
153
(
1992
).
53.
A. K.
Bhattacharjee
,
G. I.
Menon
, and
R.
Adhikari
,
J. Chem. Phys.
133
,
044112
(
2010
).
54.
D.
Frenkel
and
B.
Smith
,
Understanding Molecular Simulations-From Algorithms to Applications
(
Academic Press
,
San Diego
,
1996
).
55.
A.
Dupuis
,
D.
Marenduzzo
, and
J.
Yeomans
,
Phys. Rev. E
71
,
011703
(
2005
).
56.
I. W.
Stewart
,
The Static and Dynamic Continuum Theory of Liquid Crystals: A Mathematical Introduction
(
Taylor & Francis
,
London
,
2004
).
57.
G.
Barbero
and
R.
Barberi
,
J. Phys. (France)
44
,
609
(
1983
).
58.
J. K.
Gupta
,
J. S.
Zimmerman
,
J. J.
de Pablo
,
F.
Caruso
, and
N. L.
Abbott
,
Langmuir
25
,
9016
(
2009
).
59.
D.
Seč
,
T.
Porenta
,
M.
Ravnik
, and
S.
Žumer
,
Soft Matter
8
,
11982
(
2012
).
60.
S.
Kralj
and
S.
Žumer
,
Phys. Rev. A
45
,
2461
(
1992
).
61.
J. H.
Erdmann
,
S.
Žumer
, and
J. W.
Doane
,
Phys. Rev. Lett.
64
,
1907
(
1990
).
62.
A.
Londono-Hurtado
,
J. C.
Armas-Pérez
,
J. P.
Hernández-Ortíz
, and
J. J.
de Pablo
,
Soft Matter
11
,
5067
(
2015
).
63.
S. I.
Hernández
,
J. A.
Moreno-Razo
,
A.
Ramírez-Hernández
,
E.
Díaz-Herrera
,
J. P.
Hernández-Ortiz
, and
J. J.
de Pablo
,
Soft Matter
8
,
1443
(
2011
).
64.
A. A.
Joshi
,
J. K.
Whitmer
,
O.
Guzman
,
N. L.
Abbott
, and
J. J.
de Pablo
,
Soft Matter
10
,
882
(
2014
).
65.
T. A.
Osswald
and
J. P.
Hernández-Ortiz
,
Polymer Processing: Modeling and Simulation
(
Carl Hanser-Verlag
,
Munich
,
2006
).
66.
B. T.
Gettelfinger
,
J. A.
Moreno-Razo
,
G. M.
Koenig
, Jr.
,
J. P.
Hernández-Ortiz
,
N. L.
Abbott
, and
J. J.
de Pablo
,
Soft Matter
6
(
5
),
896
(
2010
).
67.
J. P.
Hernandez-Ortiz
,
B. T.
Gettelfinger
,
J.
Moreno-Razo
, and
J. J. D.
Pablo
,
J. Chem. Phys.
134
,
134905
(
2011
).
68.
M.
Ravnik
and
S.
Žumer
,
Liq. Cryst.
36
,
1201
(
2009
).