# Free energy pathways of a Multistable Liquid Crystal Device

###### Abstract

The planar bistable device [Tsakonas et al., Appl. Phys. Lett., 2007, 90, 111913] is known to have two distinct classes of stable equilibria: the diagonal and rotated solutions. We model this device within the two-dimensional Landau-de Gennes theory, with a surface potential and without any external fields. We systematically compute a special class of transition pathways, referred to as minimum energy pathways, between the stable equilibria that provide new information about how the equilibria are connected in the Landau-de Gennes free energy landscape. These transition pathways exhibit an intermediate transition state, which is a saddle point of the Landau-de Gennes free energy. We numerically compute the structural details of the transition states, the optimal transition pathways and the free energy barriers between the equilibria, as a function of the surface anchoring strength. For strong anchoring, the transition pathways are mediated by defects whereas we get defect-free transition pathways for moderate and weak anchoring. In the weak anchoring limit, we recover a cusp catastrophe situation for which the rotated state acts as a transition state connecting two different diagonal states.

## I Introduction

Nematic liquid crystals are complex anisotropic liquids with long-range orientational ordering de Gennes and Prost (1995). Nematics in confinement present a whole host of new theoretical and applications-oriented questions focussed on the complex inter-relationship between material properties, geometry, temperature, boundary effects and external fields in equilibrium and non-equilibrium phenomena. From a modelling point of view, we have new and challenging questions on pattern formation, defects, interfacial phenomena and rheology. From a technological point of view, nematics in micron-scale or nano-scale geometries open new doors for optical and display applications, and more recently novel biological insight Lagerwall and Scalia (2012); Blancab et al. (2013).

We re-visit the planar bistable device first reported by Tsakonas and his co-workers Tsakonas et al. (2007). The planar bistable device has been well-studied in the liquid crystal community in recent years Tsakonas et al. (2007); Luo et al. (2012); Kralj and Majumdar (2014); Anquetil-Deck and D.Cleaver (2010); Anquetil-Deck et al. (2012), partly because of its geometrical simplicity and partly because of its rich modelling landscape. It typically comprises a periodic array of shallow square or rectangular wells filled with nematic liquid crystalline material. The well surfaces are treated to induce tangent or planar boundary conditions so that the nematic molecules, in contact with the well surfaces, are constrained to be tangent to the well surfaces. Tsakonas et al. Tsakonas et al. (2007) adopt a two-dimensional Landau-de Gennes modelling approach and focus on planar nematic equilibria on the bottom square or rectangular well cross-section, on the grounds that non-planar three-dimensional equilibria are energetically expensive for shallow three-dimensional wells and hence, not physically relevant. They compute local minimizers of a two-dimensional Landau-de Gennes energy on a square or a rectangle with tangent boundary conditions and obtain two distinct classes of optically contrasting nematic equilibria: the diagonal solutions and the rotated solutions and their modelling results are supported by parallel experimental work. The tangent boundary conditions naturally create a mismatch in molecular alignments or defects at the square vertices. As the name suggests, the nematic molecules roughly align along the square diagonal for the diagonal solution. The rotated solutions have a more distorted profile as simulations and optical data suggest that the molecules rotate by about 180 degrees between a pair of parallel square edges. There are two diagonal solutions, one for each square diagonal, and four rotated solutions, related to each other by a 90 degree rotation. Luo et al. Luo et al. (2012) build on the work in Tsakonas et al. (2007) and carefully study the dependence of the diagonal and rotated solutions on the surface anchoring strength, denoted by a surface anchoring coefficient . In particular, they find that the rotated solutions only exist above a certain critical anchoring strength and the critical anchoring strongly depends on the material parameters and temperature.

