Computing equilibrium states of cholesteric liquid crystals in elliptical channels with deflation algorithms

TL;DR

This paper investigates equilibrium configurations of cholesteric liquid crystals in elliptical channels, employing deflation algorithms to identify multiple solutions and analyze how these configurations change with geometric and material parameters.

Contribution

It introduces the application of deflation and deflated continuation algorithms to study complex equilibrium states in cholesteric liquid crystals within elliptical geometries.

Findings

Multiple equilibrium configurations are identified for different aspect ratios and pitches.

Deflation algorithms effectively find diverse solutions in frustrated liquid crystal systems.

Solution branches vary significantly with geometric and material parameters.

Abstract

We study the problem of a cholesteric liquid crystal confined to an elliptical channel. The system is geometrically frustrated because the cholesteric prefers to adopt a uniform rate of twist deformation, but the elliptical domain precludes this. The frustration is resolved by deformation of the layers or introduction of defects, leading to a particularly rich family of equilibrium configurations. To identify the solution set, we adapt and apply a new family of algorithms, known as deflation methods, that iteratively modify the free energy extremisation problem by removing previously known solutions. A second algorithm, deflated continuation, is used to track solution branches as a function of the aspect ratio of the ellipse and preferred pitch of the cholesteric.

Equations

$$F$$

$$+K_3{\left|\mathbf{n}\times \nabla \times \mathbf{n}\right|}^2$$

$$M_{p,\alpha }(\mathbf{u};\mathbf{v})=\left(\frac{1}{{\Vert \mathbf{u}-\mathbf{v}\Vert }_U^p}+\alpha \right)\mathbf{I},$$

$$\mathcal{B}(\mathbf{u})=M_{p,\alpha }(\mathbf{u};\mathbf{v})\mathcal{A}(\mathbf{u})=0.$$

$$M_{p,\alpha }(\mathbf{u};{\mathbf{v}}_1,{\mathbf{v}}_2)=M_{p,\alpha }(\mathbf{u};{\mathbf{v}}_1)M_{p,\alpha }(\mathbf{u};{\mathbf{v}}_2),$$

$$M_{p,\alpha }(\mathbf{u};{\mathbf{v}}_1,{\mathbf{v}}_2,...,{\mathbf{v}}_m)=\prod _{i=1}^m{M}_{p,\alpha }(\mathbf{u};{\mathbf{v}}_i),$$

$$f(x)=\frac{1}{5^4}{x}^4-\frac{1}{5^2}{x}^2+1,$$

$$f^{\prime }(x)=\frac{4}{5^4}{x}^3-\frac{2}{5^2}x=0.$$

$$M_{p,\alpha }(x;x_0)=\frac{1}{{\left|x-x_0\right|}^p}+\alpha ,$$

$$0=f_d^{\prime }(x)$$

$$Q=\frac{1}{4\pi }\int \mathbf{n}\cdot \left(\frac{\partial \mathbf{n}}{\partial x}\times \frac{\partial \mathbf{n}}{\partial y}\right)dA,$$

Full text

I Introduction

Cholesteric liquid crystals are complex fluids that exhibit long-range orientational order, elasticity, local alignment at surfaces, optical activity and response to external stimuli Bahr and Kitzerow (2001). They are composed of chiral molecules that, in the absence of boundaries, adopt a helical structure with a preferred pitch, $q_{0}$, set by the molecular structure and the ambient temperature. There has recently been a great deal of interest in cholesterics in confined geometries because of parallels with other condensed matter systems such as chiral ferromagnets, Bose-Einstein condensates and the Quantum Hall effect. All of these systems exhibit topological defects under confinement, including skyrmions and torons, where the boundary conditions preclude adoption of the energetically preferred uniformly twisted state. Hence, they are geometrically frustrated. It was recognised some time ago that nematic liquid crystals also may potentially form skyrmions, but these are only metastable due to the lack of preferred twist Bogdanov et al. (2003).

