A theoretically informed Monte Carlo method is proposed for Monte Carlo simulation of liquid crystals on the basis of theoretical representations in terms of coarse-grained free energy functionals. The free energy functional is described in the framework of the Landau-de Gennes formalism. A piecewise finite element discretization is used to approximate the alignment field, thereby providing an excellent geometrical representation of curved interfaces and accurate integration of the free energy. The method is suitable for situations where the free energy functional includes highly non-linear terms, including chirality or high-order deformation modes. The validity of the method is established by comparing the results of Monte Carlo simulations to traditional Ginzburg-Landau minimizations of the free energy using a finite difference scheme, and its usefulness is demonstrated in the context of simulations of chiral liquid crystal droplets with and without nanoparticle inclusions.

## I. INTRODUCTION

Theoretical and computational descriptions of liquid crystalline systems are of interest in a wide variety of applications, ranging from development of biomaterials^{1} to design of sensors^{2} and optical devices.^{3–5}

Liquid crystals (LCs) exhibit intriguing structural features over multiple length scales. Defects, for example, have characteristic dimensions in the range of several nanometers, whereas elastic interactions and interfacial effects can be felt over micrometer length scales. Accordingly, a wide range of computational descriptions have been used to model LCs, including atomistic simulations,^{6–10} coarse-grained, particle based simulations,^{11–18} and continuum simulations.^{19–25}

One of the most widely used models of liquid crystal structure is that originally proposed by Landau and de Gennes,^{26–29} where the thermodynamic state of the system is described at the level of a local tensorial order parameter. In this theory, a coarse-grained free energy functional is expressed in terms of various gradients of the tensor order parameter, each corresponding to different types of deformation which are penalized according to several elastic coefficients. Such coefficients can be measured experimentally or, importantly, can be determined directly from many-particle atomistic or coarse-grained simulations,^{14} thereby providing a systematic means for describing nematic liquid crystals at multiple length scales. The minimization of free energy can be carried out through a variety of techniques, some of which involve solution of the corresponding Euler-Lagrange equations or dynamic relaxations that introduce an explicit time dependence.^{30–40} Monte Carlo (MC) simulations have also been applied to study liquid crystals.^{41–46} Recently, a theoretically informed MC approach^{47} was proposed in the particular context of the Landau-de Gennes tensorial representation, using a finite-difference (FD) scheme to integrate the underlying equations. In this work, that approach is extended to introduce a finite-element (FE) description, which provides a better description of complex nematic fields at a fraction of the computational demands incurred by traditional finite-difference techniques.

For numerical simulations, the free energy functional and its Volterra derivatives must be approximated over the domain (and its corresponding boundaries) of interest. As noted above, the most widely used scheme to do so is FDs. When a simple time integration method is applied, it results in a matrix-free algorithm.^{48–52} Dirac delta functions are used as weighting functions in the weighted residuals scheme, leading to minimization of errors over a discrete set of points or nodes. By construction, FD schemes rely on a large number of nodes to provide a good approximation of the underlying free energy. Description of complex geometries also requires a large number of nodes.^{53} This latter issue can be alleviated by resorting to radial basis function (RBF) methods, where the field approximation is implemented through a linear combination of radial basis functions, resulting in a mesh-less method. Unfortunately, the resulting algorithm includes full matrices, which restrict the applicability of the method to small systems.^{22,23,32,33}

The MC approach proposed here adopts 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.^{47} It offers the advantage that one can sample the alignment tensor parameters independently, without a need for cumbersome constraints, and the scalar order parameters, uniaxiality *S* and biaxiality *η*, as well as their corresponding director vectors, **n** and **n**′, are sampled uniformly. The integration of the free energy functional is carried out with an efficient Gaussian quadrature that is applied piecewise in discrete FEs. Common meshing techniques and algorithms can therefore be used, thereby enabling their applicability to arbitrary geometries. Also note that the FE integration scheme allows for coarser meshes than those required in finite differences, thereby permitting simulation of large systems.^{54}