Here we take the modelling work further by a systematic numerical investigation of the free energy pathways of the planar bistable device within the framework of the Landau-deGennes theory. This is a significant forward step since we not only recover the free energy minimizing states but also compute structural and energetic information about the transient states connecting pairs of distinct free energy minima. As we vary the surface anchoring strength, we find three distinct regimes: the strong anchoring, moderate anchoring and weak anchoring regimes. In the strong anchoring limit, we recover the familiar diagonal and rotated solutions as local free energy minimizers. As noted in Tsakonas et al. (2007); Luo et al. (2012), there are two different diagonal solutions and four distinct rotated solutions. We compute the transition pathways between the different stable solutions. In particular, we focus on the so-called minimum energy pathways, where every point on such a pathway is an energy minimum in all but a single distinguished direction in the phase space. There are a number of interesting findings to highlight here. Firstly, each minimum energy pathway is featured by a transition state, which is a saddle-point, an unstable critical point of the Landau-de Gennes energy. We identify the transition states as being local free energy maxima along a minimum energy pathway connecting free energy minima. In the strong anchoring limit, the transition states exhibit defects along the square edges. There may be multiple minimum energy pathways between a pair of distinct minima and these transition pathways can contain different number of defects. We interpret the optimal transition pathway as being the minimum energy pathway with the smallest free energy barrier and the least energetically expensive transition state. It is noteworthy that we do not find a direct minimum energy pathway between pairs of diagonal solutions or pairs of rotated solutions and all minimum energy pathways connect diagonal and rotated solutions. In the moderate anchoring regime, the transition pathways are not mediated by defects but are rather mediated by localized anchoring breaking along the edges, which induces a global transition between diagonal and rotated solutions. In the weak anchoring regime, the rotated solutions cease to be locally stable (consistent with the numerical findings in Luo et al. (2012)) but do exist as transition states connecting the diagonal solutions. In particular, the numerical methods in Tsakonas et al. (2007); Luo et al. (2012) cannot capture the persistence of rotated solutions, as non energy-minimizing solutions, in the weak anchoring regime.

We make some remarks on the novelty and importance of our approach. The transition states are not stable but they may be experimentally observable. Indeed the transition states in the strong anchoring limit are reminiscent of recent experimental observations of metastable states with point defects in shallow nematic chambers Lewis et al. (2014). Further, our numerical results demonstrate a subtle dependence of the optimal transition pathways on the anchoring strength. In the strong anchoring limit, the optimal transition pathway is almost independent of the anchoring strength whereas the free energy barrier decreases monotonically with anchoring strength in the moderate and weak anchoring regimes. Switching mechanisms rely on a complex interplay between anchoring strength, material properties and external fields and one could use our numerical methods to compute and analyze optimal switching pathways in the presence of external fields.

The paper is organized as follows. In Section II, we review the Landau-de Gennes theory for nematic liquid crystals. In Section III, we compute the Landau-de Gennes free energy minimizers and in Section IV, we compute the transition states and the transition pathways between free energy minima, focussing on three different anchoring regimes. We conclude in Section V with future perspectives. The computational methodologies used in this paper are described in the appendix.

## Ii Computational Model

Following the paradigm in Tsakonas et al. (2007); Luo et al. (2012), we model the planar bistable device within the two-dimensional Landau-de Gennes theory. In Landau-de Gennes theory, all the microscopic details about molecular shape and interactions are averaged out, and the nematic state is described by an macroscopic order parameter, . We take the computational domain, , to be a square in the -plane, which defines the bottom surface of the well,

(1) |

where denotes the system size. As in Luo et al. (2012), for two dimensions, the tensor can be represented by a symmetric traceless matrix,

(2) |

where is the director or the distinguished direction of molecular alignment and is a scalar order parameter that measures the degree of orientational ordering about .

We work with a simple form of the Landau-de Gennes energy as given below, with no external fields.

The first integral is the bulk free energy that drives the nematic-isotropic transition as a function of the temperature de Gennes and Prost (1995); Virga (1994). The coefficient is the re-scaled temperature; we work with temperatures below the critical nematic supercooling temperature and hence, ; the coefficients are positive material-dependent constants. For a two-dimensional as in (2), , and the bulk free energy simplifies to the familiar Ginzburg-Landau potential Luo et al. (2012). The second term is the one-constant elastic energy and is an elastic constant; there are more general quadratic elastic energy densities but we believe that the one-constant energy density suffices for qualitative purposes. The third term is a Durand-Nobili surface anchoring energy Nobili and Durand (1992) that enforces the preferred tangential anchoring on the square edges, with anchoring strength , reference configuration and boundary element . There are multiple choices for the surface energy (see Luo et al. (2012) for comparisons between three potential candidates for the tangent surface energy) but the Durand-Nobili energy has the desired numerical stability properties for all relevant ranges of .

To reduce the number of input parameters in the model, we now rewrite the free energy functional in its dimensionless form. Substituting Eq. (2) into (II) and defining , , , , , we obtain

Here, there are only two input parameters: the dimensionless elastic constant and the dimensionless surface anchoring strength . We drop the tildes in the subsequent text and all results are to be interpreted in terms of the dimensionless variables.

