Nonlinear Methods for Shape Optimization Problems in Liquid Crystal Tactoids

Absent any external forces, we specifically consider the one constant Landau–de Gennes model, where the free energy sufficient to capture the physics of nematic tactoids, $\mathcal{F}$, depends on the state variables of the system: the shape $\mathcal{M}$ and the tensor $\mathcal{Q}({\bf{x}})$. The free energy is,

where the first three bulk terms represent a Landau expansion of the free energy of ${\mathcal{Q}}$ with Landau coefficients $a$, $b$, and $c$ with units $Nm^{-2}$. These parameters, in effect, select a particular uniform value of $S$ in the bulk; by convention, $a$ is chosen to be temperature dependent, $a=a_{0}(T-T_{0})$, where $T_{0}$ is the temperature below which the isotropic phase is no longer stable. The fourth term involving $\nabla{\mathcal{Q}}$ represents elasticity, and here, we adopt the commonly used one-constant approximation with a single elastic constant $L_{1}$ with units $N$. On the boundary, $\partial\mathcal{M}$, the constant term with prefactor $\sigma$ represents the isotropic surface tension. The second term is the anisotropic surface tension with anchoring coefficient $W$ and is constructed to favor the alignment of the LC in the tangent plane of the boundary surface. The surface parameters have units $Nm^{-1}$. Here, $\mathcal{Q}_{s}$ is a $Q-$tensor with alignment in the normal direction of the surface,

where $\bm{\nu}$ is the local normal vector at the surface, $S_{0}$ is the degree of order induced by the surface, which may differ from the value set by the Landau coefficients in the bulk depending on the chemistry of the liquid crystal-host interface. The prefactor $W>0$ controls the orientation such that, combined with the negative sign in front, higher values penalize the director’s orientation toward the tangent plane. As discussed above, we impose a volume constraint,

where, $V_{0}$ is the target volume of the tactoid.

In three dimensions, $\mathcal{Q}$ can be parameterized to be symmetric and traceless,

and hence includes five independent degrees of freedom $\{q_{xx},q_{xy},q_{xz},q_{yy},q_{yz}\}$. We nondimensionalize $\mathcal{F}$ and $\mathcal{C}$ by introducing a length scale $\xi$, so that ${\bf{x}}\to\xi\bar{{\bf{x}}}$ where $\bar{{\bf{x}}}$ is non-dimensional, and divide (2.1) by $L_{1}$ to get,

with dimensionless parameters,

The volume constraint, (6), is trivially non-dimensionalized, and hence we simply choose $V_{0}$ to be dimensionless. Note that $\mathcal{Q}$ and $S$ are also non-dimensional.

Hence, to find the equilibrium order and mesh configuration, we minimize (2.1) subject to the global nonlinear constraint, (6). For the rest of the paper, we assume all quantities are suitably nondimensionalized and drop the bar notation.

To set up the discrete optimization problem, we represent $\mathcal{M}\subset\mathbb{R}^{d}$ as a simplicial complex, $M$, consisting of $p$ spatial coordinate points, ${\bf{x}}_{i}\in M$, $i=1\ldots p$, and $N$ triangular ($d=2$) or tetrahedral ($d=3$) elements. We denote by ${\bf{X}}\in\mathbb{R}^{pd}$ the collection of all ${\bf{x}}_{i}\in M$. We then use finite elements to discretize the Q-tensor, (7), over $M$, representing each of the five components as a piecewise polynomial function over the elements of $M$. Note that this enforces the symmetry and tracelessness properties of the discrete Q-tensor, which we denote as $Q$.

Finite-element approaches for Q-tensor models have been considered in [39, 42, 43], as well as for other LC frameworks [44, 45, 46]. The algorithms we develop in this paper are independent of the specific finite-element spaces chosen as long as the systems remain well-posed. However, we choose a linear piecewise-polynomial representation for simplicity, noting that higher-order representations are possible when needed. Using such linear finite elements, the degrees of freedom for the Q-tensor are defined as $Q_{i}:=\mathcal{Q}({\bf{x}}_{i})$ at each vertex point ${\bf{x}}_{i}\in M$. Thus, in three dimensions, this discretization leads to $8p$ degrees of freedom: three for each coordinate of ${\bf{x}}_{i}$, and five for the $Q_{i}$, one for each $\{q_{xx},q_{xy},q_{xz},q_{yy},q_{yz}\}$ at ${\bf{x}}_{i}$.

Next, $\mathcal{F}({\bf{x}},\mathcal{Q})$ and $\mathcal{C}({\bf{x}})$ are defined over $M$ to obtain $F({\bf{X}},Q)$ and $C({\bf{X}})$,

where $\text{Vol}(T_{i})$ is the volume of the $i^{th}$ tetrahedron, $T$, in the simplicial complex. The corresponding discrete shape optimization problem is then defined as,