The validity of the proposed method, referred to as FE-MC, is examined in the context of droplets of chiral LCs, where its results are compared to those of finite-difference approaches. We also present results from simulations of solid tetrahedral inclusions embedded in LC droplets, as a way to illustrate the method’s applicability to complex geometries.

## II. THERMODYNAMIC DESCRIPTION

As mentioned above, the thermodynamic state of the LC is described in terms of the local alignment tensor **Q**(**x**), defined in terms of the second moment 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. The tensor order parameter in Eq. (1) includes all relevant information regarding the ensemble average alignment of LC molecules, i.e., directions and level of order, as indicated by its eigenvalues and eigenvectors. Therefore, the alignment tensor

**Q**(

**x**) may be written as

^{55,56}

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 \u2212S3$. The eigenvectors, **n** and **n**′, corresponding to the maximum and second largest eigenvalues, respectively, define an orthonormal basis $n,n\u2032,n\xd7n\u2032$ for LC orientation. The scalar order parameters, *S* and *η*, are bounded by

Note that the alignment tensor is symmetric and traceless and out of its nine Cartesian components, only five are actual degrees of freedom. Thus, it can be expressed in terms of an orthonormal tensor basis^{31,57,58} of the form

where the orthonormal basis is defined by five tensors,

Here, $x\u02c6$, $y\u02c6$, and $z\u02c6$ are the canonical ℜ^{3} bases, $AST$ is the symmetric-traceless projection operation, and *δ _{ij}* is the Kronecker delta. Because the $Tm$ basis is orthonormal,

ensures that the five scalar components *a*_{j} of the alignment tensor are simple projections

and that they provide a unique and independent way of specifying the configuration of the alignment tensor field.^{47,48,54} This feature is helpful during free energy relaxation of the free energy and for including fluctuations in that relaxation.

The free energy functional consists of the sum of three terms: a short-range Landau polynomial expansion of the tensor invariants, to describe the isotropic-nematic (IN) transition, a long-range free energy density that penalizes elastic distortions, and a surface free energy contribution, given by

Note that our proposed simulation approach is not restricted to these three contributions. Additional terms, including as electromagnetic or thermal contributions, may also be included.

The Landau free energy density *f _{L}* is a phenomenological approximation obtained as a truncated expansion with respect to the invariants of

**Q**,

where parameter *A* sets an energy density scale for the model. A dimensionless parameter *U* determines the IN transition.^{30,59,60} For Eq. (9), we adopt Doi’s notation and consider a unique value for the energy scale *A*.^{26,61}

The elastic contribution is given in terms of gradients of the alignment tensor according to^{28,29,62}

where *L*_{s} are elastic constants and *ϵ _{ijk}* is the Levi-Civita pseudotensor. For a

*uniaxial*system, the elastic constants in Eq. (10) are related to the first order Frank-Oseen

^{62–67}elastic constants through