We need to prescribe a suitable form for the reference configuration, , that enforces the tangent boundary conditions on the edges. Following the formulation in Luo et al. Luo et al. (2012), we set , where , , and . The scalar has to vanish at the vertices with strong tangent anchoring, because of the mismatch in at the vertices. Hence, we define

with chosen to be . In the liquid crystal literature, is known to be typically proportional to defect core sizes de Gennes and Prost (1995); Virga (1994) and hence, we take to be approximately unity everywhere (a minimum of the bulk free energy) except for small neighbourhoods around the square vertices, where we expect to see defects.

We can use standard methods from calculus of variations to derive the Neumann boundary conditions on the edges (from the Durand-Nobili surface energy in Eq. (II)),

We use Greek indices to denote the components of the -tensor, and above.

The numerical methods for the computation of the free energy minima, transition states and minimum energy pathways are described in the appendix.

## Iii Minimum Free Energy States

We first study the minimum free energy morphologies of the system as the surface anchoring is varied. This is similar to the work reported in Luo et al. (2012) where the authors use continuation methods to compute stable equilibria of this system as is varied. Throughout this paper, we fix the value of the dimensionless elastic constant (for typical values of the elastic constant and Landau-de Gennes bulk potential parameters, this describes a micron-scale well) and use a square lattice with grid points.

In agreement with previous studies Luo et al. (2012); Lewis et al. (2014); Kusumaatmaja (2015), we find that for large anchoring strengths i.e. for , there are two classes of minima, namely diagonal and rotated states. We check the stability of the solutions by numerical computations of the Hessian of the Landau-de Gennes energy and its eigenvalues at a given solution, as is standard for numerical stability analysis. Due to the system symmetry, there are two equivalent diagonal states and four equivalent rotated states. These states are shown in Fig. 1.

For weak anchoring, , the rotated states are no longer stable minima. This is consistent with the bifurcation diagram in Luo et al. (2012) where the authors do not find rotated solutions for small . However, the rotated solutions do survive as transition states in the free energy landscape for weak anchoring, connecting the stable diagonal solutions. This is further explained in the next section and the two stable minima (diagonal states) are shown in Fig. 2.

## Iv Transition Pathways

Previous studies have been limited to free energy minima in the planar bistable device Tsakonas et al. (2007); Luo et al. (2012). In this paper we systematically compute the transition states, transition pathways and optimal transition pathways as a function of the surface anchoring parameter . To the best of our knowledge, this has not been reported elsewhere in the literature.

We distinguish three separate regimes depending on the values of : (i) strong anchoring regime, (ii) moderate anchoring regime and (iii) weak anchoring regime.

Regime I: Strong anchoring regime - The first regime is the strong anchoring limit, which we find to occur for . The optimal transition pathways are independent of in this regime whilst they are sensitive to for weaker anchoring. The computations in this regime are in fact very similar to the Dirichlet case considered in Kusumaatmaja (2015) where is fixed to be on the square edges.

We first make a few comments about the vertex defects in the diagonal and rotated solutions, both of which are local free energy minima in this regime. There are two types of point defects in these solutions: (i) splay defects wherein the director splays outward from the vertex and (ii) bend defects, for which the director bends between a pair of intersecting edges.

For a diagonal solution, there are two diagonally opposite splay defects and two diagonally opposite bend defects and for a rotated solution, the two splay defects are connected by an edge, as are the two bend defects. For a rotated to diagonal transition, we need two defect transformations along an edge: a bend to splay defect transformation and a splay to bend defect transformation via director rotation. We obtain at least three different transition pathways between a rotated and a diagonal solution, as illustrated in Figure 3 below, all of which are mediated by motion of point defects.