Liquid crystals are particularly attractive to study these defect structures, because they can be conveniently produced and imaged in three dimensions with techniques such as confocal microscopy Smalyukh et al. (2010). Cholesterics may form skyrmion lattices in two dimensions Fukuda and Žumer (2011). In three dimensions, torons resemble particulate inclusions Ackerman et al. (2014a); Leonov et al. (2014) and form chains or lattices Ackerman et al. (2015). Other more complicated defect structures called “twistions” occur in films thinner than the cholesteric pitch Ackerman and Smalyukh (2016). They also provide an elegant experimental actualisation of the Hopf fibration, a map from the 3-sphere onto the 2-sphere such that each distinct point of the 2-sphere comes from a distinct circle of the 3-sphere Chen et al. (2013). Strongly confined geometries such as micropatterned surfaces Cattaneo et al. (2016), channels Kim et al. (2015a); Guo et al. (2016a) and droplets Orlova et al. (2015) can all be used to control and order the location of skyrmions.

A key challenge in simulating these systems is that, due to the geometric frustration, they possess a particularly rich family of stationary solutions of the free energy. The ground state strongly depends on the shape of the domain and material parameters, including the cholesteric pitch. Typically, extremisation is performed from an initial guess using a relaxation algorithm or by solving a set of nonlinear Euler-Lagrange equations. In either case, having found a solution, the question remains: are there others? It is also highly desirable to track the solution set as a function of geometric and material parameters to assemble a bifurcation diagram.

A common approach to identify distinct solutions, referred to in the mathematics literature under the umbrella term of multistart methods Martí (2003), is to begin from a significant number of initial guesses. This requires extensive knowledge of the problem to devise a suitable set of initial guesses and can be inefficient as multiple guesses often converge to the same configuration. Other well-established techniques include numerical continuation Chao et al. (1975); Allgower and Georg (1993), which is particularly effective in fully resolving connected bifurcation branches but can neglect solutions if they are not homotopic with respect to the continuation parameters Farrell et al. (2016), and approaches, such as Branin’s method, that rely on numerical integration of the Davidenko differential equation corresponding to the original nonlinear problem Davidenko (1953); Branin (1972). Branin’s method is capable of systematic computation of multiple solutions but requires determinant calculations that become prohibitively expensive for large-scale problems without special structure.

In this paper, we adapt a new technique known as the deflation method Farrell et al. (2015); Adler et al. (2017) to a model problem in this class, the configuration of a cholesteric in an elliptical channel. The method is generalisable, robust and computationally efficient for large-scale applications. It has been successfully applied to a diverse set of nonlinear problems including nonlinear partial differential equations (PDEs), singularly perturbed problems, the analysis of Bose–Einstein condensates, and the computation of disconnected bifurcation diagrams Farrell et al. (2015); Chapman and Farrell (2017); Charalampidis et al. (2017); Farrell et al. (2016). This paper is structured as follows: in Section II, we briefly describe the problem; in Section III, we introduce the deflation technique with an illustrative toy example. In Section IV, we present results for the cholesteric problem. Conclusions are drawn, with prospects for further applications of the algorithm, in Section V.

II Model

We consider a cholesteric liquid crystal in a channel with elliptical cross section. Equilibrium structures are found by identifying critical points of the Frank energy,

[TABLE]

where $K_{1}$, $K_{2}$ and $K_{3}$ are the splay, twist and bend elastic constants; $\mathbf{n}$ is a headless unit vector, the director, that corresponds to the local symmetry axis of the molecular orientational distribution and $q_{0}$ is the preferred pitch for the cholesteric. Rigid anchoring (Dirichlet) boundary conditions are imposed on the boundary $\partial\Omega$, where the director is required to point perpendicular to the plane of the cross section. The energy is readily non-dimensionalised by introducing a typical length $\lambda$ and dividing through by a characteristic magnitude of the elastic constants $\tilde{K}$; henceforth, we use dimensionless parameters.

As discussed in the introduction, this problem promotes the existence of multiple local equilibria by construction. To see why, first consider the cholesteric in the absence of boundaries. As is well-known, the minimiser of (1) is a unique uniformly twisted state. Level sets of $\mathbf{n}$ form families of equally spaced planes often referred to as cholesteric “layers.” We note a valuable discussion of the limitations of this view is found in Beller et al. (2014). Variation away from this preferred structure, which is equivalent to compressing or bending the layers, implies an elastic cost. Considering a cholesteric in a disk, the minimisers of (1) are solutions where $\mathbf{n}$ rotates about the radial axis and lies everywhere perpendicular to it. The number of rotations is determined by the cholesteric pitch, $q_{0}$, which promotes a constant rotation rate, and level sets of constant orientation therefore form equally spaced concentric circles.