where *S* is the scalar order parameter and *k*_{11}, *k*_{22}, *k*_{33}, and *k*_{24} are non-vanishing elastic moduli corresponding to the independent splay, twist, bend, and saddle-splay modes.^{55} *L*_{5} describes the chirality of the system through the pitch *p*_{ch} = 2*π*/*q*_{0}.

The surface contributions to the free energy describe the interaction of the liquid crystal with boundaries or interfaces. For homeotropic anchoring, a Rapini-Papoular^{37,68–71} energy density is used

where *W* represents the strength of the anchoring, $Q\u22a5=S[\nu \nu \u221213\delta ]$ is the perpendicular tensor preferred at the surface, and ** ν**(

**x**) is the normal unit vector at any point

**x**on the surface. For degenerate planar anchoring, a 4th order Fournier-Galatola energy density is adopted

^{72}

where $Q\xaf=Q+S\delta /3$, $Q\xaf\u22a5=p\u22c5Q\xaf\u22c5p$ is the tensor’s projection on the surface, and **p** = ** δ** −

**.**

*νν*The physical properties of the LC serve to define two characteristic length scales: the nematic coherence length *ξ _{N}* and the surface extrapolation length

*ξ*. These are given by

_{S} The latter scale is also known as the Kleman-de Gennes extrapolation distance.^{26} 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 can thus be expressed in dimensionless form as

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

## III. FINITE ELEMENT INTEGRATION

A piecewise polynomial approximation scheme is used to perform the numerical integration of the free energy functional. In the finite element discretization adopted here, iso-parametric elements are used to represent the domain of interest. The free energy functional in Eq. (8) is then calculated as follows:

where *N _{E}* and

*N*are the number of discrete elements at the volume and the surface, respectively. The element alignment tensor (

_{ES}**a**

_{ı}) and its gradients (∇

**a**

_{ı}) are approximated using FE shape interpolation functions,

*M*

_{α}, according to

where *α* = 1, …, *N _{N}* and where

*N*is the total number of nodes per element. These expressions result in approximate values for the Landau ($f\u02c6L,\u0131$), elastic ($f\u02c6E,\u0131$), and surface ($f\u02c6S,\u0131$) free energy densities per element ı.

_{N}The element integrals in Eq. (20) are obtained with a Gaussian quadrature over the isoparametric elements. The resulting numerical integration, for tetrahedral elements, is given by^{53}

for the volume integrals and

for the surface integration. Here, ** ζ** represents the isoparametric coordinate system and

**J**and

**J**

_{S}are the Jacobian matrices for the volume and surface isoparametric transformations, respectively.

^{53,73}

In this manuscript, linear, six-noded, isoparametric tetrahedral elements are adopted, which provide an excellent level of approximation for the deformations considered here. Cubit (ver 14.1),^{74} used under the Argonne National Laboratory’s license, is used to generate the tetrahedral meshes. LibMesh^{75} libraries are then used to optimize and parametrize the original Cubit meshes. The FE scheme, combined with the Cubit and LibMesh libraries, provides a platform that can be easily modified to include any type of geometry, additional free energy contributions, and higher order approximations (element types). Note that in addition to the typical connectivity matrix, which lists the nodes at each element, a node-influence matrix is constructed. In this matrix, the elements that are influenced by a specific node are listed—information that we rely on for efficient *O*(1) MC algorithms.

## IV. THERMAL ANNEALING DURING FREE ENERGY RELAXATION

An annealing process can be introduced to avoid getting trapped in local minima along the underlying free energy landscape. To do so, one can introduce a fluctuating temperature $T\u02c6$, not to be confused with the thermodynamic temperature. The rationale behind annealing comes from the fact that during traditional free energy relaxations of liquid crystals, scalar order parameters, *S* and *η*, are relaxed first, followed by the “accommodation” of the corresponding eigenvectors. Including thermal fluctuations during the relaxation, following annealing, forces an **n**-to-*S* relaxation, which facilitates global minimization. In this paper, we use a previously reported theoretically informed MC relaxation that uses Gaussian integration over FE discretizations.

The idea behind the MC is to generate a Markov chain of states using Metropolis sampling of the free energy functional. At each trial configuration of the alignment tensor field **a**(**x**), a probability is assigned proportional to the Boltzmann factor

where $\beta \u22121=kBT\u02c6$ and *k _{B}* is the Boltzmann constant. The MC minimization of the free energy starts from a random distribution of the alignment tensor field in the bulk and surfaces. The field is approximated numerically on the discrete FE mesh. The free energy functional or Hamiltonian,

*F*=

*F*(

**a**, ∇

**a**), is then obtained by the FE quadrature integration. A node is selected randomly, in the bulk or surface, and a MC trial move is attempted by randomly displacing

**a**(

**x**) according to

^{47}

where “n” indicates the new, trial configuration, “o” is the original one, 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*_{μ}.

Note that a trial move only involves a local free energy difference from those elements that are influenced by the selected node. This is the place where the node-influence matrix is used to obtain an efficient *O*(1) method. The free energy minimization follows an annealing process that is controlled through the annealing temperature $T\u02c6$.^{47,76} The simulation starts from a high non-dimensional effective temperature $T=A\xi N3\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 < 1/*T* < 10^{12} over the total number of MC steps. A step is defined by the average number of accepted moves per FE nodes; we use between 200 000 and 500 000 accepted MC moves per node.

## V. RESULTS

The parameters adopted in all nematic simulations are chosen in a way that $\xi N=L1/A=10$ nm, *k*_{11} = 16 pN, and *W* = 1 × 10^{−3} J/m^{2} (moderate anchoring strength). We use the single constant approximation, i.e., *k*_{11} = *k*_{22} = *k*_{33}, and set *k*_{24} = 0. These conditions are chosen for simplicity; the FE-MC method does not require these common simplifications. Finally, *L*_{5} is changed accordingly as the chirality is increased. In the Landau-de Gennes framework, parameter *U* sets a value for the nematic scalar order parameter,

Bulk systems are isotropic phase for values of *U* < 2.8 and nematic for *U* ≥ 2.8. The results in this work are shown as a function of *U* for 2.7, 3, 3.5, and 4.0, which lead to *S*_{Bulk} = 0.333, 0.5, 0.6159, and 0.683, respectively.

To examine the applicability of the proposed method, we first consider the relaxation of chiral LC droplets of radius *R* = 250 nm. For the purposes of comparison, in addition to the free energy minimization through the proposed FE-MC, a traditional Ginzburg-Landau (GL) relaxation is also applied, using FDs for the domain discretization and field approximations. Regular meshes are used for the FD-GL simulations; additional details about this method can be found in the literature.^{47–49} Homeotropic (Rapini-Papoular) and degenerate planar (Fournier-Galatola) anchorings are imposed. The relaxations, by the two algorithms, are compared as a function of the isotropic-nematic parameter *U* and the LC chirality, which is expressed in terms of *N* = 4*R*/*p*_{ch}. Both methods converge relatively fast when homeotropic and low-chiral droplets are simulated at high *U*. In contrast, planar anchored droplets with high chirality and low *U* require longer simulations, as discussed below.

