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.

## I. INTRODUCTION

Liquid crystals (LCs) are finding a number of emerging applications in areas ranging from biomaterials^{1} to sensors^{2} 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 Hess^{20} 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 Terentjev^{21} 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-Lai^{25} 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., *k*_{11}≠*k*_{33}). For spherical geometries (e.g., droplets), Goyal and Denn presented a particularly important study in which they examined the effects of homeotropic anchoring^{30} on morphology; to do so, these authors extended the Gruhn-Hess formalism to include a saddle-splay term (*k*_{24}), 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 slit^{14} (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 energy^{31} is supplemented by a Rapini-Papoular surface contribution in the case of homeotropic anchoring^{22} 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.

## II. THERMODYNAMIC DESCRIPTION

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*),

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,

where $ \beta \u2212 1 = k B T \u02c6 $ (*k _{B}* is Boltzmann constant and $ T \u02c6 $ 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,

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

where parameter *A* sets an energy density scale for the model. A dimensionless parameter *U* determines the IN transition^{10,33,34} and the magnitude of the average scalar order parameter. For Eq. (4), we adopt Doi’s notation^{35,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

where *L*_{1} and *L*_{5} 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-Oseen

^{37–41}elastic constants

^{42}by the following expressions:

where *S* is the scalar order parameter and *k*_{11}, *k*_{22}, and *k*_{33} are non-vanishing elastic moduli corresponding to independent splay, twist, and bend modes.^{43} Here, *L*_{5} quantifies the chirality of the LC and the pitch *p*_{ch} = 2*π*/*q*_{0}. Without loss of generality, in this work, we adopt a single elastic constant approximation, i.e., *k*_{11} = *k*_{22} = *k*_{33}. We satisfy the Ericksen inequalities^{44} and the stability ranges derived by Gartland and Mkaddem^{45} and Longa *et al.*^{46}

For homeotropic anchoring, we adopt a Rapini-Papoular^{22,47–50} expression of the form

where *W _{R}* represents the strength of the homeotropic anchoring and

**Q**

_{⊥}is the perpendicular tensor orientation preferred at the surface, defined as follows:

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 expression

^{32}

where $ Q \xaf =Q+S\delta /3$, $ Q \xaf \u22a5 =p\u22c5 Q \xaf \u22c5p$ 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

*ξ*; they are defined as follows:

_{S} 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

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

## III. UNIFORM MONTE CARLO MINIMIZATION ALGORITHM

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 *P*_{prop}(*o* → *n*). Trial transitions are accepted with probability,

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 to^{43,51}

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

and

The eigenvectors, **n** and **n**′, corresponding to the maximum and second largest eigenvalues, define an orthonormal basis $ n , n \u2032 , n \xd7 n \u2032 $ 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-workers^{52} and, more recently, by Bhattacharjee *et al.*,^{11,53}

The orthonormal basis is defined by five tensors,

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,

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

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**,

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

In addition, since the characteristic polynomial of **Q** is

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

where component *μ* is selected at random from 1 to 5, *ζ* is a random number distributed uniformly between [0, 1], and $ \delta \mu \u0304 $ 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 $ \delta \mu \u0304 $ 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 \u02c6 $.^{54} The simulation starts from a high non-dimensional effective temperature $T=A \xi N 3 /\beta $ 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 < *β* < 10^{12} 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.

## IV. RESULTS

The parameters adopted in all nematic simulations are chosen to be representative of 4′-pentyl-4-cyanobiphenyl: *A* ≈ 1 × 10^{5} J/m^{3} and *k*_{11} = *k*_{22} = *k*_{33} ≈ 6 pN, resulting in *ξ _{N}* = 7.15 nm. For chiral systems, we use $ \xi N = 40 nm$ and

*L*

_{1}= 25 pN. In the Landau-de Gennes framework, parameter

*U*sets a value for the nematic scalar order parameter,

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}

### A. Hybrid nematic cell

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 *W _{R}* =

*W*=

_{F}*W*. The results presented here encompass strengths from weak (

*W*= 10

^{−7}J/m

^{2}) to strong (

*W*= 10

^{−2}J/m

^{2}) 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.

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

*ξ*, 3.0

_{N}*ξ*, and 4.0

_{N}*ξ*. The figure includes the scalar order parameter

_{N}*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

*ξ*, the value of

_{N}*S*is always higher than

*η*and

*φ*varies linearly from top to bottom; there is no UB transition. At

*d*= 3.0

*ξ*, the system shows an increment in the biaxiality, and

_{N}*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

*ξ*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

_{N}*φ*: red corresponds to orientations along the bottom wall, and blue corresponds to orientations perpendicular to the bottom wall.

### B. Spherical nano-droplet

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.

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 *L*_{1} 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.

### C. Chiral LC

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 $ \xi N = 40 nm$ and *L*_{1} = 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

*p*

_{ch}= 56.6

*ξ*(

_{N}*N*= 2

*L*/

*p*

_{ch}= 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 \u2212 $ symmetry.

^{55}

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 10^{5} 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 × 10^{5} steps. In contrast, for the GL minimization, the eigenvalues (*S* and *η*) approach their equilibrium values early in the minimization (after approximately 3 × 10^{5} 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 × 10^{5} time steps). These differences help explain why GL gets more easily trapped in metastable states.

We next proceed to consider different cell sizes, namely, *L* = 15.81*ξ _{N}* and

*L*= 31.62

*ξ*, 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.

_{N}The GL simulations considered here start from five different initial states, namely, random, cholesteric, and systems with *O*_{5}, *O*_{2}, $ O 8 \u2212 $, and $ O 8 + $ symmetries.^{55} To generate a cholesteric initial state, the *z*-component of eigenvector **n** is defined by

For *O*_{2} symmetry, the tensor order parameter is defined by

whereas for *O*_{5} by

and for *O*_{8} by

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 *O*_{2}, *O*_{5}, and $ O 8 + $, whereas *A* < 0 for $ O 8 \u2212 $.

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 *k _{B}T* 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

*ξ*and

_{N}*N*= 2.5 with

*L*= 31.62

*ξ*are included, where iso-surfaces represent values of

_{N}*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.

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.

## V. CONCLUSIONS

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.

## Acknowledgments

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.”

### APPENDIX A: ABOUT THE *L* ELASTIC CONSTANTS

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 *L _{i}* 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}

### APPENDIX B: FINITE DIFFERENCES SCHEME

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

where *f _{L}*(

**Q**

_{i}) and

*f*(∂

_{E}**Q**

_{i}) are the Landau and elastic energy density contribution at the collocation point, respectively, and

*f*(∂

_{E}**Q**

_{i}) 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

*f*

_{surface}is computed as follows:

where the surface term *f _{S}* 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,

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

and

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 × 10

^{6}node mesh.

### APPENDIX C: GINZBURG-LANDAU RELAXATION

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,

for any tensor $B x \u2208 \u211c 3 \xd7 \u211c 3 $, for any bulk point **x** ∈ *V*, and surface point **x** ∈ *σ _{S}*. The free energy must satisfy the following Euler-Lagrange equations:

for **x** ∈ *V* and

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,**

*ν* 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

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 *Q _{ij}* or following the orthonormal basis (the

**a**form, Eq. (16)). Both formulations are equivalent and lead to the same free energy minima.