where again ${\bf{X}}\in M$ and $Q$ is a symmetric and traceless 2nd-order tensor whose components are piecewise linear functions defined on $M$. Here, ${\bf{X}}^{*}$ and $Q^{*}$ represent the equilibrium configurations of the mesh points in the manifold and the Q-tensor approximation on $M$ for the minimized configuration.

In [23], the authors note that the combined $({\bf{X}},{\bf{n}})$ optimization problem exhibited stiffness due to intrinsic differences in length scale between the vertices, ${\bf{X}}$, and directors, ${\bf{n}}$, degree of freedom. To alleviate this, they use an alternating gradient descent scheme, first taking descent steps in ${\bf{n}}$ followed by ${\bf{X}}$ and repeating until convergence. Given the connection between ${\bf{n}}$ and $Q$, we use a similar approach.

Consider a vector of the spatial vertices at the $k^{th}$ iteration, ${\bf{X}}^{k}\in M^{k}$, where $M^{k}$ is the current triangulation of the grid, $M$. We first compute the gradients of $F$ and $C$ with respect to ${\bf{X}}$ on the given grid $M$, denoted by $F^{k}_{M},$ and $C^{k}_{M},$ respectively, and evaluate each gradient at ${\bf{x}}_{i}\in M^{k}.$ This trio of data, $\{{\bf{X}}^{k},F^{k}_{M},C^{k}_{M}\}$ is then projected onto the constraint’s tangent space,

Here, $\alpha_{k}$ is the $k^{th}-$iteration step-size found from performing a one-dimensional line search. The update defined in (13) is an intermediate step towards finding ${\bf{X}}^{k+1}\in M^{k+1}$ since $\widetilde{{\bf{X}}}^{k}$ only satisfies the constraint to linear order. Following [23], we perform additional reprojection steps,

In each reprojection step, $C^{k}_{M}$ is evaluated at $\widetilde{{\bf{X}}}_{i}\in\widetilde{M}^{k}.$ The prefactor, $\frac{C(\widetilde{{\bf{X}}}^{k})}{C^{k}_{M}\cdot C^{k}_{M}},$ represents the difference between the constraint value and its target; reprojection steps are repeated until the constraint is satisfied to a given tolerance. This projected value is used as the next iterate ${\bf{X}}^{k+1}\leftarrow\widetilde{{\bf{X}}}^{k}$.

To find $Q$ we follow a similar procedure but with $\{Q^{k},F^{k}_{Q}\}$,

where $\beta_{k}$ is obtained with a separate one-dimensional line search. A second projection is unnecessary as the global constraint only depends on ${\bf{X}}\in M.$ In practice, [23] starts by optimizing for the field and then for the spatial coordinates. To control mesh quality, a procedure called equiangulation is used after a few iterations of each optimization routine. This ensures that no irregular triangles (i.e., long and skinny) appear in the discrete manifold, and is crucial when the amount of spatial or order deformation increases.

The GD algorithm introduced in the previous section converges quite slowly and can stagnate under certain conditions. It also incorporates a number of metaparameters, such as the mesh quality control and the number of alternating iterations to be taken for shape and field degrees of freedom that must be hand-tuned for each problem. To address these issues, we develop a quasi-Newton (QN) based method to solve the full $({\bf{X}},Q)$ minimization problem all at once. Given the presence of the nonlinear constraint, we introduce the Lagrangian $\mathcal{L}$ of the system based on (12),

where $\lambda\in\mathbb{R}$. The necessary first-order optimality conditions for minimizing the Lagrangian are

where $F_{M}$ and $F_{Q}$ are the first-order G$\hat{\text{a}}$teaux derivatives with respect to every ${\bf{x}}_{i}\in M$ and $Q_{i}$ on $M$, respectively (i.e., the gradients as computed for GD), and $C_{M}$ is the gradient of the constraint with respect to each ${\bf{x}}_{i}$.

These equations are nonlinear, so we linearize (17)-(19) and set up the corresponding iterative scheme. Let $({\bf{X}}^{k},Q^{k},\lambda^{k})$ be the current approximations for $({\bf{X}},Q,\lambda)$, respectively. The update, $({\bf{d}}_{M},{\bf{d}}_{Q},d_{\lambda})$, to the approximation is the solution to the QN iteration,

with