Figure 1 shows the Central Processing Unit (CPU) time (where CPU time is defined in hours) associated with typical FE-MC and FD-GL relaxations as a function of the number of nodes. These results are for planar, high chirality droplets with *N* = 4 at a *U* = 2.7. For this system, free energy minimization by traditional techniques is particularly demanding. FD-GL minimization involves over 5 × 10^{6} time steps. The FE-MC minimization involves 3 × 10^{6} MC cycles. A fully explicit Euler time integration is used in the FD-GL minimization that, as mentioned earlier, leads to an efficient *matrix-free* evolution scheme. A first glance at Fig. 1 would appear to suggest that the FD-GL scheme is considerably faster than the FE-MC method. However, we need to emphasize that, in order to achieve a comparable representation of the director field, the FD scheme requires a much larger number of nodes than the FE scheme due to its pointwise weighted residual character.^{53,77,78} Note that, for example, a coarse FD discretization, which for this specific case is below 69 599 nodes, is unable to converge, thereby imposing a resolution limit for the FD-GL method (dotted line in Fig. 1). This FD-GL limit corresponds to a distance between nodes equal to *ξ _{N}*. Finer discretizations, with 198 893 and 540 113 nodes, result in mesh characteristic distances of 0.7

*ξ*and 0.5

_{N}*ξ*, respectively. On the other hand, the FE-MC is able to converge to physically meaningful phases at low discretization levels. Recall that the FE Gaussian quadrature (with four Gauss points) provides an excellent approximation to the free energy functional (integration), requiring coarser discretizations. For instance, 33 679 nodes (171 931 elements) correspond to an average distance between nodes of 2

_{N}*ξ*. Finer discretizations for 1.5

_{N}*ξ*and