For an elliptical domain, however, it is not possible to fill the ellipse with equally spaced layers, and so defects or a variable layer spacing must be introduced. The cholesteric order, which prefers a uniformly twisted state, and the shape of the domain are in competition, so the system is said to be geometrically frustrated. The frustration is resolved by adopting a compromise state, incorporating some combination of layer deformation or defects; in common with other frustrated systems there is typically more than one way to do this, leading to the possibility of more than one minimiser.

Further, we explore solutions where $K_{2}>K_{1},K_{3}$, which might occur in exotic liquid crystal systems Le et al. (2015). This choice of parameters leads to a material that is doubly frustrated because it is required to twist by the cholesteric term but the twist is relatively expensive compared to other deformation modes. As a result, the cost of modulating the cholesteric layers is reduced. The interaction of geometric and internal frustration is expected to lead to a particularly rich solution set, because they permit multiple ways of relieving the frustration: one solution might accommodate an incommensurate number of cholesteric periods by folding the layers; another might introduce a defect. These parameters therefore yield an extremisation problem that we anticipate a priori to be very challenging to explore by naive multistart methods.

III Deflation

In solving problems that possess multiple equilibrium solutions, such as that posed in Section II, a key challenge is to ensure that the true ground state has been found from the set of energetically low-lying solutions. The idea of the deflation algorithm is to successively modify the nonlinear problem under consideration to eliminate known solutions, enabling the discovery of additional distinct solutions. Consider $\mathcal{A}(\mathbf{u})=\mathbf{0}$, a set of nonlinear equations, that may admit multiple solutions. This system, for instance, could represent a set of continuous or discretised nonlinear PDEs; here we minimise the Frank energy in (1), subject to the unit-length constraint on the director, and consider the resulting system of nonlinear first-order optimality conditions. Using a known solution, $\mathbf{v}$, to the system $\mathcal{A}(\mathbf{v})=\mathbf{0}$, we define the deflation operator,

[TABLE]

where $\alpha\geq 0$ is a shift scalar, $p\in[1,\infty)$ is the deflation exponent, and $\mathbf{I}$ is the identity operator. An appropriate norm $\left\|\cdot\right\|_{U}$ must be chosen for the vector space to which the solutions belong. The deflated system is then formed by applying the deflation operator to the original nonlinear system as,

[TABLE]

Iterative techniques, such as Newton’s method, may then be applied to solve the deflated system. While these iterations are guaranteed to not converge to the known solution $\mathbf{v}$ under mild regularity conditions, the remainder of the solution space for the original system is preserved by the deflation operator. Numerical experiments have found the effectiveness of the deflation operator to be relatively insensitive to the choice of deflation parameters. However, for certain problems, performance improvements may be obtained by varying $p$ and $\alpha$ Farrell et al. (2015); Adler et al. (2017). Typical values, used everywhere in this paper, are $\alpha=1$ and $p=2$.

Having found two solutions $\mathbf{v}_{1},\mathbf{v}_{2}$, an expanded deflation operator can be constructed by composition of single-solution deflation operators,

[TABLE]

and applied to the original nonlinear system. As the set of known solutions $\{\mathbf{v}_{1},\mathbf{v}_{2},...,\mathbf{v}_{m}\}$ is expanded, the deflation operator is grown as the product of the single deflation operators for each distinct solution in the set,

[TABLE]

and its action remains multiplicative on the original system.

To provide a simple and tractable illustration of the deflation process, we apply it to the problem of locating critical values of a one dimensional objective function,

[TABLE]

displayed in Fig. 1 by solving the equation,

[TABLE]

Starting from the initial guess $x=2.2$, Newton’s method locates the first solution $x_{0}=\sqrt{5}/2$. The deflation operator is constructed following the definition in (2), as

