Solution Landscapes in the Landau-de Gennes Theory on Rectangles

TL;DR

This paper analyzes nematic liquid crystal equilibria on rectangles within the Landau-de Gennes theory, revealing how defect structures vary with domain size and shape, and exploring topological transitions.

Contribution

It provides exact limiting profiles for nematic equilibria in different size regimes and introduces the concept of topology relaxation via defects in rectangular domains.

Findings

Line defects near edges in large domains

Fractional point defects in small domains

Bifurcation diagrams as a function of size and shape

Abstract

We study nematic equilibria on rectangular domains, in a reduced two-dimensional Landau-de Gennes framework. These reduced equilibria carry over to the three-dimensional framework at a special temperature. There is one essential model variable---$\epsilon$ which is a geometry-dependent and material-dependent variable. We compute the limiting profiles exactly in two distinguished limits---the $\epsilon \to 0$ limit relevant for macroscopic domains and the $\epsilon \to \infty$ limit relevant for nano-scale domains. The limiting profile has line defects near the shorter edges in the $\epsilon \to \infty$ limit whereas we observe fractional point defects in the $\epsilon \to 0$ limit. The analytical studies are complemented by bifurcation diagrams for these reduced equilibria as a function of $\epsilon$ and the rectangular aspect ratio. We also introduce the concept of `non-trivial' topologies and the relaxation of non-trivial topologies to trivial topologies mediated via point and line defects, with potential consequences for non-equilibrium phenomena and switching dynamics.

Tables

Table 1. Table 3.1 : Boundary Condition of solutions

state	shape	$d_1$	$d_2$	$d_3$	$d_4$
D1	$╱$	$0$	$+\pi /2$	$0$	$+\pi /2$
D2	$╲$	$0$	$-\pi /2$	$0$	$-\pi /2$
R1	$\subset$	$0$	$-\pi /2$	$-\pi$	$-\pi /2$
R2	$\supset$	$0$	$+\pi /2$	$+\pi$	$+\pi /2$
R3	$\cap$	$0$	$-\pi /2$	$0$	$+\pi /2$
R4	$\cup$	$0$	$+\pi /2$	$0$	$-\pi /2$

Equations

$${\mathbf{Q}}_{\text{LdG}}=\sum _{i=1}^3{\lambda }_i{\hat{\mathbf{e}}}_i\otimes {\hat{\mathbf{e}}}_i,$$

$$\mathcal{B}=\left\{\left(x,y,z\right)\in {\mathbb{R}}^3\mid \left(x,y\right)\in {\Omega }_{aL,bL},0\le z\le h\right\},$$

$${\Omega }_{p,q}:=\left\{\left(x,y\right)\in {\mathbb{R}}^2\mid 0\le x\le p,0\le y\le q\right\}.$$

$$F_{\text{LdG}}\left[{\mathbf{Q}}_{\text{LdG}}\right]:={\int }_{\mathcal{B}}\left(\frac{K}{2}{\left|\nabla {\mathbf{Q}}_{\text{LdG}}\right|}^2+f_{\text{B}}\left({\mathbf{Q}}_{\text{LdG}}\right)\right)+{\int }_{\partial \mathcal{B}}{f}_{\text{S}}\left({\mathbf{Q}}_{\text{LdG}}\right),$$

$${\left|\nabla {\mathbf{Q}}_{\text{LdG}}\right|}^2$$

$$f_{\text{B}}\left({\mathbf{Q}}_{\text{LdG}}\right)$$

$$\displaystyle\mathbf{Q}_{\textrm{LdG}}=\left(\begin{array}[]{ccc}q_{1}-q_{3}&q_{2}&0\\ q_{2}&-q_{1}-q_{3}&0\\ 0&0&2q_{3}\end{array}\right),$$

$$\displaystyle\mathbf{Q}_{\textrm{LdG}}=\left(\begin{array}[]{cc}\mathbf{Q}&\\ &0\end{array}\right)-\frac{B}{6C}\left(\begin{array}[]{ccc}-1&&\\ &-1&\\ &&2\end{array}\right)$$

$$\frac{1}{hK}F\left[\bar{\mathbf{Q}}\right]:={\int }_{\tilde{\Omega }}\left(\frac{1}{2}{\left|\nabla \bar{\mathbf{Q}}\right|}^2+\frac{L^2}{K}\left(-\frac{B^2}{4C}\mathrm{Tr}\left({\bar{\mathbf{Q}}}^2\right)+\frac{C}{4}{\left(\mathrm{Tr}\left({\bar{\mathbf{Q}}}^2\right)\right)}^2\right)-\frac{5}{432}\frac{B^4{L}^2}{C^{3}K}\right),$$

$$\tilde{\Omega }:={\Omega }_{a,b}=\left[0,a\right]\times \left[0,b\right].$$

$$\bar{\mathbf{Q}}\left(x,y\right)=2s\left(x,y\right)\left(\hat{\mathbf{n}}\left(x,y\right)\otimes \hat{\mathbf{n}}\left(x,y\right)-{\mathbf{I}}_2/2\right)$$

$$\displaystyle\bar{\mathbf{Q}}\left(x,y\right)=\left(\begin{array}[]{cc}\bar{Q}_{11}\left(x,y\right)&\bar{Q}_{12}\left(x,y\right)\\ \bar{Q}_{12}\left(x,y\right)&-\bar{Q}_{11}\left(x,y\right)\end{array}\right),$$

$$\displaystyle\operatorname{atan2}2\left(y,x\right)=\left\{\begin{array}[]{ll}\arctan\left(y/x\right),&x>0,\\ \arctan\left(y/x\right)+\pi,&x\leq 0,y\geq 0,\\ \arctan\left(y/x\right)-\pi,&x\leq 0,y<0.\end{array}\right.$$

$$\hat{\mathbf{n}}\left(x,y\right)$$

$${\bar{Q}}_{11}\left(x,y\right)$$

$${\bar{Q}}_{12}\left(x,y\right)$$

$$\tilde{\mathbf{Q}}\left(x,y\right):=\bar{\mathbf{Q}}\left(x,y\right)/s_0,\left(x,y\right)\in \tilde{\Omega },$$

$$ϵ:=\frac{\sqrt{2CK}}{BL}.$$

$$\frac{4C^2}{B^{2}hK}\tilde{F}\left[\tilde{\mathbf{Q}}\right]:={\int }_{\tilde{\Omega }}\left({\left|\nabla {\tilde{Q}}_{11}\right|}^2+{\left|\nabla {\tilde{Q}}_{12}\right|}^2+\frac{1}{{ϵ}^2}{\left({\tilde{Q}}_{11}^2+{\tilde{Q}}_{12}^2-1\right)}^2-\frac{32}{27}\frac{1}{{ϵ}^2}\right),$$

$$\displaystyle\left\{\begin{array}[]{c}\Delta\tilde{Q}_{11}=\epsilon^{-2}\left(\tilde{Q}_{11}^{2}+\tilde{Q}_{12}^{2}-1\right)\tilde{Q}_{11},\\ \Delta\tilde{Q}_{12}=\epsilon^{-2}\left(\tilde{Q}_{11}^{2}+\tilde{Q}_{12}^{2}-1\right)\tilde{Q}_{12}.\end{array}\right.$$

$$\displaystyle\tilde{\mathbf{g}}_{d;a,b}\left(x,y\right):=\left\{\begin{array}[]{cc}\left[+T_{\frac{d}{a}}\left(\frac{x}{a}\right),0\right],&x\in\left[0,a\right],y\in\left\{0,b\right\},\\ \left[-T_{\frac{d}{b}}\left(\frac{y}{b}\right),0\right],&y\in\left[0,b\right],x\in\left\{0,a\right\}.\end{array}\right.$$

$$\displaystyle T_{d}\left(t\right)=\min\left\{t/d,1,\left(1-t\right)/d\right\}=\left\{\begin{array}[]{cc}t/d,&0\leq t\leq d,\\ 1,&d\leq t\leq 1-d,\\ \left(1-t\right)/d,&1-d\leq t\leq 1.\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{c}\Delta\tilde{Q}_{11}=0,\\ \Delta\tilde{Q}_{12}=0.\end{array}\right.$$

$$\forall i=1,2,{\Vert {\left({\tilde{Q}}_{ϵ}\right)}_{1i}-{\left({\tilde{Q}}_{\infty }\right)}_{1i}\Vert }_{L^{\infty }\left(\tilde{\Omega }\right)}\le C{ϵ}^{-2},$$

$$\forall i=1,2,$$

$$\Delta v=1,$$

$$v=0,$$

$$\forall i=1,2,{\Vert {\left({\tilde{Q}}_{ϵ}\right)}_{1i}-{\left({\tilde{Q}}_{\infty }\right)}_{1i}\Vert }_{L^{\infty }\left(\tilde{\Omega }\right)}\sim \mathcal{O}\left({ϵ}^{-2}\right).$$

$${\left({\tilde{Q}}_{\infty }\right)}_{11}\left(x,y\right)$$

$$=\sum _{k\mathrm{odd}}\frac{4\sin \left(k\pi d/a\right)}{k^2{\pi }^{2}d/a}\sin \left(\frac{k\pi x}{a}\right)\frac{\sinh \left(k\pi \left(b-y\right)/a\right)+\sinh \left(k\pi y/a\right)}{\sinh \left(k\pi b/a\right)}$$

$$-\sum _{k\mathrm{odd}}\frac{4\sin \left(k\pi d/b\right)}{k^2{\pi }^{2}d/b}\sin \left(\frac{k\pi y}{b}\right)\frac{\sinh \left(k\pi \left(a-x\right)/b\right)+\sinh \left(k\pi x/b\right)}{\sinh \left(k\pi a/b\right)}.$$

$$\Delta {\left({\tilde{Q}}_{\infty }\right)}_{11}\left(x,y\right)$$

$${\left({\tilde{Q}}_{\infty }\right)}_{11}\left(x,0\right)={\left({\tilde{Q}}_{\infty }\right)}_{11}\left(x,b\right)$$

Full text

1. 1 Introduction

2. 2 Theoretical Framework

3. 3 Two Limiting Problems in terms of $\epsilon$

4. 3.1 The $\epsilon\to\infty$ Limit

5. 3.1.1 Strong Anchoring Condition

6. 3.1.2 Weak Anchoring Condition

7. 3.2 The $\epsilon\to 0$ Limiting Problem

8. 4 Bifurcation Diagrams as a Function of $\epsilon$ and Rectangular Anisotropy

9. 5 Relaxation Mechanisms for Non-trivial Topologies

10. 6 Conclusion

11. 7 Acknowledgement

1 Introduction

Nematic liquid crystals are quintessential examples of partially ordered materials intermediate between solids and liquids [8] [31]. The constituent nematic molecules are typically asymmetric in shape and the macroscopic nematic state has a degree of long-range orientational ordering, i.e., the constituent molecules tend to align along locally preferred directions referred to as ‘directors’ in the literature. The directional nature of nematics makes them highly responsive materials and indeed, their sensitivity to light and external electric fields have made nematics the working material of choice for the multi-billion dollar liquid crystal display industry [1]. There is huge contemporary interest in understanding pattern formation in nematic systems, in the context of micro-patterned systems, thin geometries including surfaces and interfacial phenomena [11], [4], [25], [24].

There are several continuum macroscopic theories for nematic liquid crystals. We focus on the celebrated Landau-de Gennes (LdG) theory for nematics, in a reduced two-dimensional (2D) setting, as described in section 2. This approach is accepted for very thin domains, which are treated as two-dimensional domains and the nematic directors are assumed to be in the plane of the domain with tangent or in-plane boundary conditions. This reduction can be rigorously justified (see [13]) and for a special temperature, we can make an immediate connection between the reduced 2D approach and the full three-dimensional (3D) approach, as explained in detail in the next section. The 2D LdG theory is essentially the Ginzburg-Landau theory [3], phrased in terms of the 2D LdG $\mathbf{Q}$-tensor, which is a symmetric and traceless $2\times 2$ matrix with two degrees of freedom. We refer to the eigenvector with the largest positive eigenvalue as being the nematic director in this framework and the order parameter, $s$, is proportional to the corresponding eigenvalue. The corresponding equilibria are local or global minimizers of a reduced 2D LdG free energy, which is effectively the Ginzburg-Landau free energy with one rescaled parameter, $\epsilon$. This parameter encodes the geometric length scale of the domain, material properties and the nematic elasticity; more precisely at the fixed temperature under consideration, it can be interpreted as the ratio of the nematic correlation length to the macroscopic domain size [17]. The nematic correlation length is a characteristic material and temperature-dependent length scale related to the size of nematic defects [19].

The prototype problem of nematics inside square domains has been studied extensively [23], [18], [22], [17], following the experimental and modelling work presented in [32]. The reduced 2D LdG equilibria on square domains with tangent boundary conditions, which require the nematic director to be tangent or parallel to the square edges, have been well classified. In particular, in [17], the authors report a new Well Order reconstruction Solution (WORS) in this reduced setting, distinguished by a pair of mutually orthogonal defect lines along the square diagonals. The WORS is the unique nematic equilibrium for sufficiently large values of $\epsilon$. For small $\epsilon$, there are at least six different competing equilibria : two ‘diagonal’ states for which the director is aligned along one of the square diagonals and four ‘rotated’ states, for which the director rotates by $\pi$ radians between a pair of opposite edges. It is very natural to study the counterpart of this problem on rectangles with tangent boundary conditions, to systematically study the effects of geometrical anisotropy.

It is difficult to make analytic progress for the problem in general, except for two distinguished limits: the $\epsilon\to\infty$ limit which describes nano-scale geometries or ‘small’ geometries, and the $\epsilon\to 0$ limit which describes macroscopic domains with size much greater than the nematic correlation length. In the $\epsilon\to\infty$ limit, the problem effectively reduces to a Dirichlet boundary-value problem for both the components of the limiting profile, $\mathbf{Q}_{\infty}$. We show that the cross structure of the WORS is lost on all rectangles. However, the cross structure of the WORS is retained on square domains, if we replace the Dirichlet conditions with surface energies that enforce tangent conditions on the square edges, for large values of the surface anchoring parameter, $W\propto\epsilon$ for large values of $\epsilon$. In other words, the cross structure is specific to a square domain but not an artefact of Dirichlet conditions. On a rectangle, the limiting profile, $\mathbf{Q}_{\infty}$ has two nodal lines or defect lines along the shorter rectangular edges, referred to as sBD2 solutions in the remainder of the manuscript. The sBD2 solutions also survive on rectangles with surface energies and large values of the surface anchoring parameter i.e. the surface energies do not change the qualitative conclusions for large values of $W$. We use a combination of formal calculations and elegant maximum principle arguments to illustrate the effects of geometrical anisotropy in this limit.

In the $\epsilon\to 0$ limit, the limiting problem is just the Laplace equation for an angle $\theta$ in the plane of the rectangle, subject to Dirichlet conditions on the rectangular edges. In this limit, $s$ is a constant away from the rectangular vertices and the director, $\hat{\mathbf{n}}=\left(\cos\theta,\sin\theta\right)$. The tangent conditions necessarily create discontinuities in $\theta$ at the vertices and the limiting energy is effectively the Dirichlet energy of $\theta$ away from the vertices (also see [28], [10] and [2]). There are three competing nematic equilibria in this limit—the diagonal solutions, the rotated states (R1, R2) for which the director rotates between a pair of horizontal edges and the rotated (R3, R4) solutions for which the director rotates between a pair of vertical edges. The four rotated states are not energetically degenerate on a rectangle, in contrast to a square. We use simple symmetry-based arguments to compare the respective energies of these solutions and these arguments may be of wider interest. Using standard arc-continuation methods, we numerically investigate the solution landscapes on rectangles as a function of $\epsilon$ and geometrical anisotropy. The bifurcation diagrams are qualitatively similar to those reported for a square (see [30]) except for the fact that the energetically expensive R1, R2 solutions are disconnected from the connected diagonal and R3, R4 solution branches. The geometrical anisotropy affects R1, R2 solution branches, in the sense that they cannot be numerically continued to larger values of $\epsilon$ and our numerical methods cannot track the pathway from the R1, R2 solutions to the unique sBD2 solution as $\epsilon$ increases. This might have interesting experimental consequences, particularly for transition states that connect the R1, R2 stable states to the diagonal or R3, R4 states. We will investigate this in future work.

We introduce the concept of non-trivial topologies, based on the profile of $\theta$ near the rectangular vertices. This is largely a theoretical notion, inspired by super-twisted nematic devices in a planar setting [1]. The director is constrained on the rectangular edges by virtue of the tangent conditions but is free to rotate in the plane, between the edges. The minimal allowed rotation is labelled as ‘trivial’ and excess rotation is labelled as ‘non-trivial’. We speculate that non-trivial topologies can be realised by localised rotation of the nematic molecules near the vertices of a rectangular domain. Analogous remarks apply to shallow 3D wells with a rectangular cross-section. These non-trivial topologies will relax to trivial topologies when the rotation is removed. The relaxation pathways are mediated by defects of different dimensionality, depending on $\epsilon$ and the geometrical anisotropy, and can offer new optical possibilities or switching mechanisms between trivial states e.g. diagonal, R3 and R4 states. The relaxation pathways will be even more interesting with external electric fields, which are not accounted for in this manuscript.

To summarise, our paper can be regarded as a set of results for reduced LdG equilibria on rectangular domains with tangent boundary conditions. These results highlight the effects of geometrical anisotropy and to some extent, boundary conditions in certain limiting cases. The interplay between $\epsilon$ and the geometrical anisotropy has not come across yet. Such studies do not necessarily have unexpected outcomes but are much needed for systematic and structured investigations of the effects of geometry on nematic equilibria e.g. regular polygons versus polygons with sides of unequal length, dimensions of defect sets and control of defect sets and multiplicity of equilibria by tuning the geometrical parameters.

The paper is organised as follows. In section 2, we review the reduced 2D LdG theory, the free energy, the rescaling, the governing partial differential equations, and the boundary conditions. In section 3, we analytically study distinguished asymptotic limits on a rectangular domain with tangent boundary conditions. In section 4, we numerically study how the 2D LdG equilibria depend on $\epsilon$, for different values of the rectangular anisotropy. In section 5, we discuss relaxation mechanisms and their dependence on $\epsilon$ and we conclude in section 6 with some perspectives.

2 Theoretical Framework

We work within the continuum Landau-de Gennes (LdG) theory for nematic liquid crystals. The LdG theory is one of the most powerful continuum theories for nematics in the literature [8]. It describes the state of nematic anisotropy by the LdG $\mathbf{Q}_{\textrm{LdG}}$-tensor order parameter—a symmetric, traceless $3\times 3$ matrix whose eigenvectors model the directions of averaged molecular alignment in space (referred to as directors) and the eigenvalues are a measure of the degree of orientational order about the corresponding eigenvectors. From the spectral decomposition theorem, we can write $\mathbf{Q}_{\textrm{LdG}}$ as

[TABLE]

where $\left\{\hat{\mathbf{e}}_{1},\hat{\mathbf{e}}_{2},\hat{\mathbf{e}}_{3}\right\}$ constitute an orthonormal basis and $\sum_{i=1}^{3}\lambda_{i}=0$. A $\mathbf{Q}_{\textrm{LdG}}$-tensor is said to be (i) isotropic if $\mathbf{Q}_{\textrm{LdG}}=0$, (ii) uniaxial if $\mathbf{Q}_{\textrm{LdG}}$ has a pair of degenerate non-zero eigenvalues and (iii) biaxial if $\mathbf{Q}_{\textrm{LdG}}$ has three distinct eigenvalues [8]. We refer to the ‘director’ as being the eigenvector with the largest positive eigenvalue.

We study nematic equilibria on three-dimensional wells with a rectangular cross section, i.e., our domain is

[TABLE]

where $h>0$ is the height of the well, and for $p,q>0$,

[TABLE]

Here, $L>0$ is a fixed length, $a$ and $b$ measure the geometrical anisotropy and in what follows, we always fix $b=1$ and vary $a$. For the sake of simplicity, we will drop the subscripts in $\Omega_{aL,bL}$.

We impose tangent boundary conditions on the well surfaces, in the form of Dirichlet uniaxial conditions on the lateral surfaces and surface energies on the top and bottom surfaces, that enforce planar degenerate anchoring, i.e., constrain the director to be in the plane of the surface whilst respecting the tangent conditions on the edges of the rectangular cross-section. Our work builds on the substantial work done for square domains [32], [23], [17], [30], [34] and more recently on rectangular domains [36].

A particularly simple form of the LdG free energy is given by

[TABLE]

where $K>0$ is a material-dependent elastic constant, $f_{\textrm{B}}$ is the bulk energy, $f_{\textrm{S}}$ is the surface energy, and

[TABLE]

The variable $A=\alpha\left(T-T^{*}\right)$ is a rescaled temperature, $\alpha,B,C>0$ are material-dependent constants, and $T^{*}$ is the characteristic nematic supercooling temperature. Further, $\mathrm{Tr}\left({\mathbf{Q}_{\textrm{LdG}}^{2}}\right)=\sum_{i,j=1}^{3}\left(Q_{\textrm{LdG}}\right)_{ij}\left(Q_{\textrm{LdG}}\right)_{ji}$ and $\mathrm{Tr}\left({\mathbf{Q}_{\textrm{LdG}}^{3}}\right)=\sum_{i,j,k=1}^{3}\left(Q_{\textrm{LdG}}\right)_{ij}\left(Q_{\textrm{LdG}}\right)_{jk}\left(Q_{\textrm{LdG}}\right)_{ki}$. The rescaled temperature $A$ has three characteristic values [26]: (i) $A=0$, below which the isotropic phase $\mathbf{Q}_{\textrm{LdG}}\equiv 0$ loses stability, (ii) the nematic-isotropic transition temperature, $A=\frac{B^{2}}{27C}$, at which $f_{\textrm{B}}$ is minimized by the isotropic phase and a continuum of uniaxial states, and (iii) the nematic superheating temperature, $A=\frac{B^{2}}{24C}$, above which the isotropic state is the unique critical point of $f_{\textrm{B}}$. We work with $A<0$ throughout this manuscript for which $f_{\textrm{B}}$ favours an ordered uniaxial nematic phase and the physically relevant/observable nematic equilibria correspond to either global or local LdG minimizers subject to the imposed boundary conditions.

Following the methods in [12] and for certain physically relevant choices of the surface anchoring energies (that favour planar degenerate anchoring on the top and bottom surfaces) on the top and bottom surfaces, we can rigorously justify the reduction from the three-dimensional domain $\mathcal{B}$ to the two-dimensional domain $\Omega$ , in the $h\to 0$ limit or the thin film limit. In this limit, the physically relevant equilibria belong to the class of $\mathbf{Q}_{\textrm{LdG}}$-tensors that have the unit vector in the $z$-direction as a fixed eigenvector so that $\mathbf{Q}_{\textrm{LdG}}$ can be written as

[TABLE]

where the three degrees of freedom of $\mathbf{Q}_{\textrm{LdG}}$ are contained in the three scalar functions—$q_{1}$, $q_{2}$ and $q_{3}$ which only depend on the spatial variables in the cross-section. The details of this reduction can be found in [34]. In this limit, it suffices to study minimizers of the LdG free energy on the two-dimensional cross-section $\Omega$ in the class of $\mathbf{Q}_{\textrm{LdG}}$-tensors in Equation 2.10. Moreover, in [6], the authors compute rigorous bounds for $q_{3}$ as a function of $A$, building on previous work [15] in a different context. They show that for a particular choice $A=-\frac{B^{2}}{3C}$, all competing equilibria have constant $q_{3}\equiv-\frac{B}{6C}$ so that the physically relevant $\mathbf{Q}_{\textrm{LdG}}$-tensors can be written as

[TABLE]

with only two degrees of freedom. Here, $\mathbf{Q}$ is a symmetric traceless $2\times 2$ matrix with only two degrees of freedom that account for the nematic ordering in the plane of $\Omega$. Substituting Equation 2.16 into the LdG free energy and using Gamma-convergence methods in the $h\to 0$ limit as in [12], one can show that for $A=-\frac{B^{2}}{3C}$, the LdG free energy minimizers converge uniformly to minimizers of the reduced energy (after suitable rescaling $\bar{\mathbf{Q}}\left(x,y\right)=\mathbf{Q}\left(xL,yL\right)$)

[TABLE]

and

[TABLE]

In Equation 2.17, we have reduced the original three-dimensional problem to a two-dimensional problem on the re-scaled rectangle $\tilde{\Omega}$ at the fixed temperature $A=-\frac{B^{2}}{3C}$ and the remainder of this manuscript is devoted to this reduced problem; generalizations will be discussed in the Conclusions. We can equivalently write $\bar{\mathbf{Q}}$ in terms of a scalar order parameter field $s:\tilde{\Omega}\to\mathbb{R}$ and a nematic director $\hat{\mathbf{n}}:\tilde{\Omega}\to\mathbb{S}^{1}$, which is a unit-vector field that models the preferred direction of nematic alignment in the plane of $\tilde{\Omega}$ as shown below. For $\left(x,y\right)\in\tilde{\Omega}$,

[TABLE]

Alternatively,

[TABLE]

where we can relate the matrix components $\bar{Q}_{11}$ and $\bar{Q}_{12}$ to $s$ and $\hat{\mathbf{n}}$ by defining a new variable $\theta:=\operatorname{atan2}2\left(\bar{Q}_{12},\bar{Q}_{11}\right)/2$, where,

[TABLE]

When $\bar{\mathbf{Q}}\neq 0$, $\theta$ is well-defined, and we have the following important relations,

[TABLE]

We rescale the reduced energy Equation 2.17 by defining $s_{0}:=\frac{B}{2C}$ (so that $\mathrm{Tr}\left({\bar{\mathbf{Q}}^{2}}\right)=2s_{0}^{2}$ for minimizers of the bulk potential in Equation 2.17) and $\tilde{\mathbf{Q}}$ as

[TABLE]

along with

[TABLE]

Since we have fixed $A=-\frac{B^{2}}{3C}$, the ratio $\xi:=\sqrt{\frac{KC}{B^{2}}}$ is proportional to the nematic correlation length [33], which is a characteristic material-dependent length scale. Therefore, $\epsilon$ is the ratio of the nematic correlation length to the domain size (measured in terms of $L$). The rescaled energy is then equivalent to

[TABLE]

and all critical points (including global and local energy minimizers) are classical solutions of the associated Euler-Lagrange equations, for smooth boundary conditions. The Euler-Lagrange equations are given by the following system of elliptic coupled partial differential equations, also referred to as the Ginzburg-Landau equations on $\tilde{\Omega}$ in the literature [2].

[TABLE]

The last ingredient is the choice of boundary conditions. As in [32], [23], we work with tangent boundary conditions so that $\theta$ is a multiple of $\pi$ on the horizontal edges (defined by $y=0$ and $y=b$) and $\theta$ is an odd multiple of $\pi/2$ on the vertical edges (defined by $x=0$ and $x=a$). Recalling the relations Equation 2.27, this implies that $\tilde{Q}_{12}=0$ on all four edges, $\tilde{Q}_{11}=+1$ on the horizontal edges and $\tilde{Q}_{11}=-1$ on the vertical edges. However, there is a necessary mismatch at the corners; we adopt the same interpolatory approach as in [23] and define a vector field $\tilde{\mathbf{g}}_{d}$ for $0<d<\min\left\{a/2,b/2\right\}$ as,

[TABLE]

The trapezoidal function $T_{d}:\left[0,1\right]\to\mathbb{R}$ is given by,

[TABLE]

For the sake of simplicity, we will use the brief notation $\tilde{\mathbf{g}}_{d}$ instead of $\tilde{\mathbf{g}}_{d;a,b}$.

We then fix $\left(\tilde{Q}_{11},\tilde{Q}_{12}\right)=\tilde{\mathbf{g}}_{d}$ on $\partial\tilde{\Omega}$. We do not rigorously justify this choice of boundary conditions but we believe that the qualitative trends are not affected by the choice of $d$ above, provided it is sufficiently small. In what follows, we analytically and numerically study the reduced and rescaled tensor $\tilde{\mathbf{Q}}\in W^{1,2}\left(\tilde{\Omega};S_{0}\right)$ subject to the Dirichlet condition, where $S_{0}$ is the space of symmetric, traceless $2\times 2$ matrices and $W^{1,2}$ is the usual Sobolev space [9].

3 Two Limiting Problems in terms of $\epsilon$

We can make analytic progress with the solutions of the system Equation 2.35 in two distinguished limits, the $\epsilon\to\infty$ and $\epsilon\to 0$ limits respectively. The two limits have different physical interpretations. For $A=-\frac{B^{2}}{3C}$, the length scale $\xi=\sqrt{\frac{KC}{B^{2}}}$ is proportional to the nematic correlation length, which is typically of the order of tens to hundreds of nanometres and describes small-scale phenomena [33]. Hence, the first limit, $\epsilon\to\infty$, is relevant for domains with characteristic length comparable to the nematic correlation length, e.g., few hundred nanometres [17], [7]. The second limit, $\epsilon\to 0$, is relevant for domains with characteristic length $L$ much larger than the nematic correlation length, e.g., $L$ is on the micron-scale which is the usual experimentally achievable length scale although nano-scale geometries are becoming a reality with advanced nano-fabrication techniques.

3.1 The $\epsilon\to\infty$ Limit

3.1.1 Strong Anchoring Condition

In this section, we discuss solutions of Equation 2.35 in the $\epsilon\to\infty$ limit, with the Dirichlet condition $\left(\tilde{Q}_{11},\tilde{Q}_{12}\right)=\tilde{\mathbf{g}}_{d}$ defined in Equation 2.38. Dirichlet conditions are referred to as ‘strong anchoring’ in the literature. In practice, it is experimentally difficult to realise Dirichlet conditions but Dirichlet problems are relatively analytically tractable and are often, a good approximation to realistic scenarios.

One can prove there exists a critical $\epsilon_{c}>0$ such that Equation 2.35 has a unique solution for $\epsilon>\epsilon_{c}$ [20], [7]. The limit of large $\epsilon$ with strong anchoring has been studied extensively on square domains, where the authors of [17] and [7] report the new Well Order Reconstruction Solution (WORS) in this limit. The reduced WORS $\bar{\mathbf{Q}}$-tensor in Equation 2.17 vanishes along the square diagonals and the nodal lines of $\bar{\mathbf{Q}}$ partition the square into four quadrants such that the director $\hat{\mathbf{n}}$ is constant in each quadrant. In [7], the authors study the WORS in terms of solutions of the scalar Allen-Cahn equation at the special temperature $A=-\frac{B^{2}}{3C}$ and rigorously prove that the WORS exists for all $\epsilon>0$, is globally stable for $\epsilon$ large enough and loses stability as $\epsilon$ decreases.

A natural question is—is the WORS specific to a square domain with strong anchoring conditions in the $\epsilon\to\infty$ limit? We address this question in two parts: (i) by studying solutions of Equation 2.35 on a rectangle with $a>b=1$ in $\Omega$ and strong anchoring conditions and (ii) by studying solutions of Equation 2.35, in the $\epsilon\to\infty$ limit, on a square domain with weak anchoring, i.e., a weaker implementation of the Dirichlet condition Equation 2.38 as described below.

In the strong anchoring case, it is relatively straightforward to prove that $\left\|\tilde{\mathbf{Q}}\right\|_{L^{\infty}\left(\tilde{\Omega}\right)}$ is bounded independently of $\epsilon$ for all solutions of Equation 2.35 by maximum principle arguments as in [26] (here $L^{\infty}$ is the usual norm i.e. $\left\|\tilde{\mathbf{Q}}\right\|_{L^{\infty}\left(\tilde{\Omega}\right)}=\operatorname{esssup}_{\mathbf{r}\in\tilde{\Omega}}\left\|\tilde{\mathbf{Q}}\left(\mathbf{r}\right)\right\|$). Hence, in the limit $\epsilon\to\infty$, the coupled system Equation 2.35 reduces to the uncoupled Laplace equations on $\tilde{\Omega}$ as given below,

[TABLE]

Our first result shows that the unique solution of Equation 2.35, with strong anchoring and $\epsilon$ sufficiently large, is well approximated by the unique solution of the limiting problem Equation 3.3.

Proposition 3.1**.**

Let $\tilde{\mathbf{Q}}_{\epsilon}\in W^{1,2}\left(\tilde{\Omega};S_{0}\right)$ be a solution of the LdG Euler-Lagrange equations Equation 2.35 on the rectangular domain $\tilde{\Omega}$ for $\epsilon>0$, subject to the boundary condition $\left(\left(\tilde{Q}_{\epsilon}\right)_{11},\left(\tilde{Q}_{\epsilon}\right)_{12}\right)=\tilde{\mathbf{g}}_{d}$ on $\partial\tilde{\Omega}$ for a given $0<d<\min\left\{a/4,b/4\right\}$. Then $\tilde{\mathbf{Q}}_{\epsilon}$ converges to the unique solution $\tilde{\mathbf{Q}}_{\infty}$ of Equation 3.3 as $\epsilon\to\infty$, subject to the same Dirichlet condition with error estimates

[TABLE]

for a positive constant $C$ independent of $\epsilon$.

Proof of Proposition 3.1.

By [28] (proposition 13), we have that $\tilde{\mathbf{Q}}_{\epsilon}\in C^{\infty}\left(\tilde{\Omega};S_{0}\right)$. Similarly, we have that $\tilde{\mathbf{Q}}_{\infty}\in C^{\infty}\left(\tilde{\Omega};S_{0}\right)$. Since $\left\|\left(\tilde{g}_{d}\right)_{i}\right\|_{L^{\infty}\left(\partial\tilde{\Omega}\right)}\leq 1$ ($i=1,2$), we have $\left\|\left(\tilde{Q}_{\epsilon}\right)_{1i}\right\|_{L^{\infty}\left(\tilde{\Omega}\right)}\leq 1$ ($i=1,2$) [2] (proposition 2). Therefore, comparing Equation 2.35 and Equation 3.3, we have,

[TABLE]

Let $v\in C^{\infty}\left(\tilde{\Omega};\mathbb{R}\right)$ be a solution of,

[TABLE]

The solution $v$ only depends the domain $\tilde{\Omega}$ which is fixed. We can check that $-\epsilon^{-2}v$ and $\epsilon^{-2}v$ are the sub-solution and the super-solution for both matrix components of $\left(\tilde{\mathbf{Q}}_{\epsilon}-\tilde{\mathbf{Q}}_{\infty}\right)$ and the result then follows, i.e.,

[TABLE]

∎

We can compute the solution of the limiting problem Equation 3.3 exactly for a given Dirichlet boundary condition as shown below.

Proposition 3.2**.**

Let $\tilde{\mathbf{Q}}_{\infty}$ be the unique solution of Equation 3.3, on the rectangular domain $\tilde{\Omega}$ defined above, subject to the boundary condition $\left(\left(\tilde{Q}_{\infty}\right)_{11},\left(\tilde{Q}_{\infty}\right)_{12}\right)=\tilde{\mathbf{g}}_{d}$, on $\partial\tilde{\Omega}$ for given $0<d<\min\left\{a/4,b/4\right\}$. Then we have $\left(\tilde{Q}_{\infty}\right)_{12}\equiv 0$ in $\tilde{\Omega}$ and, for all $\left(x,y\right)\in\tilde{\Omega}$,

[TABLE]

Proof of Proposition 3.2.

The existence of a unique solution for the Laplace equation with Dirichlet conditions, as in Equation 3.3 is standard [9]. Let $\tilde{\mathbf{Q}}_{\infty}$ be the corresponding unique solution of Equation 3.3 with boundary condition $\left(\left(\tilde{Q}_{\infty}\right)_{11},\left(\tilde{Q}_{\infty}\right)_{12}\right)=\tilde{\mathbf{g}}_{d}$ on $\partial\tilde{\Omega}$, where $\left(\tilde{Q}_{\infty}\right)_{12}\equiv 0$ on $\partial\tilde{\Omega}$ from our choice of $\tilde{\mathbf{g}}_{d}$ above. Hence, we immediately have $\left(\tilde{Q}_{\infty}\right)_{12}\equiv 0$ in $\mathrm{Int}\left({\tilde{\Omega}}\right)$. To compute $\left(\tilde{Q}_{\infty}\right)_{11}$, we solve the following boundary-value problem for $x\in\left(0,a\right)$ and $y\in\left(0,b\right)$,

[TABLE]

Consider the scalar function $f\left(x,y;a,b\right)$ which satisfies the following boundary-value problem,

[TABLE]

We solve for $f\left(x,y;a,b\right)$ by separation of variables (also see [21]). A standard computation shows, using boundary conditions at $x=0$ and $x=a$,

[TABLE]

where $c_{k}^{\prime}$ and $c_{k}^{\prime\prime}$ are constants depend on $k$. We use the remaining boundary conditions at $y=0$ and $y=b$ and Fourier series methods to obtain

[TABLE]

Then, it is straightforward to verify that

[TABLE]

and the result Equation 3.8 now follows. ∎

On a square with $a=b$, the WORS solution is distinguished by $\tilde{\mathbf{Q}}=0$ along the square diagonals, in particular $\tilde{\mathbf{Q}}\left(\frac{a}{2},\frac{a}{2}\right)=0$ at the square centre. It is reasonable to ask if $\tilde{\mathbf{Q}}_{\infty}$ retains these nodal diagonal lines for $a\neq b$. A simple test is to check if $\tilde{\mathbf{Q}}_{\infty}\left(\frac{a}{2},\frac{b}{2}\right)\neq 0$ which would demonstrate the loss of the WORS structure for $a\neq b$. In Figure 3.1, we plot $\tilde{\mathbf{Q}}_{\infty}$ on a square ($a=b=1$) which is simply the WORS and $\tilde{\mathbf{Q}}_{\infty}$ on a rectangle ($a=1.5$, $b=1$) to illustrate the differences to the reader.

Proposition 3.3**.**

Let $\tilde{\mathbf{Q}}_{\infty}$ be the unique solution of Equation 3.3 on the rectangular domain $\tilde{\Omega}$, subject to the boundary condition $\left(\left(\tilde{Q}_{\infty}\right)_{11},\left(\tilde{Q}_{\infty}\right)_{12}\right)=\tilde{\mathbf{g}}_{d}$, for a given $0<d<\min\left\{a/4,b/4\right\}$. For $a>b=1$, we have

[TABLE]

so that we lose the cross structure of the WORS on all rectangles (similar arguments also apply to $0<a<b=1$).

Proof of Proposition 3.3.

Let $\tilde{\mathbf{Q}}_{\infty}^{\prime}$ denote the WORS solution of Equation 3.3 on a square domain $\Omega_{b,b}$ with boundary condition $\left(\left(\tilde{Q}_{\infty}^{\prime}\right)_{11},\left(\tilde{Q}_{\infty}^{\prime}\right)_{12}\right)=\tilde{\mathbf{g}}_{d;b,b}$ on $\partial\Omega_{b,b}$. To prove the result above, we show that $\left(\tilde{Q}_{\infty}\right)_{11}\left(a/2,b/2\right)>\left(\tilde{Q}_{\infty}^{\prime}\right)_{11}\left(b/2,b/2\right)=0$. We construct five auxiliary boundary value problems below and we illustrate the boundary conditions in Figure 3.2, using the boundary function in Equation 2.42. Let $I_{p}=\left[-p/2,p/2\right]$ for $p>0$.

[TABLE]

Using the maximum principle [9], we have $u_{1}\geq 0$ in $\Omega_{1}$ and $u_{2}\geq 0$ in $\Omega_{2}$, so that, $u_{1}>u_{2}$ on the vertical lines $x=\pm b/2$. Hence,

[TABLE]

By the strong maximum principle, $u_{1}\left(0,0\right)>u_{2}\left(0,0\right)$.

Similarly, we have that $u_{3}\leq 0$ on $\Omega_{3}$ and $u_{4}\leq 0$ in $\Omega_{4}$, so that

[TABLE]

By the strong maximum principle, $u_{3}\left(0,0\right)>u_{4}\left(0,0\right)$. We can also check that $u_{4}\left(x,y\right)=u_{5}\left(xb/a,yb/a\right)$, for $\left(x,y\right)\in\Omega_{4}$, since the Laplace operator is invariant with respect to uniform scaling and hence, $u_{4}\left(0,0\right)=u_{5}\left(0,0\right)$.

By superposition of boundary conditions, we see that the functions $\left(\tilde{Q}_{\infty}\right)_{11}$ and $u_{6}:=u_{1}+u_{3}$ only differ by a translation,

[TABLE]

and we have the similar relationship for $\left(\tilde{Q}_{\infty}\right)_{11}^{\prime}$ and $u_{7}:=u_{2}+u_{5}$,

[TABLE]

The result immediately follows since

[TABLE]

∎

Remark 3.1**.**

In the limit $d\to 0$, the function $\left(\tilde{Q}_{\infty}\right)_{11}$ converges pointwise to

[TABLE]

which is still smooth in $\mathrm{Int}\left({\Omega}\right)$, so that the maximum principle is valid. Thus, the result in Proposition 3.3 holds in the $d\to 0$ limit for $a>b$.

3.1.2 Weak Anchoring Condition

In Proposition 3.3, we show the perfect cross structure of the WORS on a square domain is lost as soon as we introduce geometrical anisotropy with $a\neq b$, at least with Dirichlet boundary conditions or strong anchoring. It is equally natural to ask if the WORS is an artefact of the strong anchoring condition i.e. do we lose the cross structure as soon as we relax the boundary conditions on the edges of the square or do we have a continuous transition from the WORS to a cross-free profile, as the anchoring is relaxed by means of a suitably defined anchoring parameter? The ‘cross’ structure refers to the fact that the WORS vanishes along the square diagonals and we refer to the nodal diagonal lines as a ‘cross’.

We work with $\mathbf{Q}$-tensors (defined in Equation 2.16) that are independent of the vertical dimension, on a square well $\Omega$ with $a=b=1$ and $L>0$. The key difference is that we impose a Durand-Nobili surface energy [29] on the lateral surfaces $\partial\Omega\times\left[0,h\right]$ (instead of a Dirichlet uniaxial tangent condition) with

[TABLE]

for an anchoring coefficient $W>0$, and $\mathbf{g}_{d}\left(x,y\right)=\tilde{\mathbf{g}}_{d}\left(x/L,y/L\right)$, $\left(x,y\right)\in\partial\Omega$. In the limit $W\to\infty$, we recover the strong anchoring condition in Equation 2.38. Using the standard scalings, the rescaled LdG free energy is

[TABLE]

Defining $\tilde{\mathbf{Q}}$ as in Equation 2.30 and $\epsilon$ as in Equation 2.31, we work with a reduced and rescaled LdG free energy of the form

[TABLE]

where the reduced anchoring coefficient is

[TABLE]

We treat $B$, $C$, $K$ to be fixed material constants in this manuscript and hence, we have $\tilde{W}\propto W\epsilon^{-1}$. Our aim is to investigate the effects of a mild relaxation of the strong anchoring conditions and hence, we assume that $W=\mathcal{O}\left({\epsilon}\right)$ as $\epsilon\to\infty$ so that $\tilde{W}=\alpha$ for some positive constant $\alpha$ independent of $\epsilon$, for large values of $\epsilon$. We can readily prove existence of minimizers of Equation 3.43 under these assumptions from the direct methods in the calculus of variations and for sufficiently large $\epsilon$, we can follow the arguments in [20] and [7] to demonstrate the uniqueness of critical points of Equation 3.43 in this limit. It is straightforward to show that minimizers of Equation 3.43 satisfy the following system of partial differential equations:

[TABLE]

where $\mathbf{\nu}_{\tilde{\Omega}}$ is the normal to the edges of the re-scaled rectangular domain, $\tilde{\Omega}$. Following the methods in [7], we can prove the existence of a Well Order Reconstruction Solution (WORS) for the system Equation 3.45 and Equation 3.46 with $\tilde{Q}_{12}=0$ everywhere and $\tilde{Q}_{11}=0$ on the square diagonals and the cross-structure is preserved.

Under these assumptions, we can also check that minimizers of Equation 3.43 converge uniformly to minimizers of the following limiting energy as $\epsilon\to\infty$:

[TABLE]

We can compute the limiting profile, $\tilde{\mathbf{Q}}_{\infty}$, exactly in terms of its matrix components $\left(\tilde{Q}_{\infty}\right)_{11}$ and $\left(\tilde{Q}_{\infty}\right)_{12}$, as shown below. We define the following function $J$ for $q\in W^{1,2}\left(\tilde{\Omega}\right)\cap C^{\infty}\left(\mathrm{Int}\left({\tilde{\Omega}}\right)\right)$, $\tau\in\mathbb{R}$, $g\in C^{\infty}\left(\partial\tilde{\Omega}\right)$,

[TABLE]

If $q$ is a critical point (e.g., local minimizer) of the energy $J$, then we necessarily have,

[TABLE]

This is simply the Laplace equation on a rectangular domain with Robin boundary conditions, which can be solved by means of separation of variables. It follows that

[TABLE]

For $\tau>0$, the boundary-value problem Equation 3.51 has a unique solution which can be computed exactly.

Remark 3.2** (Cross structure of $\tilde{\mathbf{Q}}_{\infty}$ with $\tilde{W}>0$ for a square with $a=b=1$).**

We have a unique solution $\left(\tilde{Q}_{\infty}\right)_{11}$ for Equation 3.52 and $\left(\tilde{Q}_{\infty}\right)_{12}$ for Equation 3.53 for $\tau>0$. Since $\left(\tilde{g}_{d}\right)_{2}\equiv 0$, this yields $\left(\tilde{Q}_{\infty}\right)_{12}\equiv 0$ on $\tilde{\Omega}$. Using the antisymmetry of $\left(\tilde{g}_{d}\right)_{1}$, we have

[TABLE]

which implies that

[TABLE]

Hence, we obtain

[TABLE]

from these relations quite easily. This simple argument demonstrates that $\tilde{\mathbf{Q}}_{\infty}$ vanishes along the two square diagonals $y=x$ and $y=1-x$ and the cross structure is retained in the limit too, under the assumption that $\tilde{W}$ remains bounded as $\epsilon\to\infty$.

Remark 3.3** (Explicit limiting solution with weak anchoring).**

Next, we analytically solve the limiting problem Equation 3.51 for $\left(\tilde{Q}_{\infty}\right)_{11}$. In the $d\to 0$ limit, the boundary function $\left(g_{d}\right)_{1}$ is either $+1$ or $-1$ on all four edges. Using ideas similar to Proposition 3.2, we can write $\left(\tilde{Q}_{\infty}\right)_{11}$ as

[TABLE]

where the function $f\in W^{1,2}\left(\tilde{\Omega}\right)\cap C^{\infty}\left(\mathrm{Int}\left({\tilde{\Omega}}\right)\right)$ is a solution of the following boundary-value problem,

[TABLE]

We can solve for $f$ by separation of variables. Skipping all the technical details, we have

[TABLE]

where $\left\{p_{k}\right\}$ is the set of solutions of the transcendental equation,

[TABLE]

(Also see [36] for similar calculations.)

To summarise, we investigate the effects of both geometrical anisotropy ($a\neq b$) and relaxed but strong anchoring on the WORS in this section. We lose the cross structure of the WORS for $a\neq b$. However, for $W=\mathcal{O}\left({\epsilon}\right)$ and for $\epsilon$ sufficiently large, we retain the cross structure of the WORS so that the WORS is not an artefact of the Dirichlet conditions. These properties are preserved in the limit since we can compute the limiting profiles, $\tilde{\mathbf{Q}}_{\infty}$ in both cases. Of course, as $W$ decreases and $\tilde{W}\to 0$ for $W\ll\epsilon$, then $\tilde{\mathbf{Q}}_{\infty}$ can be any constant matrix as expected, since there are no boundary constraints in this limit. In Figure 3.3, we plot $\tilde{\mathbf{Q}}_{\infty}$ on a square with $\tilde{W}=3$ and $\tilde{W}=10$ respectively. The cross structure is necessarily more diffuse with smaller values of $\tilde{W}$. In Figure 3.4, we illustrate similar numerical results for a rectangle with $a=1.5$, $b=1$; the cross structure disappears as we see nodal lines of $\tilde{\mathbf{Q}}_{\infty}$ along the shorter rectangular edges. We refer to such solutions as BD solutions in the rest of the text, by analogy with similar terminology in [34].

3.2 The $\epsilon\to 0$ Limiting Problem

Let $\tilde{\mathbf{Q}}_{\epsilon}$ be a global minimizer of the rescaled LdG free energy (Equation 2.32) subject to the fixed strong anchoring boundary conditions $\left(\left(\tilde{Q}_{\epsilon}\right)_{11},\left(\tilde{Q}_{\epsilon}\right)_{12}\right)=\tilde{\mathbf{g}}_{d}$ on $\partial\tilde{\Omega}$ for given $d>0$. For $\epsilon$ sufficiently small, in the interior, we have to leading order (more precisely, under suitable conditions, the error is $\mathcal{O}\left({\epsilon^{2}}\right)$ in $L^{\infty}$ norm [2]),

[TABLE]

where $\hat{\mathbf{n}}=\left(\cos\theta,\sin\theta\right)^{T}$ and $\theta$ is the solution of the Laplace equation

[TABLE]

subject to appropriately defined Dirichlet conditions [21], [28].

Hence, we can obtain good approximations to locally stable minimizers of Equation 2.32 for sufficiently small $\epsilon$, by studying Dirichlet boundary-value problems for the angle $\theta$ on a rectangle as given below (we set $d=0$ in the Dirichlet condition $\tilde{\mathbf{g}}_{d}$ for the sake of this computation),

[TABLE]

where $d_{1},d_{2},d_{3},d_{4}$ are arbitrary constants accounting for the constant boundary conditions on the edges of the rectangle. This is analogous to the approach in [22]. In the case of a rectangle with $a>b$, it is known that there are three competing equilibria—the diagonal states for which $\hat{\mathbf{n}}$ is aligned along a diagonal of the rectangle, the rotated R1 and R2 states for which $\hat{\mathbf{n}}$ rotates by $\pi$ radians between a pair of parallel horizontal edges and the rotated R3, R4 states for which $\hat{\mathbf{n}}$ rotates by $\pi$ radians between a pair of parallel vertical edges. For $a>b$, the R3, R4 states have lower energies than the R1, R2 states (see [32], [22] for details). We can obtain good approximations to the diagonal and rotated solutions by studying solutions of Equation 3.67 subject to the boundary conditions enumerated in Table 3.1, which are consistent with [23]. These boundary conditions enforce $\pm 1/4$ topological defects in the vector field, $\hat{\mathbf{n}}$, at the corners of $\tilde{\Omega}$. To be more precise, we say that a vertex defect has charge $t/4$ ($t$ is an odd integer) when the director rotates by $2\pi t/4$ radians between a pair of adjacent edges; the sign being determined by the sense of the rotation. With this convention, the $\pm 1/4$ defects at the vertex corners are the weakest admissible defects and are referred to as trivial cases in section 5.

For arbitrary constants $d_{1},\cdots,d_{4}$, the solution of Equation 3.67 and Equation 3.68 can be written in terms of a function $f$, where $f$ is a solution of the following boundary-value problem,

[TABLE]

The solution can be computed from Equation 3.12 in the limit $d/a\to 0$,

[TABLE]

with the symmetry $f\left(x,y;a,b\right)=f\left(a-x,y;a,b\right)$. Then

[TABLE]

The D1 state corresponds to

[TABLE]

with the symmetries

[TABLE]

Referring to Table 3.1, the R3 solution corresponds to

[TABLE]

with the symmetries

[TABLE]

In both cases, it suffices to consider the solution of Equation 3.67 on the quadrant $\left[0,a/2\right]\times\left[0,b/2\right]$. The D1 solution is subject to the boundary conditions,

[TABLE]

whereas the R3 solution is subject to,

[TABLE]

Thus, the R3 solution is a minimizer of the Dirichlet energy

[TABLE]

in a smaller space than the D1 solution, which is a minimizer of the Dirichlet energy in a larger space, and hence the D1 solutions cannot have higher energies than the competing R3 solutions. This gives a simple explanation of the fact that diagonal states are observed more frequently in experiments (see [22], [32]) since they have lower energies than their rotated counterparts, without the need for any elaborate calculations.

4 Bifurcation Diagrams as a Function of $\epsilon$ and Rectangular Anisotropy

In [30], the authors extensively discuss the reduced LdG equilibria on a square domain for low temperatures, as a function of the square length $L$. They numerically demonstrate that the WORS (featured by a vanishing $\mathbf{Q}$-tensor along the square diagonals; see Equation 2.19) is the unique solution for $L$ small enough, its bifurcation into stable diagonal branches as $L$ increases, a further bifurcation into two unstable BD branches (defined as solutions with $Q_{12}=0$ but without the property of $Q_{11}=Q_{12}=0$ at the square centre) and then a tracking of how the unstable BD branches bifurcate into unstable rotated solutions and the rotated states gain stability as $L$ increases. For $L$ large enough, there are six distinct stable solution branches—the D1 and D2 solutions corresponding to the two different diagonal solutions and four energetically degenerate R1, R2, R3, R4 solutions on a square. The rotated solutions are related to each other by a $\pi/2$ rotation on a square and are hence, energetically degenerate. As expected, we lose the degeneracy as soon as we break the symmetry of the square.

We follow the same paradigm on rectangular domains $\tilde{\Omega}$, to track solution branches as a function of the rescaled parameter $\epsilon$ (which is inversely proportional to the rectangular size) and the geometrical anisotropy, measured in terms of $a$ for fixed $b=1$. We use arc-continuation methods [16] to track different solution branches and it is possible that this method does not locate all solution branches, or misses some high energy unstable branches. However, this is a standard and well-accepted method for computing bifurcation diagrams for solutions of systems of partial differential equations.

We numerically find eight different kinds of solutions on a rectangle —the D1, D2, R1, R2, R3, R4 solutions as enumerated in the Equation 2.38 along with two solutions (only relevant for large $\epsilon$), labelled as BD1 and BD2. These solutions are described by their corresponding $\tilde{\mathbf{Q}}$-tensors, e.g., $\tilde{\mathbf{Q}}_{\mathrm{D1}}$, $\tilde{\mathbf{Q}}_{\mathrm{D2}}$, $\ldots$ for a given value of $\epsilon$. The BD1 and BD2 solutions are special since $\left(\tilde{Q}_{\mathrm{BD1}}\right)_{12}=\left(\tilde{Q}_{\mathrm{BD2}}\right)_{12}\equiv 0$ which implies that the corresponding $\hat{\mathbf{n}}$ is either $\hat{\mathbf{n}}=\left(1,0\right)$ or $\hat{\mathbf{n}}=\left(0,1\right)$ (see Equation 2.19) everywhere in the rectangular interior. BD1 refer to BD solutions that have $\hat{\mathbf{n}}=\left(0,1\right)$ at the centre of the rectangle with $\hat{\mathbf{n}}=\left(1,0\right)$ near the horizontal edges, $y=0$ and $y=b$. BD2 refer to BD solutions that have $\hat{\mathbf{n}}=\left(1,0\right)$ at the centre of the domain with $\hat{\mathbf{n}}=\left(0,1\right)$ near the vertical edges, $x=0$ and $x=a$. For $a=b$, the solutions BD1 and BD2 are energetically degenerate whereas for $a>b$, one can heuristically see that BD1 is unfavourable compared to BD2 since the transition layers for BD1 are located along the longer edges, $y=0$ and $y=b$. We point out that BD2 is the limiting profile, $\tilde{\mathbf{Q}}_{\infty}$ on a rectangle $\tilde{\Omega}$ with $a>b$ discussed in subsubsection 3.1.1. On similar grounds, we lose the degeneracy between the rotated states and the R1 and R2 states have higher energies than the R3, R4 states for $a>b$. We also study the stability of the solutions, and use prefixes ‘s’ and ‘u’ for stable and unstable solutions respectively. We test the stability ($L^{2}$-stability) of a solution by computing the smallest eigenvalue of the second variation of the discrete LdG free energy functional [35], and if the smallest eigenvalue is positive, the solution is locally stable.

We use the notation $E_{\textrm{c}}\left(\epsilon,a\right)$ to denote the Landau-de Gennes energy of the solution in class ‘c’ at given $\epsilon$ and geometrical aspect ratio $\delta=\frac{a}{b}=a$ (since $b=1$ in our simulations). We plot the energies of the solutions of the Euler-Lagrange equation with respect to the parameter $\epsilon$, see 4(a), 4(c) and 4(e). For small values of $\epsilon$ and $a>b$, there are three competing stable solutions, sD1, sR2 and sR3, and the energies are ordered as $E_{\textrm{sD1}}\left(\epsilon,a\right)<E_{\textrm{sR3}}\left(\epsilon,a\right)\leq E_{\textrm{sR2}}\left(\epsilon,a\right)$, which is consistent with the results in [22].

We track three different solution pathways, with three distinct initial conditions sD1, sR3 and sR2 respectively, for small values of $\epsilon$. We discuss each pathway separately. As $\epsilon$ increases, the solution sD1 transitions to sBD2, which is the unique limiting solution described in Proposition 3.1; we do not observe any unstable solutions in this pathway and the qualitative features seem to be independent of the geometrical anisotropy $a$. We define the parameter, $\epsilon_{\textrm{sD1}\to\textrm{sBD2}}\left(a\right)$, be the value of $\epsilon$ where this transition occurs, defined by $\left(\tilde{Q}_{\textrm{sBD2}}\right)_{12}\equiv 0$. As $a$ increases from $1$, $\epsilon_{\textrm{sD1}\to\textrm{sBD2}}\left(a\right)$ decreases. We illustrate this transition in Figure 4.1.

The second pathway has the competing stable rotated solution, sR3, as initial condition, and we gradually increase $\epsilon$ and track the corresponding changes. The second pathway is different from the first pathway described above. We observe the following sequence of transitions from sR3 to uR3, to uBD2, and finally to sBD2 and the qualitative features are independent of $a$. We can similarly define the parameters $\epsilon_{\textrm{sR3}\to\textrm{uR3}}\left(a\right)<\epsilon_{\textrm{uR3}\to\textrm{uBD2}}\left(a\right)<\epsilon_{\textrm{uBD2}\to\textrm{sBD2}}\left(a\right)$. As $a$ increases from $1$, $\epsilon_{\textrm{sR3}\to\textrm{uR3}}\left(a\right)$ and $\epsilon_{\textrm{uR3}\to\textrm{uBD2}}\left(a\right)$ increase, and $\epsilon_{\textrm{uBD2}\to\textrm{sBD2}}\left(a\right)\equiv\epsilon_{\textrm{sD1}\to\textrm{sBD2}}\left(a\right)$ decreases. As $a\to\infty$, the differences among these three critical parameters tend to zero. We illustrate this sequence of structural transitions in Figure 4.2.

The third pathway is quite different, for which we use sR2 as initial condition for $\epsilon$ small enough and then gradually increase the parameter $\epsilon$. In fact, the qualitative features of the pathway also depend on the value of $a>1$. At small $a$, e.g. $a=1.25$, the solution sR2 changes to uR2, then to uBD1, and finally the solution branch disappears at an end point. We can similarly define the parameters $\epsilon_{\textrm{sR2}\to\textrm{uR2}}\left(a\right)<\epsilon_{\textrm{uR2}\to\textrm{uBD1}}\left(a\right)<\epsilon_{\textrm{end}}\left(a\right)$, such that we cannot find further solutions on this branch by arc-continuation methods for $\epsilon>\epsilon_{\textrm{end}}\left(a\right)$. At intermediate $a$, e.g., $a=1.5$, the solution sR2 transitions to uR2 and the solution branch disappears on increasing $\epsilon$, and we have $\epsilon_{\textrm{sR2}\to\textrm{uR2}}\left(a\right)<\epsilon_{\textrm{end}}\left(a\right)$. There is no intermediate transition to uBD1 for intermediate values of $a>1$. At larger $a$, we observe the termination of the branch at sR2. As $a$ increases from $1$, $\epsilon_{\textrm{end}}\left(a\right)<\epsilon_{\textrm{sD1}\to\textrm{sBD2}}\left(a\right)$ decreases. We illustrate these structural transition pathways in Figure 4.3, for different values of $a$.

We plot bifurcation diagrams for the reduced LdG equilibria for three values of $a=1,1.25,1.5$ respectively, following the same measures as in [30]. The figures are labelled as 4(b), 4(d) and 4(f). These plots illustrate the distinct solution pathways and changes in the number of solutions as $\epsilon$ decreases, and the sensitivity of the solution landscape to the aspect ratio $a$, although more exhaustive studies are needed to this effect. We use the quantities $\tilde{\Omega}^{-1}\int_{\tilde{\Omega}}\tilde{Q}_{11}^{2}$, $\tilde{\Omega}^{-1}\int_{\tilde{\Omega}}\tilde{Q}_{12}^{2}$ and $\tilde{\Omega}^{-1}\int_{\tilde{\Omega}}s^{2}$ in Figure 4.5, from which we deduce that the structural transition from sD1 to sBD2 is a first order structural transition i.e. discontinuous as a function of $\epsilon$.

The pertinent question is—how are the solution landscapes on a rectangle different from the solution landscapes on a square? There are no striking differences except perhaps that the sR2 solution pathway is not connected to the sD1 and sR3 pathways directly for $a>1$. Similarly, we do not find a sBD1 state for a rectangle even for smaller values of $a$, for which the energetic penalty of the transition layers along the longer horizontal edges is lower. It is interesting that sD1 and uBD2 transition to sBD2 at the same value of $\epsilon$ on two distinct solution pathways. It is worth noting that BD states are always unstable on a square and this is an illuminating example of how geometrical anisotropy can stabilise BD states on a rectangle. Further, the solution pathways do not appear to be highly sensitive to $a$, although this could be a limitation of the numerical methods. Further, as $a$ increases from $1$, and for small $\epsilon$, $E_{\textrm{sD1}}\left(\epsilon,a\right)$ and $E_{\textrm{sR2}}\left(\epsilon,a\right)$ increase, while $E_{\textrm{sR3}}\left(\epsilon,a\right)$ decreases. When $a=1$, $E_{\textrm{sR3}}\left(\epsilon,a\right)=E_{\textrm{sR2}}\left(\epsilon,a\right)$. As $a\to\infty$, $E_{\textrm{sR3}}\left(\epsilon,a\right)\to E_{\textrm{sD1}}\left(\epsilon,a\right)$, for any fixed and small $\epsilon$, which agrees with the results in [22]. An open question pertains to how a system would transition from the sR2 solution to the sBD2 solution as $\epsilon$ increases, since sBD2 is the unique solution for $\epsilon$ large enough?

5 Relaxation Mechanisms for Non-trivial Topologies

We can classify nematic equilibria on a rectangular domain with tangent boundary conditions in terms of their structural details near the rectangular vertices. More precisely, a nematic state, in our framework, is described by a 2D LdG $\mathbf{Q}$-tensor. In terms of the corresponding director, $\bf n=\left(\cos\theta,\sin\theta\right)$, we necessarily need that $\theta$ is some multiple of $\pi$ on the horizontal edges and that $\theta$ is an odd multiple of $\pi/2$ on the vertical edges. In the simplest case, the director rotates by $\pm\pi/2$ radians between a pair of adjacent edges, as is the case with the diagonal and rotated solutions described above. We refer to these solutions as having ‘trivial’ topologies. However, it is admissible for the director to rotate by $n\pi/2$ radians between a pair of adjacent edges, for an odd integer $\left|n\right|\geq 3$, and such equilibria are deemed to have ‘non-trivial’ topologies. The degree of a vertex is defined to be $\omega=\frac{1}{2\pi}\frac{n\pi}{2}=\frac{n}{4}$ for some odd integer $n$, which is positive if the rotation is anticlockwise along an anticlockwise-oriented arc connecting the two intersecting edges at the vertex in question, and negative otherwise.

It is natural to ask if non-trivial topologies can be realised in practice? One can conceive a three-dimensional set-up of a shallow well $\mathcal{B}$ with a rectangular cross-section $\Omega$, with Dirichlet uniaxial tangent conditions on the lateral surfaces and suitable surface anchoring energies on the top and bottom surfaces. As argued in previous sections, under physically reasonable assumptions, we can reduce to a rescaled two-dimensional problem for which we study planar profiles on a rectangle $\tilde{\Omega}=\left[0,a\right]\times\left[0,b\right]$ with tangent conditions for the directors on the edges. The liquid crystal molecules could be locally rotated near chosen vertices, inducing non-trivial behaviour or non-trivial topologies. We expect these non-trivial topologies to relax to trivial topologies once the local rotation is removed and the relaxation pathways will be sensitive to $\epsilon$ and $\delta=\frac{a}{b}$, with interesting consequences for the corresponding optical properties and the relaxation pathways could even offer new switching mechanisms between diagonal and rotated states. Whilst this is largely a proof of concept at this stage, we explore some relaxation mechanisms from non-trivial to trivial topologies in this section as described below.

The first step is to construct an initial condition with a non-trivial topology for our numerical solvers. As a concrete example, we numerically solve

[TABLE]

Let $\theta_{0}$ denote the solution of the boundary-value problem above; then $\theta_{0}$ has a $5/4$-degree defect at $\left(0,0\right)$, a $-3/4$-degree defect at $\left(a,b\right)$, and two $-1/4$-degree defects at $\left(0,b\right)$ and $\left(a,0\right)$. We refer to $\left(0,0\right)$ and $\left(a,b\right)$ as non-trivial vertices and to $\left(0,b\right)$ and $\left(a,0\right)$ as being trivial vertices. We use $\theta_{0}$ to define a $\tilde{\mathbf{Q}}_{0}$-tensor field on the rescaled rectangle $\tilde{\Omega}$, by setting $\left(\tilde{Q}_{0}\right)_{11}=\cos\left(2\theta_{0}\right)$, $\left(\tilde{Q}_{0}\right)_{12}=\sin\left(2\theta_{0}\right)$. The boundary condition $\left(\left(\tilde{Q}_{0}\right)_{11},\left(\tilde{Q}_{0}\right)_{12}\right)=\tilde{\mathbf{g}}_{d}$ on $\partial\tilde{\Omega}$ for a given $d>0$ (sufficiently small). Then we can compute the 2D LdG free energy minimizer with $\tilde{\mathbf{Q}}_{0}$ as initial condition, using the steepest descent method [5]. It is clear that the numerical procedure will converge to either sBD2 (when $a>b$) for large $\epsilon$ or one of the diagonal or rotated states for small $\epsilon$ and we record the intermediate states during the energy minimization procedure to study the relaxation mechanisms as a function of $\epsilon$ and the aspect ratio $\delta$.

In the numerical implementation, we fix $d=0.03$ for the boundary condition, and $b=1$ . For $a=1$, i.e. a square domain, we find three different relaxation modes for different ranges of $\epsilon$ in the 2D LdG model (see Equation 2.30 to Equation 2.42). For small $\epsilon$ (for example, $\epsilon=0.05$), the non-trivial defects with $\left|\omega\right|\geq 3/4$, split into multiple defects and pairs of these newly created defects merge near vertices to give us the diagonal state that connects the two non-trivial vertices (D1). More precisely, we see two $1/2$-defects emanate from the $5/4$-vertex and two $-1/2$-defects emanate from the $-3/4$-vertex, and these two pairs attract each other and annihilate near the two trivial vertices. These leaves us with two $1/4$-degree vertices at $\left(0,0\right)$ and $\left(a,b\right)$ and the D1 state connecting them. For intermediate values of $\epsilon$ (for example, $\epsilon=0.09$), the $-3/4$-defect expels a $-1/2$-defect into the interior and the $5/4$-defect expels three $+1/2$-defects into the interior. A $\pm 1/2$-pair annihilates in the bulk and the remaining two $+1/2$-defects move along the shorter edges to the trivial vertices. In other words, we observe defect annihilation in the bulk and get the diagonal solution with the opposite diagonal orientation (D2); for large $\epsilon$ (for example, $\epsilon=0.4$), the initial state relaxes to the sWORS state as expected, through almost isotropic states in the bulk. See Figure 5.1 for more details. For $a=5$, i.e., a rectangular domain, we find three analogous relaxation modes for different values of $\epsilon$. For small $\epsilon$ (for example, $\epsilon=0.05$), the non-trivial defects (with degree $\left|\omega\right|\geq 3/4$) split into multiple defects and some of these defects move towards each other in the bulk and others migrate to the nearest vertex to yield the rotated state (R3); for intermediate $\epsilon$ (for example, $\epsilon=0.09$), we observe the defect splitting near the non-trivial vertices and the expelled defects move along the long edges to either annihilate each other or change the topology of a vertex, yielding the diagonal solution (D2); for large values of $\epsilon$ (for example, $\epsilon=0.4$), the initial state relaxes to the sBD2 state through almost isotropic states (with $\tilde{\mathbf{Q}}=0$) in the interior. See Figure 5.2.

These examples are certainly not exhaustive. However, they do highlight certain generic concepts—the non-trivial vertices will relax into trivial vertices by expelling the excess charge, $\eta=\left(|n|-1\right)/2$, as a combination of $\pm 1/2$ defects in the director field. This combination is not unique, and we cannot predict the combination as a function of $\epsilon$ and $\delta$ at present. For small values of $\epsilon$ and for $\delta=1,5$, the transient states are ordered (with relatively high scalar order parameter) and the different relaxation mechanisms will offer different optical properties. Intermediate values of $\epsilon$ exhibit transient states with almost line defects i.e. lines of low order along the rectangular edges or in the interior, surrounded by regions of high order. As $\epsilon$ increases further, the transient states become increasingly disordered and again, these could be interesting examples of ‘wetted’ nematic states. Further, it seems that non-trivial topologies can relax to either the diagonal or rotated states and there appears to be no clear selection mechanism. Therefore, these numerical experiments offer interesting new routes for exploiting non-trivial topologies— to create new defects, control defect dynamics, create new wetted states and even offer new pathways between diagonal and rotated states for switching dynamics. The role of $\epsilon$ is relatively clear, both for the size of the disordered regions in the transient states (the size increases linearly with $\epsilon$ as expected and hence, defects can move more easily for smaller values of $\epsilon$) and for the choices of the final relaxed state but the role of the geometrical anisotropy is less clear from these simulations.

6 Conclusion

We have systematically studied reduced 2D nematic equilibria on rectangles, as a function of a material-dependent and geometry-dependent parameter, $\epsilon$. Our results carry over immediately to the three-dimensional framework at the fixed temperature $A=-\frac{B^{2}}{3C}$. Following the methods in [34], the qualitative conclusions will hold for all $A<0$, at least in the limiting cases.

We compare our results on a rectangle with those on a square, with analytic results in two asymptotic limits: the $\epsilon\to\infty$ and the $\epsilon\to 0$ limits respectively. Our first result concerns the loss of the WORS (featured by two isotropic mutually perpendicular defect lines connecting two pairs of diagonally opposite vertices) on a rectangle with Dirichlet or infinitely strong anchoring on the edges. On a rectangle, we get two isotropic defect lines connecting two distinct pairs of adjacent vertices, for nano-scale rectangles in the $\epsilon\to\infty$ limit. Interestingly, the WORS structure survives on squares with relaxed but sufficiently strong anchoring in the $\epsilon\to\infty$ limit, as can be deduced from a uniqueness and symmetry argument. Our arguments are quite generic and we can extend them to arbitrary 2D polygons. In [14], the authors study reduced 2D LdG equilibria on regular polygons and in the $\epsilon\to\infty$ limit, they report that the limiting profiles have a unique isotropic point at the centre of the regular polygon, with the exception of the square. We can work with irregular polygons too and we speculate that the location of the zeroes or isotropic points are determined by the geometrical anisotropies. In fact, we conjecture that we may only have isotropic interior point defects on generic 2D polygons with tangent boundary conditions in the $\epsilon\to\infty$ limit and the special symmetries of the square and rectangle yield higher-dimensional line defects in this limit. It remains to be seen if we can obtain multiple interior point defects in this distinguished limit by exploiting the geometrical features in 2D.

In the $\epsilon\to 0$ limit for micron-scale or larger geometries, we study the interplay between three competing stable equilibria: the diagonal states and two different types of rotated states: the R1, R2 states and the R3, R4 states. Since the rectangle has lesser symmetry than a square and for $a>b$, the R3 and R4 states have lower energies than the R1, R2 states. We use a simple symmetry based argument to deduce that the diagonal states have lower energy than the competing R3, R4 states, in the limiting scenario. The asymptotic results are complemented by numerical studies of the solution bifurcations on rectangles, as a function of $\epsilon$ and the geometrical anisotropy. We find eight different solutions—the two diagonal solutions, the R1, R2, R3, R4 solutions as enumerated above and the BD1 and BD2 solutions respectively. The BD2 solutions are precisely the limiting profiles in the $\epsilon\to\infty$ limit. Interestingly, the BD solutions with transition layers/nodal (zero) lines near a pair of opposite edges are always unstable on a square. In the case of a rectangle with $a>b$, the BD2 solutions are stable whilst the BD1 solutions are unstable. The geometrical anisotropy essentially disconnects the R1, R2 solutions from the remainder of the bifurcation diagram. The effects of the geometrical anisotropy are more pronounced on the R1, R2 solution pathway in the sense that the structural transitions along this pathway depend on the geometrical parameter $a$. In fact, the BD1 solution (although unstable) is only observed for small values of $a$. This is an interesting example of the effects of $a$, or geometrical anisotropy in our framework.

In the last part, we study the relaxation mechanisms for tangent states with non-trivial topologies into the stable diagonal or rotated states, or the sBD2 state, for different values of $\epsilon$ and the rectangular aspect ratio. In [27], the authors perform related numerical experiments in an Oseen-Frank framework which can be identified with our 2D LdG framework under the restriction of constant $s$; the authors predict that the relaxation mechanisms should either concentrate near the edges or take place in the interior according to the aspect ratio. Our numerical experiments illustrate the effects of $\epsilon$; we find that the relaxation mechanisms are dominated by ejections of bubbles of ‘non-trivial topology’ and the size of these bubbles depends on $\epsilon$. We do not observe clear trends with respect to the aspect ratio possibly because our numerical examples are not exhaustive, the relaxation mechanisms are not unique and we have not been able to locate the aspect-ratio dependent relaxation mechanisms or because the Oseen-Frank studies in [27] are restrictive and lose relevance in the more comprehensive Landau-de Gennes framework. These relaxation mechanisms are of practical relevance for transition pathways between equilibria, switching mechanisms etc. and we will study these mechanisms more carefully in the future. To conclude, our model examples on a rectangle shed useful insight into the consequences of geometrical anisotropy, how anisotropy can distort and displace defects in stable equilibria, stabilise states, disconnect solution branches and in doing so, introduce new possibilities for multi-stability in tailor-made nematic-filled geometries.

7 Acknowledgement

The authors acknowledge support from a Royal Society Newton Advanced Fellowship. LF acknowledges support from the Visiting Postgraduate Scholar programme at University of Bath. LF thanks Yucen Han for comments. LF and LZ were partially supported by National Natural Science Foundations of China (NSFC 11871339, 11861131, 11571314). AM acknowledges support from Shanghai Jiao Tong University as a Visiting Professor. AM was supported by an EPSRC Career Acceleration Fellowship EP/J001686/1 and EP/J001686/2 and is supported by an OCIAM Visiting Fellowship and the Keble Advanced Studies Centre. The authors thank Giacomo Canevari for very useful discussions on the weak anchoring section.

Reference