Here, ${\bf{B}}^{k}$ is an approximation to the Hessian of the $({\bf{X}},Q)$ portion of the Lagrangian, denoted by $\nabla^{2}\mathcal{L}^{k}:=\begin{bmatrix}\mathcal{L}_{MM}^{k}&\mathcal{L}_{MQ}^{k}\\ \mathcal{L}_{QM}^{k}&\mathcal{L}_{QQ}^{k}\end{bmatrix}$, where the entries are second-order G$\hat{\text{a}}$teaux derivatives with respect to the mesh and the $Q-$tensor at $({\bf{X}}^{k},Q^{k})$. Additionally, ${\bf{A}}^{k}$ contains the gradient of the constraint function at ${\bf{X}}^{k}$, and ${\bf{R}}^{k}$ is the nonlinear residual for the $({\bf{X}},Q)$ portion of the system. Since $\mathcal{A}$ is a saddle-point matrix, which poses challenges for building efficient solvers, we choose a BFGS approximation and use explicit formulae to compute ${\bf{H}}^{k}={({\bf{B}}^{k})}^{-1}$ and $\mathcal{A}^{-1}$ [47]. This allows us to find a closed form for computing the updates $({\bf{d}}_{M},{\bf{d}}_{Q},d_{\lambda})$ which we use to set up a recursive matrix-free approach of implicitly updating the approximation of ${\bf{H}}^{k}.$

Next, we find the step size for the field, $\alpha_{Q}$, and the step size for the spatial coordinates $\alpha_{M}.$ The step size $\alpha_{Q}$ is found with backtracking line search such that it satisfies,

where $\eta\in(0,1).$ Then, we update the field accordingly,

The step size $\alpha_{M}$ for the spatial coordinates is found using an $\ell^{1}$ monitor function [48] defined as

where $\mu\in\mathbb{R}$ is a penalty term for the constraint. This monitor function decides what factor of ${\bf{d}}_{M}$ and $d_{\lambda}$ should be accepted in order to decrease $F({\bf{X}}^{k+1},Q^{k})$. Similar to finding $\alpha_{Q}$, we use backtracking until the following inequality is satisfied,

where $\eta\in(0,1)$, $\mu_{k}>\|\lambda_{k}\|_{\infty}+\rho$ ($\rho$ is a positive constant), and

Finally, we update the spatial coordinates and Lagrange multiplier,

Given that the constraint is prescribed only to the mesh, we allow for its respective Lagrange multiplier, $\lambda^{k}$, to have the same step size $\alpha_{M}$ as the mesh update.

Finally, we note that quasi-Newton methods have been extensively studied, and several convergence results have been established. For the BFGS approach we consider, a well-known result [25, Thm. 18.6] guarantees superlinear convergence of the method under certain assumptions. Though this result does not consider line search, for all of the tests presented here, the step sizes approach one as the method converges. The assumptions of the theorem include conditions on the continuity and differentiability of the objective and constraint function that can be shown by direct calculation. Moreover, there are conditions on the linear independence of the gradients of the constraints in order to guarantee that the first-order optimality conditions, (17)-(19), are satisfied. Since we only have one constraint here, this is trivially true. We also need to assume that the Hessian, $\nabla^{2}\mathcal{L}$, at the local solution, $({\bf{X}}^{*},Q^{*})$ and its initial approximation, ${\bf{B}}^{0}$, are symmetric and positive definite. The latter is straightforward since ${\bf{B}}^{0}=\delta\mathcal{I}$ with $\delta>0$, while the former requires several assumptions on the physical parameters. The numerical experiments we present indicate that we are within the regime where this is satisfied. Finally, as with all Newton-based methods, we need the initial guess to be “good enough” in order for the iterates to remain in the basin of attraction. To guarantee this, we introduce the notion of nested iteration.

Both QN and GD, as described above, are sensitive to initial guesses. In particular, as the anisotropic surface tension parameters increase, we expect to see the nematic LC tactoids elongating towards the characteristic eye shape [14, 2]. Numerically, this means that by starting from the same initial guess for every parameter value, we may not be close to the basin of attraction for some regions of parameter space. This can lead the method to stagnate at locally optimal solutions, not the expected global minima. To remedy this issue, we wrap each method with nested iteration.

Nested iteration [27, 28] begins with an initial coarse grid, denoted by $M_{i}$ with only a few vertex points, $p$, that represents the problem domain. This coarse grid may not capture all the details of the solution, but solving the nonlinear system on this grid is computationally inexpensive. Thus, (12) is solved on the coarse grid with either QN or GD until a preferred convergence criterion is reached. The coarse grid, $M_{i}$, is then subdivided into smaller elements. In this work, we consider uniform refinement such that each element is divided into four (for triangular elements) or eight (for tetrahedral elements) smaller elements. The coarse solution from $M_{i}$ is then linearly interpolated onto the finer grid $M_{i+1}$ and used as an initial guess for solving the problem now represented on $M_{i+1}$ of size $8p$ (in three-dimensions) and $4p$ (in two-dimensions). This initial solution on $M_{i+1}$ should then be a more accurate guess as it came from solving a similar problem on $M_{i}$, and thus fewer iterations to converge are expected.