The first mechanism, as shown in Fig. 3(a-c), is mediated by a defect. The defect is created near the top left corner, propagates towards the bottom left corner and subsequently settles there to yield a bend defect. The numerical results suggest that the singular bend behaviour at the top left corner of the initial rotated solution is propelled into the square interior leading to the creation of a defect. As the defect propagates downwards, it rotates simultaneously and induces a global director rotation within the square domain. Finally, the point defect settles at the bottom left corner to become a bend defect and the system relaxes into a diagonal solution. The second mechanism, Fig. 3(d-f), is dominated by a defect created near the bottom left corner. The rotated solution has a splay defect at the bottom left corner; this defect is propelled into the square interior, it moves upwards along the square edge and whilst moving, rotates and induces a global director rotation. It finally settles at the top left corner as a splay defect and the system relaxes into a diagonal solution. These two mechanisms have identical free energy profiles, as shown in Figure 5. Further analysis using our numerical methods also shows that there can be other transition pathways with higher free energy barriers. We interpret the optimal transition pathway as being the minimum energy pathway with the smallest barrier and the least energetically expensive transition state. An example of a transition pathway with a higher free energy barrier is shown in Fig. 3 (g-i), where the pathway is mediated by the creation and propagation of a pair of defects. We also note that the transition pathways between other rotated and diagonal states can be related by symmetry to the specific example shown in Fig. 3. In Fig. 4, we plot the order parameter along the transition pathways, as shown in Fig. 3. The defects at the four corners, along with the bulk defects (mediating the transitions), are clearly identified by small order i.e. by small values for such near the defects.

The free energy profiles for the three transition mechanisms discussed in Fig. 3 are shown in Fig. 5(a). The -axis corresponds to the “normalised path length” between configuration along the pathway and the initial rotated state,

(5) |

The normalization constant is chosen such that the total path length between the rotated and diagonal states is 1. For the first two mechanisms described above, the profiles are virtually indistinguishable and the maxima in the free energy profiles correspond to the saddle points (transition states) in the free energy landscape of the system. The transition state configurations are depicted in Figs. 3(b), (e), and (h). It is worth noting that we do not find any direct pathway between two diagonal or two rotated states. A transition between two diagonal states or two rotated states is always composed of a sequence of rotated-diagonal transitions.

Regime II: Moderate Anchoring - We label the range as the medium anchoring regime. There are three key features of this regime: (i) the diagonal and rotated solutions survive as local free energy minima, (ii) the free energy barrier (along minimum energy pathways) decreases monotonically with and (iii) the optimal transition pathways do not feature defects.

This regime is further divided into two sub-regimes (a) and (b) . For case (a), the transition from a rotated to a diagonal solution is due to a “sequential” director rotation on one side of the square device, without any transient defects. Heuristically, it is energetically preferable to break the tangential anchoring along a square edge leading to director rotation along an edge as opposed to the creation of defects. As in the strong anchoring regime, we find that there are two mechanisms with degenerate free energy profiles (including the free energy barrier).

In the first mechanism, the director starts to rotate clockwise from the bottom left corner, breaking the tangential anchoring along the left vertical edge. The rotation then propagates upwards, as shown in Fig. 6(a-c) and this rotation suffices to transform the splay defect at the bottom left corner into a bend defect and conversely, the bend defect at the top left corner into a splay defect. The final state is a diagonal solution. The second mechanism is shown in Fig. 6(d-f). The director rotates clockwise from the top left corner and the sweeping motion propagates downwards, inducing a global director rotation and a transition from a rotated solution to a diagonal solution. The configurations in panels (b) and (e) correspond to the transition states and the two mechanisms have identical free energy profiles.

For weaker anchoring, for , the penalty for breaking tangential anchoring is weak and the rotated to diagonal transition follows from a “global” director rotation along the left vertical edge. As shown in Fig. 7(a-d), the director rotates simultaneously along the entire length of the left vertical edge and this rotation induces a global director distortion within the square, leading to a rotated to diagonal transition.

Regime III: Weak Anchoring - As the surface anchoring strength further decreases, for , the rotated solutions are no longer stable. The disappearance of the rotated states correspond to cusp catastrophes Wales (2001), where typically a minimum and two transition states coalesce to a single point. This point, reminiscent of the rotated state (as such, there are four of them), is a saddle point (transition state) between the two stable diagonal solutions. The transition state in panel Fig. 8 (c) resembles a rotated solution, with a uniform director profile in the square interior accompanied by two transition layers near the vertical edges. The corresponding free energy profile can be seen from Fig. 5(b) for . In this regime, the minimum energy pathways may therefore be represented by the diagram in Fig. 2 where the double arrow represents four equivalent transition pathways.

The movies for all the possible transition pathways described above are available as supporting information.

## V Conclusions

We have systematically studied the free energy landscape of a multistable nematic liquid crystal device as we vary the strength of the surface anchoring potential. Previous work focuses on the free energy minima and we take this work further by computing transition states, transition pathways and optimal transition pathways between minima, along with the free energy barriers.