_{N}*ξ*require 70 363 nodes (368 063 elements) and 214 237 nodes (1 165 753 elements), respectively. Also note that MC relaxations that are started from random LC configurations converge to the correct solution, whereas FD-GL needs to be started from educated guesses for the initial configurations in order to avoid getting trapped in local minima.

_{N}^{47,48}A quantitative comparison between the minima identified by the two methods is shown in Figure 2, which summarizes the relative error (or relative difference) between FD-GL (198 893 nodes) and FE-MC (33 679 nodes) simulations. The FD-GL values were used as the reference value for the relative error, which is plotted as a function of the isotropic-nematic parameter

*U*and the LC chirality

*N*. The difference between both methods is always lower than 7% for all simulated conditions. Higher differences are obtained for high chirality and low

*U*, which is a consequence of the fact that, as the chirality is increased, finer resolutions are required for both methods. In addition, low

*U*values give rise to complex blue phases that also require finer discretizations. We also performed comparisons (not shown) between FD-GL with 540 113 nodes and FE-MC with 70 363 and 214 237 nodes. As the number of degrees of freedom increases, both methods converge progressively to the same limiting result; the relative differences between the FD-GL and the FE-MC decrease to a maximum of 3% and 1%.

As mentioned above, contrary to the FD-GL approach, the proposed FE-MC scheme is able to find physically meaningful results for coarse discretizations. These “coarse” final phases lack the LC orientational details that are resolved in finer discretizations; however, they are computationally inexpensive and provide the proper relaxation “learning” paths for more demanding simulations. Figure 3 presents a convergence study of the FE-MC scheme, where the relative error (or relative difference) between the free energy functional with 70 363 nodes and coarser descriptions is shown. The error is shown as a function of *U* and *N*. At low chirality, the differences between fine and coarse discretizations are lower than 3% for every *U*. As the chirality increases, the relative error increases, indicating the need for high-resolution descriptions in order to capture the chiral characteristic length scale (or pitch). In addition, at high chirality, there are significant differences in the errors between low and high *U*. Recall that blue phases appear at low *U*, which require finer discretizations. In Fig. 3, the convergence of FD-GL simulations with *N* = 2 is included for completeness. A FD-GL simulation with 540 113 nodes is used as the reference value. Notice that the maximum differences between FD discretizations are 3%; these relative differences are of the same order of magnitude as FD-GL vs. FE-MC, where the number of nodes in the FE-MC is one order of magnitude lower than the FD-GL. These observations serve to underscore an important advantage to the FE-MC, given that larger systems can be explored (from the perspective of CPU memory requirements). In addition, inspection of Fig. 1 also indicates that CPU times for the FE-MC scheme with ten times less nodes than the FD-GL are almost equivalent.

Figures 4 and 5 present the final, equilibrium states for homeotropic and planar LC droplets, respectively. These states are obtained after FE-MC and FD-GL relaxations for *R* = 250 nm, *U* = 3 and a moderate anchoring strength *W* = 1 × 10^{−3} J/m^{2}. Aside from the small differences in the free energy functional between the FE-MC and the FD-GL (Fig. 2), the final phases are completely equivalent. The only visible difference is the roughness that results from the FD scheme, because it uses a regular mesh to represent a spherical domain. This point is one of the key advantages provided by the FE integration method: complex geometries are described accurately and the approximations for volume and surface integrals are excellent. In the figures, blue iso-surfaces for the scalar order parameter *S* = 0.4 are shown. These iso-surfaces represent the location of LC defects. In addition, surface contours are representative of LC orientations at the droplet surface: red is for a complete perpendicular orientation, while blue is for a planar orientation. These orientations are measured with the dot product between the LC principal eigenvector **n** and the surface normal direction ** ν**. In Fig. 5, a cross section of the LC director field

**n**is included; the color scheme follows the same nomenclature as the surface contours.

For both anchoring orientations, cholesteric phases are observed at low chirality for *U* = 3. As the chirality is increased, blue phases start to appear. For homeotropic droplets (Fig. 4), an interesting transition is observed at *N* = 2, where both methods are trapped in a local cholesteric minimum. The FD-GL relaxation relies on a “good” initial guess to escape this local minimum and find the global blue-phase minimum. The FE-MC, on the other hand, passes through