We begin by modeling two-dimensional nematic tactoids, where the computational domain represents a slice of an object that is infinitely extended in the perpendicular direction. We investigate solutions where the nematic director lies in-plane on the computational domain, and hence $\mathcal{Q}$ can be described by a reduced parameterization,

that includes only two independent degrees of freedom $q_{xx}$ and $q_{xy}$. We construct the preferred surface value of the Q-tensor ${\mathcal{Q}}_{s}$ in (2.1) from the local tangent vector ${\bf{t}}$ on the boundary,

To favor alignment with the local tangent vector, we must prefactor the anchoring term in (2.2) with a positive sign,

The volume constraint becomes a cross-sectional area constraint for a target area $A_{0}$,

Consequently, the constraint has the explicit discrete form,

We use the same parameters as those listed in (9) for the two-dimensional geometry. In this two-dimensional setting, we envision our domain as a $2D$ slab embedded in a three-dimensional space. Therefore, the energy functional (2.1) now has dimensions of energy per unit length or $Jm^{-1}=N$. Setting $\xi=1\times 10^{-7}~{}m$, then, the constants corresponding to those in the non-dimensionalized free energy, (2.2), are

For the numerical experiments below, we let the surface tension, $\overline{\tau}$, vary from $1$ to $100$ and the surface anchoring, $\overline{\omega}$ vary from $0.01$ to $1.$ This corresponds to a dimensionalized surface tension, $\sigma$, ranging from $10^{-4}-10^{-2}\,N$ and anchoring values, $W$, ranging from $10^{-6}-10^{-2}\,N$. Moreover, we artificially scale the Landau coefficients so that the defect size is $\sim 10$ times its true value. The scalar order parameter $S_{0}$ for the tangential anchoring is initialized as,

For the rest of the section, we refer to nondimensionalized quantities only and drop the bar notation.

For most test problems, the initial guess on the coarsest grid is defined on a circle of area $1$ with equally distributed vertices representing the degrees of freedom, ${\bf{X}}\in M$. The initial director field is aligned parallel to the $x-$axis. We start by solving the discrete problem on the coarsest grid of size $|M_{i}|=16$, where the notation $|M_{i}|$ denotes the number of vertices on grid $M_{i}$, and propagate the solution across nine levels of refinement until the finest grid of size $|M_{i}|=591,361$. For both methods described in this work, after each level of refinement, we perform the equiangulation procedure mentioned above, that improves the quality of the mesh, ensuring there are no irregular elements. The preferred convergence criterion is the relative change in the energy $F^{k}$, with the tolerance set at $10^{-6}$. We compare the results from the NI algorithm to those obtained running solely on the finest grid, $M_{9}$. We compare the final energy, $F^{k}$, the runtime in seconds, and the number of iterations on each grid.

We conclude the numerical results section by simulating the spatial and orientational configuration of a challenging three-dimensional problem involving the formation of a nematic tactoid. We use the same material constants as above [41], with $\xi=1\times 10^{-7}~{}m$ and, as before, we report the distribution of the scalar order parameter, $S$, the corresponding director field, the converged energy, $F^{k}$, the number of iterations for each grid level of the QN with NI scheme, and the runtime in seconds. For the numerical experiments below, we fix the surface tension, $\tau=10$, and vary the surface anchoring, $\omega$ from $0.01$ to $0.18$. The preferred convergence criterion is similarly the relative change in the energy $F^{k}$, this time with the tolerance set at $10^{-4}.$

Recall that $|M_{i}|$ denotes the number of vertices on grid $M_{i}$. Here in three dimensions, we consider 6 different grids of increasing size, ranging from $|M_{1}|=27$ to $|M_{6}|=275,793$. Note that there are eight degrees of freedom attached to each grid point, and the actual number of vertices on each level might vary depending on the aspect ratio of the current configuration.

For the first value of $\omega$ ($\omega=0.01$) we consider, the initial guess is defined on a sphere with volume $1$ with equally distributed vertices representing the degrees of freedom, ${\bf{X}}\in M_{1}$. The initial director field, ${\bf{n}}=(1,0,0)$, is aligned parallel to the $x-$axis. However, as we see numerically, the shape’s aspect ratio is linearly $\omega-$dependent, much stronger than the logarithmic $\omega-$dependence from the two-dimensional tactoid shapes seen earlier. Thus, the initial guess of a sphere is not adequate for convergence as $\omega$ is increased. To mitigate this, we incorporate continuation with respect to $\omega$ on the coarsest grid, $M_{1}$. From there, nested iteration is implemented as before. For example, the final equilibrium state on grid $M_{1}$ using $\omega=0.01$ is used as the initial guess for the simulation of $\omega=0.02$ on grid $M_{1}$. The final state for $\omega=0.02$ on grid $M_{1}$ is then used for the $\omega=0.04$ simulation on the coarsest grid and so on.