[TABLE]

where $\left|\cdot\right|$ denotes the standard absolute value. Applying this to (7) yields the deflated optimality condition,

[TABLE]

to be solved for $x$. The function $f^{\prime}_{d}(x)$ is also plotted in Fig. 1 for deflation parameters $\alpha=1$ and $p=2$. Notice that $x_{0}$ is not a solution to $f^{\prime}_{d}(x)=0$, while the remaining solutions to $f^{\prime}(x)=0$ persist as solutions to $f^{\prime}_{d}(x)=0$. Use of the deflation operator precludes convergence of certain iterative techniques, such as Newton’s method, to $x_{0}$ while facilitating convergence to additional distinct solutions from the same initial guess. With the deflation parameters chosen previously and the same initial point, $x=2.2$, Newton’s method converges to the critical point $x_{1}=0.0$. Thus, two solutions are obtained from the same initial guess. The process may then be repeated by constructing a multi-deflation operator, incorporating both known roots, to enable the discovery of the third distinct solution at $x_{2}=-\sqrt{5}/2$ to (7) and hence identify all extremal values of $f(x)$.

While deflation is a useful device for finding distinct solutions, the number of solutions discovered may still depend on the analyst supplying a suitable set of initial guesses. A systematic way to generate the set of initial guesses to use is provided by continuation. Suppose that the problem incorporates some set of parameters $k$, which for the liquid crystal problem includes the elastic constants and preferred pitch $q_{0}$. Given a set of solutions for some initial value of these parameters $k_{0}=\bar{k}_{0}$, we use each of these solutions as an initial guess for Newton’s method at a nearby parameter $k_{0}=\bar{k}_{0}+\delta k$, deflating each solution as we find it. Subsequently, we use the full power of the deflation approach to search for new solution branches that, if discovered, can be extended to other values of $k_{0}$ using standard continuation techniques. The solutions at $k_{0}=\bar{k}_{0}+\delta k$ can in turn be used as initial guesses to find the solutions at $k_{0}=\bar{k}_{0}+2\delta k$, and so on. This combination of deflation and continuation is referred to as deflated continuation and is an even more powerful algorithm than standard deflation applied to a single nonlinear problem Farrell et al. (2016). It can be interpreted physically as computing a bifurcation diagram, a portrait of the solutions of an equation as a parameter varies. We shall use both the deflation and the deflated continuation approach to resolve ground state solutions of the cholesteric problem in the next section.

IV Results

We apply the deflation algorithm described above to the cholesteric problem in an ellipse, varying the aspect ratio of the domain, $\mu$, and preferred cholesteric pitch, $q_{0}$. In each simulation, the director is held fixed on the boundary such that the director points out of the plane, i.e. $\mathbf{n}=(0,0,1)$. Elastic constants are chosen to be $K_{1}=1$, $K_{2}=3.2$ and $K_{3}=1.1$, corresponding to the exotic splay-bend cholesteric described above. The computational domain is centered on the origin with major axis parallel to the $x$-axis. The area of each bounding ellipse is held fixed at $\frac{3\pi}{2}$.

Our multilevel finite-element code used to compute stationary points of the free energy (1) is thoroughly described elsewhere Adler et al. (2015a, b, 2016, 2017). Briefly, the code uses the Cartesian representation of the director $\mathbf{n}=(n_{x},n_{y},n_{z})$ and directly finds equilibrium points of the Frank energy (1) by applying Newton linearisation to the first-order optimality conditions in variational form, resulting from the constrained minimisation. The code is based on deal.II Bangerth et al. (2007) and features mesh refinement and nested iteration Starke (2000), such that the problem is discretised and solved first on coarse meshes where computation is cheap, resolving major solution features, and then interpolated to more refined meshes to determine finer attributes of the computed approximation. Nested iteration has been shown to significantly improve computational efficiency in a wide variety of problems including liquid crystal simulations Adler et al. (2016, 2010); Manteuffel et al. (2006). Here, we use a nested iteration mesh hierarchy of four mesh levels of refinement with $4,884$ degrees of freedom on the coarsest grid and ending with 297,988 degrees of freedom at the finest level. Finally, a damping factor, $\omega$, is applied to each Newton step for both the undeflated and deflated systems. This damped Newton stepping is combined with an increased step size, $\bar{\omega}$, when the nonlinear residual drops below 0.1. The accelerated Newton stepping is applied to increase the rate of convergence when a candidate solution begins to closely satisfy the optimality conditions.