We have classified the system behaviour into three categories, according to surface anchoring strength and qualitative properties of transition pathways. In the strong anchoring regime, , we obtain two distinct classes of free energy minima: diagonal and rotated solutions, as expected from previous work. The optimal transition pathways and the corresponding transition states feature either a or a defect, localized near an edge. It is worth noting that these transition states bear strong resemblance to experimental observations of states with internal defects reported in Lewis et al. (2014), suggesting possible metastability of transition states. Further, there are multiple transition pathways between free energy minima, which can be mediated by multiple defects but with higher free energy barriers than the optimal pathways.

The existence of a transition state or a saddle point-type critical point of the Landau-de Gennes energy (II), connecting the stable diagonal and rotated solutions in the strong anchoring regime, can be theoretically justified from the celebrated Mountain Pass Theorem Almeida (1999). In the limit of infinite anchoring, one can define a topological classification scheme for two-dimensional director fields (see Robbins and Zyskin (2004) for a similar three-dimensional topological classification), for which the diagonal and rotated solutions belong to different topological sectors. This can also be seen more intuitively from the arrangement of the splay and bend point defects in the diagonal and rotated solutions respectively. The Mountain Pass Theorem guarantees the existence of at least one critical point of the Landau-de Gennes energy connecting the topologically distinct diagonal and rotated solutions. In our framework, Mountain Pass solutions correspond to transition states along optimal transition pathways.

In the moderate anchoring regime, for , we do not observe any defect along the optimal transition pathways. The rotated-to-diagonal transition is driven by clockwise director rotation (and anti-clockwise for the opposite diagonal-to-rotated transition) localized along an edge. For , the rotation starts from one corner and propagates to another corner, along an edge. For weaker anchoring strengths (), it is preferable for the director to rotate along an entire edge, since the energy penalty for violating tangent anchoring is relatively weak.

For , we capture the cusp catastrophe events as the rotated states become unstable and act as transition states for the transition pathways connecting two stable diagonal solutions. This was previously not noted in the literature.

We note that the same effect i.e. loss of stability of the rotated solutions, can also be achieved by continuously decreasing the square size whilst maintaining strong anchoring on the square edges. In M.Robinson et al. (2015), the authors study the critical points of the two-dimensional Landau-de Gennes energy as a function of the square size, with strong anchoring on the edges. For micron-sized wells, they recover the diagonal and rotated solutions as we do. For a critical well size, the rotated solutions lose stability (but exist as non-minimizing solutions) whilst the diagonal solutions retain stability. For nano-scale wells, there is a unique critical point, referred to as the order-reconstruction solution, by analogy with similar findings in Kralj and Majumdar (2014).

For moderate and weak anchoring, the height of the free energy barrier decreases monotonically with whilst the height of the free energy barrier is independent of in the strong anchoring regime. Secondly, the configurations of the transition states and the free energy barriers depend strongly on the anchoring strength, for moderate and weak anchoring. These results can therefore be exploited to fine-tune the desirable surface properties i.e. use an anchoring strength that yields the optimal free energy barrier for device applications. For example, a very low barrier compromises the stability of the minimum free energy configurations and the device performance. A very high barrier is an impedance to realistic switching, necessitating large power input or strong external electric fields.

We will now discuss how the dimensionless free energy can be converted to SI unit (Joules). Note that the free energy defined in Eq.(II) is in two dimensions (i.e. free energy per unit thickness), while all practical experiments are, of course, in three dimensions. As such, we need to account for the thickness of the device, , such that , where is the dimensionless free energy defined in Eq. (II). The parameters and are experimentally measurable quantities. As an example, for the common liquid crystal material MBBA Mkaddem and Jr. (2000), the characteristic bulk constants are of order and . To estimate , we have assumed that the temperature is below the critical nematic transition temperature (C ). From Fig. 5, the typical free energy barrier in dimensionless units is of order . Assuming the thickness of the device is 1 micron and the area of the device is (100 micron), the typical free energy barrier in SI unit is therefore of order J. This is, of course, much larger than . On the other hand, if the device dimension can be miniaturized to 1 micron x 1 micron x 100 nm, then the typical free energy barrier becomes of order J, comparable to . The free energy barriers for other materials and device dimensions can be computed in an analogous way.