the cholesteric phase, which is sustained during several MC cycles, and then reaches the final blue phase. For the planar case (Fig. 4), the chiral transition is observed at *N* = 2.5 where, even though both methods predict the same surface orientations, they predict slightly different blue phases. We believe that the FE-MC final phase is the correct one, given its ability to avoid local-minimum traps.

We now consider the applicability of the FE-MC method in a situation involving a complex geometry. A regular tetrahedral particle with an edge of *L* = 200 nm (height of 163 nm) is embedded in a spherical planar droplet. The anchoring strength for both, droplet and particle, is *W* = 1 × 10^{−3} J/m^{2}. Homeotropic and degenerate planar particles are considered with *U* = 3.165 (*S* = 0.548). Figure 6 shows two representative meshes for a particle at the center (*d*/*D* = 0) and close to the droplet surface (*d*/*D* = 0.8). They are generated using Cubit.^{74} The mesh resolution and node distribution were selected using the same criteria that are used in traditional computational fluid mechanics calculations. For instance, the mesh was setup by ensuring that the surface to volume ratio is between 0.2 and 0.25 and by ensuring that the smallest features (e.g., defect cores) be resolved. Similarly, the mesh resolution is chosen to be highly close to surfaces (for drops and particles), thereby providing a good description of any boundary and its vicinity. In addition, LibMesh^{75} tools were used to optimize the mesh conditioning. Detailed consistency and convergence tests were performed. During every free energy minimization, the value of the energy minima and the final LC orientation were followed as a function of the number of volumetric and superficial elements. In this way, proper discretizations were selected from convergence plots similar to those in Fig. 3. For situations involving high chirality or complex geometries, we found that these “convergence” rules were not generalizable and each system required an individual analysis. Therefore, simulations were performed from coarse to fine resolutions, including multiple surface/volume ratios; once a desired level of convergence was obtained, the full scale simulation was performed. We found that low discretizations were also helpful to select the “annealing” process for each system, which motivated even further the coarse-to-fine process.