As a first example, in Figure 2 we display the results of a typical run for aspect ratio $\mu=1.5$ and $q_{0}=8$. The algorithm is initialised from three initial guesses (shown in the left column of Fig. 2). As anticipated, deflation enables the discovery of several solutions for each value by successively removing them with the deflation operator. Several of these solutions possess energetically degenerate partners that are obtained by simple reflection about an axis of the ellipse. These are also found by deflation, even though knowledge of the symmetry of the problem is not explicitly built into the algorithm.

It is important to note that the solutions to which the code converges are stationary points of the Lagrangian (Frank energy plus unit-length constraint), not necessarily minimisers of the energy. It is therefore highly desirable to characterise the nature of each solution as it is found, i.e. determine if it is a local minimum, a local maximum, or a saddle point. To do this, we must verify the second-order sufficiency conditions: a stationary point is a local minimum if the reduced Hessian of the energy (the Hessian projected onto the nullspace of the linearised constraints, i.e. restricted to feasible perturbations) is positive-definite. One approach would be to assemble the linearised constraint Jacobian, compute a (dense) basis for its nullspace using the singular value decomposition, construct the (dense) reduced Hessian, and compute its eigenvalues; however, this would be unaffordably expensive. A better way is to exploit the fact that the eigenvalues of the reduced Hessian of the energy are related to the eigenvalues of the (sparse) Hessian of the Lagrangian: by counting the number of negative eigenvalues of the Hessian of the Lagrangian, and comparing it to the dimension of the constraint space, we can determine the inertia (the number of positive, zero, and negative eigenvalues) of the associated reduced Hessian of the energy (Nocedal and Wright, 2006, Thm. 16.3). This allows for the characterisation of the nature of each solution found using a single sparse $LDL^{T}$ decomposition, computed using the FEniCS, PETSc and MUMPS software packages Logg et al. (2011); Balay et al. (2015); Amestoy et al. (2001).

In Figure 3, we show the computed ground state (lowest-energy solution) as a function of $\mu$ and $q_{0}$. For each solution, we display the value of the energy functional and also compute the skyrmion number Bogdanov et al. (2003),

[TABLE]

a topological index that represents the number of times $\mathbf{n}$ covers the unit sphere. Such indices help identify topologically distinct solutions: as the parameters $\mu$ and $q_{0}$ are slowly varied in Fig. 3, the ground state mostly changes smoothly. However, between certain values, e.g. $q_{0}=4$ and $5$ with $\mu=1$, a transition to a new solution as the ground state occurs; this is accompanied by a change in the skyrmion number. Some solutions that are quite different in appearance e.g. $\mu=1.85$ between $q_{0}=4$ and $5$ or $q_{0}=8$ and $9$ have identical $Q$ because they are linked by a continuous deformation of the director field. While deflated continuation enables us to find intermediate solutions between chosen values, and resolve transitions that take place, the skyrmion number provides a partial classification of the distinct branches that arise.

For certain values of $\mu$ and $q_{0}$, the solution set discovered by deflation alone failed to include a minimal energy solution that was stable. These values are indicated in Figure 3 with an asterisk. For these values, we applied the deflated continuation technique described above to identify the stable ground state shown in Figure 3. For instance, consider the case $\mu=1.15$ and $q_{0}=7$. The lowest energy solution found using deflation is displayed as an inset indicated by an asterisk in Fig. 4, but possesses unstable directions according to the Hessian analysis described above. We therefore use the $\mu=1.15$ and $q_{0}=6$ solution set and continue in $q_{0}$ from $6$ to $7$. The energies of intermediate solutions obtained in the process are plotted in Fig. 4 as a bifurcation diagram; the initial and final solution set recovered in this process are displayed in Fig. 4 as insets. A new pair of solutions, not in the $q_{0}=6$ solution set, emerges through a fold bifurcation at approximately $q_{0}\approx 6.67$ and represents the stable ground state at higher values of $q_{0}$; the prior ground state becomes higher in energy and is now unstable.