Our methods can be extended to three dimensions and to include external fields. For example, we could use similar methods to study transition pathways in three-dimensional set-ups, as studied in Kralj and Majumdar (2014), where the authors employ a three-dimensional Landau-de Gennes modelling approach and find a new biaxial order-reconstruction pattern for shallow nano-scale rectangular wells with strong anchoring. We speculate that the biaxial order-reconstruction pattern loses stability for larger wells but exists as a transition state connecting diagonal and rotated solutions for large micron-scale wells with strong anchoring.

Similarly, we can use these numerical methods to describe more complex geometries and boundary conditions, for example the Zenithally Bistable Device studied in Raisch and Majumdar (2014). They are also compatible with other simulation techniques tailored to the dynamics of liquid crystalline systems or generally complex fluids problems (such as finite element, lattice Boltzmann method, etc). These numerical methods promise to be a valuable design or optimization tool, allowing us to obtain a detailed picture of the free energy pathways of a complex system, and subsequently how the pathways evolves as the system parameters are varied.

## Vi Appendix: Computational Techniques

Here we summarise the most salient features of the computational techniques. A more detailed description and implementation of the methods can be found in Kusumaatmaja (2015).

### vi.1 Discretization of the Free Energy Functional

We discretize the Landau-de Gennes free energy on a square grid. Specific stencils are needed to approximate the derivatives and it is important that they are at least second order accurate. We have used

where the superscripts and label the lattice points in two dimensions. We have used a 150x150 grid based on the results of Luo et al. Luo et al. (2012), where they carried out convergence tests on the energy-minimizing solutions in the Landau-de Gennes framework.

Minimum Free Energy States - The first step in a survey of the free energy landscape is to compute the majority, if not all, of the possible minima in the system. We follow a stochastic approach, where each step consists of a trial move followed by an energy minimization. It is identical to a basin-hopping algorithm Li and Scheraga (1987); Wales and Doye (1997) at infinite temperature, such that every minimum state found is recorded. The simplest trial move consists of random perturbations of the lattice field values, . Here is a random number between -1.0 and 1.0, and is the amplitude of the perturbation, which we usually take to be . The energy minimization is carried out using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm Nocedal (1980); Liu and Nocedal (1989). In this paper, we typically take 500 (basin-hopping) steps to sample the minima of the system.

### vi.2 Transition States and Minimum Energy Pathways

Given the minimum free energy states, we numerically compute the transition states (saddle points in the free energy landscapes) using a combination of the doubly-nudged elastic band (DNEB) method Trygubenko and Wales (2004); Rühle et al. (2013) and the hybrid eigenvector-following technique Wales (2003). Transition states are special saddle points or critical points in the free energy landscapes where the energy gradients are zero in all eigendirections and one of the eigenvalues of the Hessian is negative. This is to be contrasted to minima, where the gradients are zero and the eigenvalues are positive in all eigendirections.

The DNEB method is a double-ended search algorithm for finding transition states and minimum energy pathways between any two pair of minima. Minimum energy pathways have an important feature such that every point on such a pathway is a minimum in all but a single distinguished direction in the phase space. As a consequence, a maximum along the minimum energy pathway corresponds to a transition state.

In DNEB, a set of images, , are placed between the two endpoints (minima), and . The symbol represents all the degrees of freedom in the system, i.e. . We typically use 30 images, and they are initialised by taking a linear interpolation between the two endpoints.

The images are relaxed using Landau-de Gennes energy gradients that include contributions from two components. The first contribution, , comes from the derivatives of the free energy functional with respect to the lattice degrees of freedom, . The second contribution is a spring force, , which is required to keep the images roughly equidistant. The following spring potential is applied on every image

(6) |

and are respectively the distances to the left and right images, . The coefficient is the spring constant.

Using the true and spring gradients, without any projection, often results in two issues Henkelman and Jónsson (2000, 2000): (i) corner-cutting, where images are pulled away from the minimum energy path; and (ii) sliding-down problems, where images slide down from barrier regions. In DNEB, we use two projections. First, we only retain the components of the true gradient that are perpendicular to the unit tangent vector ,

(7) |

Secondly, the following component of the spring constant gradient is retained,

(8) |

where . The unit tangent is defined as in Henkelman and Jónsson (2000). As a convergence criterion, we stop the DNEB run when the root mean square of the gradients as defined in Eqs. (7) and (8) in the manuscript is , which in our experience is more than enough to ensure the calculations have converged.