Simulations are carried out in which the particle is moved along one of its principal axes parallel to the droplet center line. The particle location is measured by the relative position *d*/*D*, where *D* is the droplet diameter *D* = 2*R* = 1 *μ*m and *d* is the location of the tetrahedron center of mass (see Fig. 6). In this study, the particle was never allowed to touch the droplet surface, thus *d*/*D* ∈ [ − 0.8, 0.8]. Figure 7 shows a potential of mean force *ω*, defined as the free energy difference between the minimum free energy state (*d*/*D* = − 0.8 for both, homeotropic and planar, particles) and other locations,

Recall that a particle-free planar droplet undergoes a bipolar phase characterized by the presence of two boojum defects at opposite poles.^{11,33,49,79} It is then not surprising that the particle prefers to be close to the droplet surface where its induced defect merges with the boojums, thereby minimizing the elastic distortions in the LC. In this free energy minimum, the particle face is parallel to the droplet surface for both anchoring alignments. Figure 8 illustrates zoomed LC configurations surrounding the particle for the marked sections in Fig. 7. Sections “aa” and “dd” depict the lower free energy states for the homeotropic and planar particles. The magenta iso-surfaces represent the defect locations of low scalar order parameter *S* = 0.4.

As the particle moves towards the center of the droplet (Fig. 7), there are additional free energy penalizations due to the extra defects around it. These are ∼700 k_{B}T and ∼1600 k_{B}T for the homeotropic and planar particles, respectively (measured at *d*/*D* = 0, sections “bb” and “ee”). Interestingly, the additional energy penalizations are higher for the planar particle than for the homeotropic one; however, note that for the planar particle, there is an entire face where surface penalizations are large, whereas for the homeotropic particle, only the edges present defects. The other striking result is the different shape that the potential of mean force has between the two anchoring orientations. The planar particle shows a maximum in the free energy difference at *d*/*D* = 0 and then decreases monotonically to a local minimum at *d*/*D* = 0.8 (where the particle point is closer to the droplet surface, see section “ff”). The homeotropic particle potential of mean force continues to increase up to *d*/*D* = 0.6. This additional increase in the free energy difference is generated by the fact that the particle “forces” the boojum to be displaced (section “cc”), thereby inducing additional elastic distortions in the LC. At *d*/*D* > 0.6, the boojum is arrested by the particle and *ω* decreases strongly until the local minimum at *d*/*D* = 0.8 is reached (section “ff”). This phenomenon has also been observed for spherical nano-particles embedded in LC droplets.^{15,48} According to the potential of mean force in Fig. 7, *d*/*D* = 0.8 corresponds to a metastable state, i.e., there is a large energy penalization to move (or rotate) the particle to find the most stable state at *d*/*D* = − 0.8.

## VI. CONCLUSION

A theoretically informed 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. The free energy functional is integrated following a finite element piecewise discretization and a Gaussian quadrature. Such a combination of numerical techniques results in an efficient method applicable to complex geometries and large systems. It 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 validation of the method was performed through one-on-one comparisons to results from a traditional finite-difference based on Ginzburg-Landau relaxation technique in the context of several chiral LC droplets. It was found that the FE-MC method required approximately ten times fewer nodes to provide satisfactory approximations of the systems considered here; both methods converged to unique energy states as the discretization level was increased. An important advantage from the proposed MC approach is related to its capacity to reach such global minimum states without resorting to educated initial guesses of the equilibrium result.

To illustrate the applicability of the method, a solid tetrahedral particle was embedded in a LC planar droplet. The free energy of the system was calculated as a function of the particle position within the droplet. Homeotropic and planar tetrahedrons were considered. Similar to spherical nanoparticles, the solids were found to prefer locations near the surface of the droplets. The preferred tetrahedral orientation has one of the faces parallel to the droplet surface. A planar tetrahedron generates more expansive defects than homeotropic ones, because an entire face generates free energy penalizations.

## Acknowledgments

The study of nanoparticles suspended in liquid crystals reported in this work is supported by the National Science Foundation through Grant No. DMR-1410674. The study of chiral liquid crystal droplets, for which results are also reported here, is supported by the Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering Division, Grant No. DE SC0004025. 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 the Postdoctoral Fellowship under 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 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: FLUCTUATIONS IN A GINZBURG-LANDAU RELAXATION

During a GL relaxation, the symmetric and traceless projector operator is used to minimize the free energy following an Euler-Lagrange scheme, i.e.,

for the bulk and surfaces (boundary conditions), respectively, and where ** ν** is the unit outward normal to the surface. The Volterra derivatives are defined by

The solution to these equations is found by allowing the tensor order parameter **Q** to evolve towards equilibrium, according to an isothermal relaxation dynamics as follows:

with boundary conditions given in Eq. (A1), and where *γ* is a rotational viscosity (or diffusion) coefficient.^{22,23,37} The GL evolution equation may be written as follows:

with boundary conditions

where subscripts V and VS represent the Volterra derivatives and the corresponding boundary operator, i.e.,

and

where Einstein index notation is used.

For the GL relaxation, the time is scaled by the nematic relaxation time $\tau =\gamma \xi N2/L1$, while the elastic constants are scaled by *L*_{1}, i.e., $L\u02c65=L5/L1$. For the boundary conditions, *θ* = 1 imposes homeotropic anchoring, while *θ* = 0 is for degenerate planar anchoring. The anchoring strength is scaled as $W\u02c6=W\xi N/L1$. Similar to the MC, the GL relaxation is carried out using the orthonormal basis tensor set, where all the expressions containing *Q _{ij}*, ∇

*Q*, and ∇

_{ij}^{2}

*Q*are transformed into

_{ij}**a**, ∇

**a**, and ∇

^{2}

**a**.

The GL evolution is done following a Galerkin FE method in the weighted residual approach. When applied to Eq. (A4) results in

where **M** = *M*_{α} are the FE shape functions. After numerical integration and using a semi-implicit Euler scheme for the time derivative,^{53} a linear system of equations are obtained for the time evolution as follows:

where

The resulting coefficient matrix **A**, similar to any FE method, is sparse. Therefore, to gain computational efficiency, SuperLU libraries are used to decompose and back-substitute the linear system in time.^{80–82}

Fluctuations are included defining a rotational diffusivity *D _{R}* and by recognizing that thermal noise can be included directly into the GL relaxation following a Langevin-like scheme,

^{31,58}i.e.,

and

where $T\u02c6$ is the fluctuating temperature, *ζ _{R}* is the LC rotational drag coefficient, and

**𝒬**is a symmetric and traceless Gaussian white noise, which satisfies the fluctuation-dissipation theorem by enforcing

where 〈〉 denotes the average over the probability distribution of the noise. The Langevin-like noise **𝒬** may also be expressed in the orthonormal basis set as **𝒬** = ∑𝔮_{j}**T**^{j}, which enforces the symmetric and traceless character of the tensor and reducing the fluctuation-dissipation theorem to

An Îto integration^{83,84} is performed for **𝒬** into the GL numerical scheme, resulting in

where *d*𝔮_{j} have zero mean and variance $dt$. Similar to the MC annealing, *D _{R}* is decreased during the GL relaxation until no fluctuations are included and the global minima are reached.

Figure 9 shows the evolution of the free energy during minimization of a droplet with homeotropic anchoring using Ginzburg-Landau and theoretically informed Monte Carlo relaxations. Two different approximation methods are considered: finite differences and finite elements. The droplet with *N* = 0, *R* = 250 nm, *W* = 1 × 10^{−3} J/m^{2}, and *U* = 5 was selected because, for this simple case, all methods should arrive at the same energy minimum. Recall that high chiral systems require MC, while complex geometries require FE. The figure includes iso-surfaces representing low values of the scalar order parameter *S*. As shown in Fig. 9, both relaxations, GL and MC, and both approximation methods, FD and FE, reach the same energy minima with the corresponding LC phase within the droplet. Each method follows a different minimization path, but they are able to reach the minima at similar times, during a GL, and similar MC steps, during MC. For reasons discussed in the main manuscript, FE methods require fewer nodes than FD. In this particular case, less than half-FE nodes were enough to reach the same energy minima obtained by FD. An important difference between FE-GL and FD-GL relaxations is the fact that a fully explicit time integration (Euler) will result in a “matrix-free” FD-GL, whereas the FE-GL must resolve an explicit matrix, which needs a direct or iterative solver. For the results in Fig. 9, a sparse LU decomposition solver was used.^{80–82} We found that this minimization can be accelerated by the addition of the fluctuating terms, i.e., the rotational diffusivity. On average, the computational time of a non-fluctuating FE-GL was ten times longer than that of a fluctuating relaxation. Therefore, during the FE-GL, an initial rotational diffusivity of *D _{R}* = 0.1 was used and decreased, at intervals of 0.01 every 5000 time steps, until zero was reached (no fluctuations). This is why the FD-GL and FE-GL relaxations differ in shape in Fig. 9; the FD-GL evolves smoothly to the minima, whereas the FE-GL decreases the energy by a stepping process that corresponds to the decrease in the rotational diffusivity. As discussed in our previous work,

^{47}an important difference between GL and MC relaxations is the evolution, towards the energy minima, of the scalar,

*S*and

*η*, and vectorial,

**n**and

**n**′, order parameters. A typical GL minimization occurs in a

*S*→

**n**process, where the scalar order parameter is relaxed first followed by a slow director vector relaxation. Conversely, a MC minimization follows a

**n**→

*S*process, which results in a faster and un-trapped free energy relaxation. Including fluctuations in the FE-GL forces the method to follow an

**n**→

*S*process as well and helps to explain the source of the accelerated minimization.