The same procedure was applied to all the problematic cases in Fig. 3, where the solutions shown are the lowest energy solutions found and are all verified as stable. This example illustrates the power of deflated continuation to track distinct branches in the solution set and identify solutions very different from the initial guesses provided. While it remains possible that the true ground state remains elusive for some values in $(\mu,q_{0})$ space, it is clear that deflation and deflated continuation are powerful tools to assist in the assembly of phase diagrams.

The solutions found in Fig. 3 catalogue the interesting interplay of the elastic constants, cholesteric parameter, and confining geometry. For $\mu=1$, a circular domain, increasing $q_{0}$ initially leads only to the incorporation of additional rotations as expected. Around a critical value of $q_{0}=7$, the contours of constant orientation are greatly deformed as the number of radial rotations in the channel increases from $\pi$ at $q_{0}=6$ to $3\pi/2$ for $q_{0}=8$. A similarly deformed structure is visible at $q_{0}=10$, which is apparently close to a jump from $3\pi/2$ rotations to $2\pi$.

For higher aspect ratios, the contours of constant orientation can be deformed by the geometry of the channel. For example, for aspect ratio $\mu=1.5$ and $q_{0}\leq 7$, the ground state consists of the director rotating by $2\pi$ about the radial direction from the center; above this value an extra $\pi$ rotation is incorporated. Comparing the shape of contours of constant tilt, notice that the interior “layer” is approximately circular for $q_{0}=5$, but becomes more elongated with increasing $q_{0}$. For $q_{0}=8$, the ground state is strikingly different: a highly deformed interior layer is accommodated within one contiguous outer layer. The ground state for $q_{0}=9$ reverts to the expected pattern, simply incorporating an additional twist. Furthermore, as $\mu$ increases, the transition points between states with different amounts of twist occur at higher values of $q_{0}$, and are typically proceeded by substantial layer deformation. Therefore, the confining geometry plays a role in deterring or encouraging the addition of layers.

We note that deflation uncovers a particularly large number of solutions for $q_{0}=8$ and that many of the solutions have relatively low energy compared to the ground state. For other values of $q_{0}$, only a few of the solutions are close to the ground state in energy. We speculate that this phenomenon is due to commensurability, with some values of $q_{0}$ being more compatible with the shape of the domain than others. For a circular domain, commensurate solutions exist where $q_{0}$ happens to allow an integer multiple of $\pi$ rotations from the center. Maximally strained solutions occur between these special values, potentially inducing deformation of the layers to relieve the frustration. This is clearly visible at $\mu=1.5$ and $q_{0}=8$, or at $\mu=1.85$ and $q_{0}\geq 9$, where the inner layer in both cases is highly tortuous to fill the interior of the domain.

To resolve the sequence of transitions that occurs around one of the maximally strained solutions, we visualise a bifurcation diagram in Figure 5 for an elliptical domain with aspect ratio $\mu=1.5$ and with the preferred pitch ranging from $q_{0}=5$ to $q_{0}=9$. We initialise the computation with the solutions found by our previous analysis. The diagram shown in 5A displays a relatively small solution set for $q_{0}<7$, but above $q_{0}\approx 7.4$, a dense thicket of additional solutions appears. As is visible in the expanded region depicted in 5B, this consists of a rapid succession of fold and pitchfork bifurcations. Moreover, the region between $q_{0}=7.4$ and $9.0$ contains a notable number of stable solutions with free energies in relatively close proximity to the ground state. The corresponding solutions at several values of $q_{0}$ are displayed in 5C, illustrating the striking complexity of solution sets that can be uncovered using deflated continuation.