A maximum in the DNEB path corresponds to a transition state candidate, which we then refine by applying the hybrid eigenvector-following technique. We use a Rayleigh-Ritz approach Wales (2003), based on gradients of the Landau-de Gennes energy, to compute the transition states since this method is more efficient for large systems.

Once the transition states are found, small displacements are applied in the two downhill directions. Energy minimizations are again carried out using the LBFGS algorithm Nocedal (1980); Liu and Nocedal (1989). Given a minima - transition state - minima triplet, this run yields the transition pathway for this triplet and the corresponding free energy barrier. We check that this pathway is consistent with the one obtained from DNEB calculations. The free energy barriers could be interpreted as the free energy differences between the transition state and the free energy minima. There are typically multiple transition pathways between a pair of free energy minima. The optimal pathway is the pathway with the smallest free energy barrier.

## References

- de Gennes and Prost (1995) P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, New York, 1995.
- Lagerwall and Scalia (2012) J. P. F. Lagerwall and G. Scalia, Curr. Appl. Phys., 2012, 12, 1387–1412.
- Blancab et al. (2013) C. Blancab, D. Coursault and E. Lacaze, Liquid Crystals Reviews, 2013, 1, 83 –109.
- Tsakonas et al. (2007) C. Tsakonas, A. J. Davidson, C. V. Brown and N. J. Mottram, Appl. Phys. Lett., 2007, 90, 111913.
- Luo et al. (2012) C. Luo, A. Majumdar and R. Erban, Phys. Rev. E, 2012, 85, 061702.
- Kralj and Majumdar (2014) S. Kralj and A. Majumdar, Proc. Royal Soc. A, 2014, 470, 20140276.
- Anquetil-Deck and D.Cleaver (2010) C. Anquetil-Deck and D.Cleaver, Phys. Rev. E, 2010, 82, 031709.
- Anquetil-Deck et al. (2012) C. Anquetil-Deck, D. Cleaver and T. J. Atherton, Phys. Rev. E, 2012, 86, 041707.
- Lewis et al. (2014) A. Lewis, I. Garlea, J. Alvarado, O. Dammone, P. Howell, A. Majumdar, B. Mulder, M. P. Lettinga, G. Koenderink and D. Aarts, Soft Matter, 2014, 39, 7865–7873.
- Virga (1994) E. G. Virga, Variational Theories for Liquid Crystals, Chapman & Hall, London, 1994.
- Nobili and Durand (1992) M. Nobili and G. Durand, Phys. Rev. A, 1992, 46, R6174.
- Kusumaatmaja (2015) H. Kusumaatmaja, The Journal of Chemical Physics, 2015, 142, 124112.
- Wales (2001) D. J. Wales, Science, 2001, 293, 2067.
- Almeida (1999) L. Almeida, Revista Matemática Iberoamericana, 1999, 15, 487–545.
- Robbins and Zyskin (2004) J. M. Robbins and M. Zyskin, J. Phys. A, 2004, 37, 10609 – 10623.
- M.Robinson et al. (2015) M.Robinson, C.Luo, A.Majumdar and R.Erban, In preparation, 2015.
- Mkaddem and Jr. (2000) S. Mkaddem and E. C. G. Jr., Phys. Rev. E, 2000, 62, 6694–6705.
- Raisch and Majumdar (2014) A. Raisch and A. Majumdar, Europhys. Lett., 2014, 107, 16002.
- Li and Scheraga (1987) Z. Li and H. A. Scheraga, Proc. Natl. Acad. Sci. USA, 1987, 84, 6611.
- Wales and Doye (1997) D. J. Wales and J. P. K. Doye, J. Phys. Chem. A, 1997, 101, 5111.
- Nocedal (1980) J. Nocedal, Math. Comput., 1980, 35, 773.
- Liu and Nocedal (1989) D. Liu and J. Nocedal, Math. Program., 1989, 45, 503.
- Trygubenko and Wales (2004) S. A. Trygubenko and D. J. Wales, J. Chem. Phys., 2004, 120, 2082–2094.
- Rühle et al. (2013) V. Rühle, H. Kusumaatmaja, D. Chakrabarti and D. J. Wales, J. Chem. Theory Comput., 2013, 9, 4026–4034.
- Wales (2003) D. J. Wales, Energy Landscapes, Cambridge University Press, Cambridge, 2003.
- Henkelman and Jónsson (2000) G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
- Henkelman and Jónsson (2000) G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.