Because deflation suppresses solutions close to the initial guess, even very distant solutions can be recovered. A particularly important result is that, while our initial guesses are smooth functions, the algorithm was able to spontaneously identify solutions with disclinations. Confined cholesterics in rectangular domains or channels have been experimentally and numerically shown to exhibit structures with disclinations Kim et al. (2015b); Ackerman et al. (2014b); Guo et al. (2016b), as discussed in the introduction. The existence, type, and number of the disclinations were shown to be modulated by changes to the depth and width of the channel, as well as the cholesteric pitch. Figure 6 displays three examples for $\mu=1.5$ and varying $q_{0}$ with computed director fields and associated elastic energy densities from which disclinations are readily identified. For $q_{0}=5$, the solution shown in Fig. 6A has four defect points, arranged in a diamond pattern near the center of the domain. The solutions in Fig. 6B and C were found for $q_{0}=7$ and possess two and one disclinations, respectively. Comparing the solutions’ free energies, it is clear that the single defect structure is energetically preferred. We note that the energy of the solution in 6C corresponds to the third lowest energy (and first unstable state) shown for $q_{0}=7$ in the bifurcation diagram in Figure 5A. Thus, our results suggest that the propensity of the cholesteric to forming metastable structures with disclinations in elliptical channels, perhaps upon quenching from the isotropic phase, strongly depends on the aspect ratio of the channel boundary. Moreover, our numerical experiments indicate that multiple equilibrium configurations with distinct disclination patterns may exist for the same geometry and material parameters.

The deflation technique plays a central role in the discovery of these disclination arrangements. For instance, numerical simulations in Guo et al. (2016b) relied on a priori knowledge, gained from experimental observations, that disclination structures should be present in order to initialise the Newton iterations within a basin of attraction. We emphasise that here the simulations are initialised with smooth director fields. Multiple solutions and the emergence of disclinations occurred as spontaneous discoveries enabled by the deflation computations. In situations where experimental and analytical information is limited, such numerical capabilities facilitate a more robust and thorough exploration of the admissible solution space of a given problem.

V Conclusion

In this paper, we present a new technique, deflation, for recovering equilibrium solutions of the free energy of a liquid crystal. The utility of the method is shown on a toy example and then used to determine the structure of a cholesteric in an elliptical domain. The ground state is identified for a range of aspect ratios $\mu$ and preferred pitches $q_{0}$, showing gradual deformation of the solutions as a function of these parameters and transitions to different solutions at critical values. For selected values of $\mu$ and $q_{0}$, we compute the bifurcation diagram, finding remarkably dense solution sets near the transition points. In future work, we will apply the method to characterise the solution set of more complex geometries involving cholesterics, such as the rich 3D structures observed in Smalyukh et al. (2010).

The deflation methodology significantly enhances the utility of Newton iterations applied to nonlinear systems by enabling Newton’s method to converge to multiple solutions from the same initial guess. Applying the deflation operator ensures that subsequent Newton iterations do not converge to previously discovered solutions and effectively modifies the basin of attraction to include unknown solutions. While convergence to multiple solutions is not guaranteed through use of the deflation operator, sufficient conditions for convergence of deflated iterations are constructed in Farrell et al. (2016) based on a generalisation of the Rall-Rheinboldt theorem. Our results illuminate the effect of deflation, but also highlight the case that the use of multiple initial guesses to discover all solutions is not entirely eliminated with deflation. However, deflated continuation systematically provides a sequence of good initial guesses and increases overall efficiency and reliability of multiple solution discovery through more systematic exploration of the solution space.

More generally, deflation and deflated continuation therefore allow theorists to recover energetically low-lying solutions in which a liquid crystal may become kinetically trapped. It can also be used to track different solution branches when a system exhibits a bifurcation with respect to some external parameter, e.g. the applied field in a Freedericksz transition. The method is very general and can be readily adapted to different representations, e.g. the Q-tensor, and parametrisations; it is likely to be most useful in systems where little analytical guidance or experimental imaging are available.

Acknowledgements.

The work of SPM was partially supported by an NSERC discovery grant. TJA is supported by a Cottrell Award from the Research Corporation for Science Advancement, and a CAREER award from the National Science Foundation (award number DMR-CMMT-1654283). PEF is supported by an EPSRC Early Career Research Fellowship (award number EP/K030930/1). The authors contributed to the work as follows: SPM devised the toy example in Section III; DBE, PEF, and JHA performed the simulations in Section IV and all authors contributed to their analysis; TJA wrote the paper, which was revised by all authors.