## Abstract

We numerically study structural transitions inside shallow sub-micrometre scale wells with square cross section, filled with nematic liquid crystal material. We model the wells within the Landau–de Gennes theory. We obtain two qualitatively different states: (i) the diagonal state with defects for relatively large wells with lateral dimension greater than a critical threshold and (ii) a new, two-dimensional star-like biaxial order reconstruction pattern called the well order-reconstruction structure (WORS), for wells smaller than the critical threshold. The WORS is defined by an uniaxial cross connecting the four vertices of the square cross section. We numerically compute the critical threshold in terms of the bare biaxial correlation length and study its dependence on the temperature and on the anchoring strength on the lateral well surfaces.

## 1. Introduction

Liquid crystals (LCs) are mesogenic phases of matter with physical properties intermediate between those of the conventional solid and liquid phases of matter [1]. Nematic LCs are the simplest type of LCs; nematics are complex anisotropic liquids with a certain degree of long-range orientational ordering. The structural properties of a nematic LC can be easily manipulated by incident light, electric or magnetic fields and mechanical effects [2]. As a consequence, they exhibit a wealth of diverse physical phenomena and are often used as a testing ground for fundamental physics in condensed matter systems [3]. Nematic LCs now have widespread applications across science and technology, e.g. electro-optic devices, sensors and notably, the huge and thriving LC display (LCD) industry around the globe [4–6].

Nematics in confinement have generated tremendous interest in recent years [2,7–9]. Such confined LC systems typically exhibit topological defects or localized material imperfections [10,6]. Defect core structures are poorly understood; experimental and numerical studies indicate that defect cores can be associated with uniaxial–biaxial structural transitions, i.e. structural transitions from a phase with a single distinguished direction (uniaxial) to a phase with two distinguished directions (biaxial) [2,9]. Such structural transitions were first reported in the benchmark paper [9]. In [9], the authors numerically analyse the interior of a defect core; the nematic LC is in a uniaxial phase at the defect core and in a uniaxial phase away from the defect core. The corresponding directors/distinguished directions (2.2) are perpendicular to one another and the LC material mediates between the two conflicting uniaxial states via a continuum of biaxial states. These transitions, often referred to as ‘order reconstruction’ (OR) phenomena were subsequently reported in a batch of numerical papers [11–15], which indicate that OR phenomena typically occur in severely confined systems with comparable material and geometrical length scales or/and under the action of intense electric fields. For example, in [16], the authors study the OR phenomenon for a twisted hybrid nematic cell. The authors numerically find two classes of solutions: (i) the conventional purely uniaxial twisted solution wherein the molecules twist uniformly between the two boundary plates and (ii) the OR untwisted solution characterized by a biaxial band around the centre of the cell. The OR solution is the only observable solution in nano-confined geometries where the cell gap is comparable to a material-dependent length scale, known as the biaxial correlation length [11]. Experimental studies of the OR phenomenon have been reported in, for example, [12,17] and such studies are particularly relevant for new LC systems with intrinsic biaxiality and their potential applications in nano-science and technology.

Our work is motivated in part by the results in [4,18] with a view to characterize new two-dimensional OR patterns in prototype geometries. This is an interesting problem since classical OR patterns are typically one-dimensional with the structural characteristics varying along one coordinate direction, e.g. along the radial direction in a spherical droplet [9], along the vertical direction in a cylinder [13,14] or along the normal direction to parallel bounding plates [16], etc., and two-dimensional OR patterns offer new scientific possibilities. In [4,18], the authors study a model LC device comprising a periodic array of shallow micrometre-sized square wells. These square wells are filled with nematic LC; the surfaces are treated so as to induce tangential or planar boundary conditions [4,5], i.e. on the well surfaces, the LC molecules are constrained to be in the plane of the surfaces. This relatively simple geometry is bistable or multistable in the sense that there are multiple stable LC textures with contrasting optical properties [18]. In [4,18], the authors work within the Landau–de Gennes (LdG) theory for nematic LCs [3]. Given that the wells are shallow, they argue that it suffices to model two-dimensional structural variations within the plane of the square and work with a strictly two-dimensional model, i.e. study two-dimensional configurations on the square. They work with micrometre-size wells so that the square size, denoted by *R*, is much greater than a material length scale known as the bare biaxial correlation length [13,14]. The bare biaxial correlation length is a measure of the correlation in nematic order, at the nematic-isotropic transition temperature, and is typically a few tens of nanometres [8,12,19]. In this limit, the stable LC states are effectively uniaxial almost everywhere, except for small localized defects around the square vertices [4]. The authors find two stable LC states: (i) the diagonal state with the average uniaxial alignment along one of the square diagonals and (ii) the rotated state where the average direction of alignment rotates by 180° between a pair of opposite square edges.

We carefully examine the effects of severe confinement on the same model system. In particular, we work with a three-dimensional model, under the assumption that all structural details are independent of the vertical coordinate. Our three-dimensional modelling incorporates biaxiality, which is outside the scope of the two-dimensional model employed in [4,18]. We measure the square cross-section size, *R*, in units of the bare biaxial correlation length, *ξ*^{(0)}_{b}, and obtain two qualitatively different LC textures as a function of *R*: (i) diagonal structure with defects (DSD) for large *R* and (ii) the WORS for small *R*. The DSD is qualitatively similar to the diagonal solution in [4,18], and our three-dimensional modelling illustrates the biaxial structure of the defect cores around the square vertices. The WORS is a new two-dimensional OR pattern connecting the four square vertices characterized by a star-shaped rim of maximal biaxiality connecting the four square vertices. The biaxial rim separates a negatively ordered uniaxial state along the square diagonal from the positively ordered uniaxial states near the square edges. In other words, the system mediates between different uniaxial states (the positively ordered boundary states and the negatively ordered diagonal state) via a two-dimensional star-shaped ring of intermediate biaxial states. The WORS is the only stable LC pattern for small wells, with *R*≤*R*_{c}, where *R*_{c} is a critical threshold length. In §4, we compute *R*_{c} in units of *ξ*^{(0)}_{b}, systematically investigate how *R*_{c} depends on temperature and anchoring strength on the lateral boundaries, to obtain quantitative information about the stability of the WORS.

This paper is organized as follows. In §2, we review the LdG theory for nematic LCs. In §3, we outline the model geometry, the relevant parametrizations, the material length scales and the governing Euler–Lagrange equations. In §4, we present our numerical results and summarize our conclusions and future outlook in §5.

## 2. Theoretical background

Let **Q**-tensor parameter [1,3]. The **Q**-tensor order parameter can be viewed as a macroscopic measure of the system anisotropy, i.e. a measure of the anisotropy in dielectric response to electric fields or anisotropy in magnetic susceptibility or response to magnetic fields. Mathematically, **Q** is a symmetric, traceless 3×3 matrix and from the spectral decomposition theorem the **Q**-tensor can be written in terms of its eigenvalues, {*s*_{i}}, and eigenvectors, {*e*_{i}}, as shown follows [2]:
**Q**=0, (ii) uniaxial phase if **Q** has a pair of degenerate non-zero eigenvalues, and (iii) the biaxial phase if **Q** has three distinct eigenvalues [3]. In the uniaxial case, **Q** can be concisely expressed as
*director field* *e*_{1}, and the scalar order parameter *S*. In other words, *e*_{1}, is the eigenvector with the non-degenerate eigenvalue and all directions perpendicular to *e*_{1} are physically equivalent. In the biaxial phase, **Q** can be written as
*e*_{i}} are the eigenvectors and (*S*,*U*) are two scalar order parameters that measure the degree of orientational ordering about the different eigenvectors [20]. A quantitative measure of the degree of biaxiality is defined to be
*β*^{2}∈[0,1] and *β*^{2}=0 for all uniaxial states [5]. The maximal value, *β*^{2}=1, is attained if and only if one of the eigenvalues, *s*_{i}, vanishes.

The LdG theory is a variational theory and the stable physically observable LC states correspond to minimizers of an appropriately defined LdG energy functional [5,20]. We work with the following LdG free energy
*d*^{3}** r** and

*d*

^{2}

**denote the volume and area measures, respectively [2,21]. The condensation (**

*r**f*

_{c}), elastic (

*f*

_{e}) and surface (

*f*

^{(i)}

_{s}) free energy densities are given by [2]

*A*

_{0},

*B*and

*C*are material constants,

*T*denotes the temperature and

*T*

_{*}is a characteristic material-dependent supercooling temperature. One can readily check that all stationary points of

*f*

_{c}are necessarily uniaxial and that there are three characteristic temperatures: (i)

*T*=

*T*

_{*}below which the isotropic phase is unstable, (ii)

*T*

_{IN}=

*T*

_{*}+

*B*

^{2}/27

*A*

_{0}

*C*associated with a first-order isotropic–nematic phase transition, and (iii) the nematic superheating temperature,

*T*

_{**}=

*T*

_{*}+

*B*

^{2}/24

*A*

_{0}

*C*, above which ordered nematic states do not exist. For temperatures

*T*<

*T*

_{IN}, the condensation enforces a uniaxial scalar order parameter,

*T*>

*T*

_{IN},

*f*

_{c}pushes the system towards the disordered isotropic phase. The elastic energy density penalizes any spatial inhomogeneities in the system and we work within the simplest one-constant approximation in (2.7). The elastic constant is denoted by

*L*>0 in (2.7). The surface anchoring energy density,

*f*

^{(i)}

_{s}, enforces a preferred orientation, encoded by

**Q**

^{(i)}

_{s}, on the

*i*th lateral surface. The anchoring coefficient,

*w*

^{(i)}, is a measure of the strength of the surface anchoring and the

*S*=

*S*

_{eq}(

*T*). However, the boundary conditions typically require that both conditions cannot be simultaneously satisfied in confined geometries, resulting in inhomogeneous pattern formation. In what follows, we compare the competing effects of the elastic contributions, condensation term and the surface effects by means of two different length scales: a material-dependent length scale known as the bare biaxial correlation length and the surface extrapolation length [12,13].

## 3. Modelling framework

### (a) Geometry

We study nematic samples inside a three-dimensional well with square cross section, *h* is the well height, *R* is the lateral dimension and *h*≪*R*, consistent with the assumptions in [4,18]. We take *h*=*R*/10.^{1} The coordinate unit-vectors are denoted by (*e*_{x},*e*_{y},*e*_{z}). We assume that the structural characteristics only depend on the (*x*,*y*)-spatial coordinates and are independent of the vertical *z*-coordinate; an assumption often used for shallow wells [4,18]. The boundary plates are located at {*x*=0,*R*}, {*y*=0,*R*} and {*z*=0,*h*}, and all plates are treated to induce tangent or planar degenerate boundary conditions. The tangent boundary conditions are implemented using a combination of strong and weak anchoring and natural boundary conditions.^{2} On the lateral surfaces, *x*={0,*R*} and *y*={0,*R*}, the tangent conditions are implemented via Dirichlet conditions/strong anchoring or weak anchoring/surface energies. In the strong anchoring case, the corresponding **Q**-tensor is prescribed to be strictly uniaxial, with director parallel to *e*_{x} and *e*_{y}, on the *xz*- and *yz*-surfaces, respectively. We impose natural boundary conditions on the top and bottom faces, *z*=0 and *z*=*h*. All numerical simulations are initialized using planar initial conditions. The planar initial conditions, combined with the natural boundary conditions, bias the solutions to respect the planar degenerate conditions on the top and bottom plates so that no particular direction in the (*x*,*y*)-plane is distinguished or singled out.

### (b) Parametrization

We adopt the following three-dimensional parametrization of the **Q**-tensor order parameter, as given in [13,22]:
*q*_{1}, *q*_{2} and *q*_{3} are independent of *z* and only depend on *x* and *y*. We assume *e*_{3}=*e*_{z} is always an eigenvector of **Q**. The remaining two eigenvectors are allowed to rotate with respect to the reference frame, (*e*_{x},*e*_{y}), by an angle *φ*∈[0,*π*). We have
*q*_{2} is a measure of the departure of the eigenframe, (*e*_{1},*e*_{2},*e*_{3}), from the Cartesian frame, (*e*_{x},*e*_{y},*e*_{z}), and the two frames coincide when *q*_{2}=0. The eigenvalues {*s*_{i}} are explicitly given by *s*_{3}=−2*q*_{3}. In particular, the biaxiality parameter is defined to be

We briefly compare the parametrization (3.2) with the parametrization in [4] where the authors study strictly two-dimensional LC configurations on a square domain; the **Q**-tensor is then given by
**Q**^{3}=0, the biaxiality parameter in (3.5) has no physical meaning and all defects correspond to isotropic regions or locally melted regions [4]. The two-dimensional parametrization (3.6) is a special case of (3.2) with *q*_{3}=0. The representation (3.2) contains more information than (3.6) as it allows us to investigate structural variations in both *q*_{3} and *β*^{2}, both of which are outside the scope of the two-dimensional representation (3.6). In particular, the parametrization (3.2) can resolve the rich biaxial structure of defect cores.

On the faces, *x*=0 and *x*=*R*, we set the preferred **Q**-tensor to be
*y*=0 (bottom) and *y*=*R* (top), we set
**Q**=**Q**^{(y)}_{s} on the respective lateral surfaces. For systems with finite anchoring, we use the prescribed **Q**-tensors in (3.7) and (3.8) as the preferred state, **Q**^{(i)}_{s}, in the surface anchoring energy density in (2.8). On the surfaces *z*=0 and *z*=*h*, we impose natural boundary conditions i.e. require that
*ν*=±*e*_{z} is the outward normal to these faces in the (*x*,*y*)-plane. This is consistent with the assumed invariance in the *z*-direction and when combined with planar initial conditions for numerical simulations, ensures that tangent boundary conditions are preferentially imposed on these surfaces.^{3}

### (c) Scaling

Let *t* be a dimensionless temperature defined to be *t*=(*T*−*T*_{*})/(*T*_{**}−*T*_{*}). Therefore, *t*(*T*=*T*_{*})=0, *t*(*T*=*T*_{IN})=8/9 and *t*(*T*=*T*_{**})=1. We define
*S*_{**}=*S*_{eq}(*T*_{**})=*B*/4*C*.

The structural characteristics of static LC equilibria are dictated by a complex interplay between material properties and surface anchoring effects. The key material length scale is the biaxial correlation length defined to be
*ξ*^{(0)}_{b} as a characteristic temperature-independent material length scale and for conventional nematics, *ξ*^{(0)}_{b} is of the order of a few tens of nanometre [19,23]. The strength of the surface interactions is traditionally described in terms of the surface extrapolation length [14,21].
*L* are around 10^{−12}*N* and *w*^{(i)}∈(10^{−8},10^{−3}) *N*/*m* [5].

We further define the following dimensionless quantities: *F* is given as in (2.5), and the corresponding dimensionless free energy densities are defined to be

In the expression for *φ*_{s} is either *φ*_{s}=0 (plates at *y*=0 and *y*=*R*) or *φ*_{s}=*π*/2 (plates at *x*=0 and *x*=*R*). The strong anchoring limit, *τ* is large compared with the scaled elastic constant, (*ξ*^{(0)}_{b}/*R*)^{2}, then the solution will predominantly minimize the condensation energy, i.e. be largely uniaxial with constant order parameter, *S*=*S*_{eq}(*T*) (at least away from defects). On the other hand, if *τ* and (*ξ*^{(0)}_{b}/*R*)^{2} are of comparable magnitude, then elastic distortions and deviations from the condensation energy minima (uniaxial phases with *S*=*S*_{eq}(*T*)), e.g. biaxiality, have comparable energetic costs and, hence, OR patterns are energetically viable. Similarly, if *d*, we compute phase diagrams for uniaxial–biaxial structural transitions, as a function of *R*, anchoring strength and temperature, and these phase diagrams corroborate our heuristic insights.

### (d) Euler–Lagrange equations

We use standard methods in calculus of variations to compute the Euler–Lagrange equations for extremal points of the LdG energy functional [24]. The Euler–Lagrange equations are a coupled system of elliptic partial differential equations for (*q*_{1},*q*_{2},*q*_{3}) as shown follows:
_{⊥}=∂^{2}/∂*x*^{2}+∂^{2}/∂*y*^{2}.

The boundary conditions on the plates *x*=0 and *x*=*R* are
*d*^{(x)}_{e} is a measure of the anchoring strength, *w*^{(x)}_{s} on each plate.

Similarly, we have the following boundary conditions on the plates *y*=0 and *y*=*R*:
*d*^{(y)}_{e} is a measure of the relative anchoring strength. The corresponding strong anchoring conditions are {*q*_{1}=−*τ*/2, *q*_{2}=0 and *q*_{3}=*τ*/6} on *y*={0,*R*} and {*q*_{1}=*τ*/2, *q*_{2}=0 and *q*_{3}=*τ*/6} on *x*={0,*R*}. The strong anchoring or Dirichlet conditions can be recovered from (3.17)–(3.19) and (3.20)–(3.22) in the limit, *d*^{(x)}_{e}→0 and *d*^{(y)}_{e}→0, respectively.

We solve the Euler–Lagrange equations (3.14)–(3.16) and the boundary constraints (3.17)–(3.22) using relaxation methods that have been successful in the study of static LC textures with topological defects [13,24]. These methods compute the static solutions, (*q*_{1},*q*_{2},*q*_{3}), by mimicking a dynamic gradient-flow-like procedure along which the total free energy continuously decreases till the equilibrium is attained, for an explicitly prescribed initial condition. We use three different kinds of initial conditions: (i) bulk uniaxial alignment along *e*_{x}, (ii) bulk uniaxial alignment along *e*_{x}+*e*_{y}, or (iii) the isotropic phase with **Q**=0. The **Q**-tensor is recovered from the solution, (*q*_{1},*q*_{2},*q*_{3}), by using the parametrization (3.2). The solutions are robust with respect to different choices of initial conditions and we thus deduce that they are numerically stable.

## 4. Results

Macroscopic wells with *R*≫*ξ*^{(0)}_{b}, e.g. micrometre-sized wells, have already been modelled in detail in [4,18]. We focus on wells with *R* comparable to *ξ*^{(0)}_{b} in this paper. We use fixed values of *L*,*B*,*C* and, hence, the bare biaxial correlation length *η*=*R*/*ξ*^{(0)}_{b}. In simulations, we work with temperatures below *T*_{*}, corresponding to *τ*>2 (see equation (3.10)). Previous work indicates that biaxial textures are more likely in this low-temperature regime, compared with the high-temperature regime with temperatures above the nematic–isotropic transition temperature [22,25]. For example, we carry out illustrative simulations with *τ*=4. The same temperature was chosen in [22], where the authors study OR patterns in ‘classical’ hybrid planar cells.

### (a) Diagonal structure with defects

In figures 1, 2*a* and 3, we consider the specific example of a LC well with *R*=4.5*ξ*^{(0)}_{b}, at temperature *τ*=4 (with *t*=−8 in (3.10)), with strong anchoring conditions on the lateral surfaces. In dimensional terms, this would correspond to a well with *R*∼120–150 nm [19]. The Dirichlet conditions (see (3.7) and (3.8)) induce an alignment mismatch along the four vertical edges and, consequently, we obtain four line defects along the four vertical edges in the *z*-direction. In figures 1 and 2*a*, we plot the spatial profile of *β*^{2}, and in figure 3 we plot the leading eigenvector (with the largest positive eigenvalue) of the computed **Q**-tensor on the cross section *z*=0. This eigenvector evidently has a diagonal profile across the square (*x*,*y*) cross section and is strongly reminiscent of the diagonal solution reported in [4,18]. For a purely uniaxial solution with two degenerate eigenvalues, this eigenvector is simply the director in (2.2). From figure 3, it is clear that the leading eigenvector has four defects at the four square vertices and each defect can be regarded as a quarter of +1 or −1-degree defect, i.e. the defects at the top left and bottom right are a quarter of +1-degree radial defects and the defects at the bottom left and top right are a quarter of −1-degree defects. In other words, the leading eigenvector rotates by 90° between a pair of coincident square edges.

In figures 1 and 2*a*, we study the structural characteristics of each defect core in terms of *β*^{2}; the profile is somewhat similar to the classical OR phenomenon studied in [9], etc. We have two intersecting faces at each vertical edge and these faces enforce strictly uniaxial alignments with mutually perpendicular directors, along *e*_{x} and *e*_{y}, respectively, and positive scalar order parameter. As is evident from the plots, the LC state mediates between the two uniaxial ‘boundary’ states by means of an intermediate biaxial pear-shaped lobe. This biaxial lobe has a rim of maximal biaxiality (*β*^{2}=1) and the interior of the lobe has suppressed biaxiality.

In figure 4, we look at *β*^{2} along a square diagonal on the surface *z*=0. Our model assumes invariance in the *z*-direction and, hence, it suffices to look at structural characteristics on *z*=0. We see that *β*^{2} monotonically increases from the vertex to the rim of the biaxial lobe, where *β*^{2}≈1. Here, *β*^{2} has a local maximum and this local maximum is followed by a local minimum near the centre of the diagonal. At the centre, we are relatively far away from all four square vertices. Then, *β*^{2} increases as we approach the diagonally opposite vertex, attains a local maximum at the rim of the diagonally opposite biaxial lobe and finally decreases to zero as we hit the opposite vertex. The *β*^{2}-profile is characterized by two maxima, one for each square vertex, separated by a local minimum. The biaxial defect lobes are asymmetric in shape, in response to the neighbouring defect cores at the different vertical edges. One could define the characteristic linear defect core size to be the distance from the vertex to the nearest rim of maximal biaxiality, measured along the square diagonal. This length is roughly given by *ξ*_{b}, as is expected from previously reported results in the literature [13,14,21,22]. We refer to this diagonal profile, with biaxial defect cores near the vertical edges, as the DSD in the subsequent text.

### (b) Well order-reconstruction structure

In this section, we study structural transitions induced by gradually decreasing the ratio *η*=*R*/*ξ*^{(0)}_{b} at a fixed temperature *τ*=4. In figure 2*a*–*d*, we demonstrate a sequence of *β*^{2}(*x*,*y*) textures obtained by gradually decreasing the ratio, *η*=*R*/*ξ*^{(0)}_{b}, with strong anchoring conditions on the lateral surfaces. As *η* decreases, the biaxial defect lobes become larger (relative to domain size), overlap and eventually connect to yield a two-dimensional star-shaped rim of maximal biaxiality connecting the four vertices. In particular, in figure 2*d*, the biaxial ring separates uniaxial states with negative order parameter along the square diagonals from positively ordered uniaxial states near the square edges (consistent with the boundary conditions).^{4} In this sense, this is a fully two-dimensional OR pattern connecting the four vertices.

If we follow *β*^{2} along a small circular arc, centred around a square vertex, we expect to see the following sequence of characteristic states: (i) positive uniaxiality (positive scalar order parameter) with *β*^{2}=0 on the square edge, (ii) maximal biaxiality with *β*^{2}=1, (iii) negative uniaxiality (negative scalar order parameter) with *β*^{2}=0, (iv) maximal biaxiality with *β*^{2}=1, and (v) positive uniaxiality with *β*^{2}=0 on co-incident edge. By symmetry, we expect to see state (iii) along the square diagonal. In figure 4, we monitor the structural characteristics along a square diagonal. In figure 4*a*, we plot *β*^{2}(*η*), measured along the square diagonal, for different values of *η*. For relatively large *η*, say *η*=4.5, the *β*^{2}-profile has two distinct maxima. This corresponds to the DSD structure reported above. We refer to the two maxima in the regions *x*/*R*<0.5 and *x*/*R*>0.5 as the *left peak* and *right peak*, respectively. As *η* decreases, the maxima approach each other and coalesce at *η*=3.39. For smaller values of *η*, say *η*=3.3, the *β*^{2}-maximum is suppressed in magnitude, i.e. is smaller than unity. At *η*=*η*_{c}=3.28±0.01, *β*^{2} vanishes along the square diagonal (within numerical resolution) and *β*^{2}=0 along the diagonal for all *η*≤*η*_{c}.

In figure 4*b*, we compute 〈*β*^{2}〉_{d} and *β*^{2}_{m} as a function of *η*, corresponding to the average value and the maximal value of *β*^{2} along the square diagonal, respectively. Additionally, we also plot the position of the left peak, denoted by *x*_{m} in the region *x*∈[0,*R*/2], and the spatial average of *β*^{2} across the square cross section, denoted by 〈*β*^{2}〉, as a function of *η*. From figure 4*b*, we see that the global OR profile is established at *η*=*η*_{c}=3.28±0.01, for which *η*≤*η*_{c}. For *η*≤*η*_{c}, we have an uniaxial LC state, with negative scalar order parameter, along the square diagonal. The corresponding critical well size *β*^{2}〉(*η*)≈0.4 for *η*≤*η*_{c}. As *η*≤*η*_{c}, *x*_{m} is not defined for *η*<*η*_{c}. In what follows, we refer to the star-shaped biaxial OR pattern in figure 2*d* as the well order-reconstruction (WORS). We define the onset of the WORS by the critical value of *η*=*η*_{c} for which *z*=0, such that *β*^{2}_{m}=0 for all *η*<*η*_{c}. Our simulations suggest that the DSD–WORS transition is continuously achieved by decreasing *η* and hence this is a second-order transition.

### (c) Phase diagrams

In figure 5*a*, we numerically study the dependence of the critical well size, *t*=(*T*−*T*_{*})/(*T*_{**}−*T*_{*}), with strong anchoring conditions on the lateral surfaces. As *t* increases towards *T*=*T** (the temperature at which the isotropic phase loses stability), *R*_{c} increases steeply. This steep ascent can be explained on the grounds that the condensation energy penalty for deviations from the uniaxial state with *S*=*S*_{eq}(*T*), decreases as *t* increases. Hence, biaxial OR patterns (such as the WORS pattern) become increasingly energetically viable, and hence observable as *t* increases.

In figure 5*b*, we impose equal anchoring strength on all four lateral surfaces and plot *ξ*^{(0)}_{b}/*d*_{e}, where *t*=−8 and *t*=0. These temperatures correspond to qualitatively different temperature regimes (i.e. *T*≪*T*_{*} and *T*=*T*_{*}<*T*_{IN}). We observe a slowly increasing profile of *R*_{c} versus 1/*d*_{e} for both temperatures. From the reasoning in figure 5*a*, the critical threshold *R*_{c} is always higher for higher temperatures. Further, *R*_{c} rapidly saturates to a limiting value as 1/*d*_{e} increases; this limiting value is simply the strong anchoring limit. For *t*=−8 or equivalently for *τ*=4, this limiting value is *d*_{e} decreases or equivalently as the anchoring strength decreases, the uniaxial boundary conditions are relatively weakly implemented. This effectively increases the lateral well size as the system has anchoring freedom on the lateral surfaces and hence does not need to strictly match the preferred uniaxial states on intersecting edges. The WORS profile ceases to exist for anchoring strengths below a certain critical threshold. For such cases, the anchoring conditions are not strong enough to support geometrically imposed frustration.

### (d) Comparison with classical order reconstruction structure

In figure 6*a*–*c*, we compare the WORS with the classical OR pattern reported in a batch of numerical papers [16,21,22]. In figure 6*a*, we plot *β*^{2}(*x*,*y*) across the square cross section and observe a distinct star-shaped rim of maximal biaxiality with *β*^{2}(*x*,*y*)=1. In figure 6*b*, we plot −*β*^{2}(*x*,*y*) to emphasize the cross-shaped uniaxial ordering (i.e. *β*^{2}=0) along diagonals. In figure 6*b*, the local eigenframe of **Q** coincides with the {*e*_{x},*e*_{y},*e*_{z}} laboratory frame and *q*_{2}=0 throughout the square domain in (3.2).

We can reproduce the classical OR phenomenon in our framework by enforcing conflicting planar and homoeotropic anchoring on the plates at *y*=0 and *y*=*R*_{y}, respectively, while imposing natural boundary conditions on *x*=0 and *x*=*R*_{x}. Equivalently, we require that
*R*_{y}, while keeping *R*_{x} fixed, at constant temperature *τ*=4. At a critical value, *η*^{(CL)}=2.45 (the superscript ^{(CL)} marks the *classical* OR transition), we observe a structural transition into the classical OR pattern and the classical OR pattern is observed for all *η*<*η*^{(CL)}. We plot *β*^{2} for the classical OR case in figure 6*c* and see two sheets of maximal biaxiality (*β*^{2}=1) enclosing a sheet of negative uniaxiality (*β*^{2}=0). On comparing figure 6*a*,*c*, it is clear that the WORS is qualitatively different to the classical OR pattern; the WORS arises from mutually perpendicular uniaxial alignments on four pairs of coincident edges, whereas the classical OR phenomenon arises from mutually perpendicular uniaxial alignments on pairs of parallel edges. While the WORS is observed for *η*<3.28, the classical OR pattern is observed for *η*<2.45, i.e. the WORS has a broader window of stability compared with the classical OR case. We do not have rigorous explanations for these effects but elevated stability thresholds can be potentially attributed to the topological defects present in the WORS. Similar effects have been observed in hybrid cylindrical cells with a boojums on bounding plates [21].

### (e) Interior defects and the WORS pattern

We demonstrate a DSD–WORS structural transition induced by the introduction of locally melted regions within the DSD structure. This could, in practice, be realized by electronic-beam (e-beam) lithography (see http://en.wikipedia.org/wiki/Electron-beamlithography) [15]. E-beam lithography techniques can fabricate sub-10 nm cavities/indentations within square or cuboid geometries. Alternatively, one could introduce a nanoparticle (NP) treated to impose localized isotropic melting at the NP–LC interface. We numerically introduce a melted square cavity of lateral size, *r*_{m}, at the centre of the well. This requires **Q**(*R*/2±*r*_{m}/2, *R*/2±*r*_{m}/2,*z*)=0. We take *τ*=4, *R*/*ξ*^{(0)}_{b}=4.5, with strong anchoring conditions as in (3.7) and (3.8) and no external fields. In dimensional terms, this describes a well with lateral dimensions between 120 and 150 nm [19]. In figure 7*a*–*c*, we numerically compute LC equilibria in the absence and presence of locally melted regions, respectively. We recover the familiar DSD in figure 7*a*. However, in figure 7*b*, we observe a star-shaped rim of maximal biaxiality, characteristic of the WORS. In figure 7*c*, we displace the melted region from the centre towards the lower plate at *y*=0 and lose the WORS. This suggests that the WORS can be efficiently stabilized if the locally melted region is compatible with the structural symmetry of the WORS.

In figure 8, we study the effect of the size of the melted region, *r*_{m}, on the DSD–WORS structural transition. We follow 〈*β*^{2}〉 as a function of *η*, for different values of *r*_{m}. The curves are almost stationary till they hit the critical value of *η*=*η*^{(m)}_{c} and the WORS is stable for all *η*≤*η*^{(m)}_{c}. It is evident that interior melted regions substantially increase the value of *η*^{(m)}_{c}. For example, we obtain *η*_{c}=3.28±0.01 for *r*_{m}=0 (no melted region) and *η*^{(m)}_{c}=4.67±0.01 for *r*_{m}/*R*=0.1. Hence, the threshold value is increased by roughly 42%.

## 5. Conclusion

We model and numerically analyse LC textures inside shallow square wells, within the LdG theory for nematic LCs. These wells were first reported in [18] and modelled in the papers [4,18]. In [4,18], the authors study large micrometre-sized wells and perform a strictly two-dimensional study on a square or a rectangle. Our aim is to quantify the effects of confinement, temperature, boundary conditions on the experimentally observed well textures, with special emphasis on new biaxial textures. Our work purely focuses on static LC textures and we do not study field-induced or dynamic effects.

We study shallow square wells with cross-sectional dimension *R*, where *R* is measured in units of the bare biaxial correlation length, *ξ*^{(0)}_{b}. We adopt a three-dimensional parametrization of the LdG **Q**-tensor (as opposed to two-dimensional parametrization used in previous work) and assume that all variables only depend on the spatial coordinates, (*x*,*y*), and are independent of the *z*-coordinate. The boundary conditions on the lateral surfaces, *x*={0,*R*} and *y*={0,*R*}, are imposed either via Dirichlet conditions (as in (3.7) and (3.8)) or via a surface anchoring energy as in (2.8). The Dirichlet conditions in (3.7) and (3.8) create line defects along the vertical edges. For wells with *R*≫*ξ*^{(0)}_{b}, the static structures are effectively uniaxial away from the edges and are well described by the *diagonal* and *rotated* solutions in [4,18]. We study wells with *R*=*O*(*ξ*^{(0)}_{b}) in this paper and obtain two distinct LC states: the DSD and WORS patterns. The DSD pattern is qualitatively similar to the diagonal solution reported in [4,18] and is only obtained for relatively large wells with *R* greater than a temperature-dependent critical length, denoted by *R*_{c}. The DSD is largely uniaxial away from the vertical line defects, the principal eigenvector has a diagonal profile on the bottom surface (*z*=0) and there are localized pear-shaped biaxial neighbourhoods near each vertical edge, the rims of which exhibit maximal biaxiality. We believe that the macroscopic properties of the DSD solution are qualitatively similar to the diagonal solution in [4,18] and our three-dimensional modelling recovers the biaxial structure of the defect cores, which is expected but not reported in previous work.

We progressively decrease the ratio *η*=*R*/*ξ*^{(0)}_{b} and at a critical value *η*=*η*_{c}(*τ*) (e.g. *η*_{c}(4)=3.28±0.01), we observe a global WORS with a star-shaped ring of maximal biaxiality connecting the four vertices of the square cross section on *z*=0. This ring is repeated throughout the height of the cell. The biaxial ring separates uniaxial states with negative order parameter along the square diagonals from uniaxial states with positive order parameter outside the ring, respectively. The critical value depends on the temperature and on the anchoring strength on the lateral surfaces. We define a quantitative criterion for the onset of the WORS pattern, in terms of the maximum value of *β*^{2} (*η*, denoted by *η*_{c}, for which *η*≤*η*_{c}. The DSD–WORS transition is a continuous transition and we interpret stable states with *η*>*η*_{c} as being a variant of the DSD state. Further, we study the dependence of *η*_{c} on the temperature and the anchoring strength on the lateral surface, with no external fields in figure 5. Numerical simulations show that interior locally melted regions, do not substantially alter the symmetry of WORS but substantially increase the stability window of the WORS.

Next, we discuss the reliability of the mesoscopic LdG theory for sub-micrometre and nano-scale systems, such as the model systems studied here. Firstly, the LdG theory has been used, with success, to study nano-scale defect structures [9] and the OR phenomenon in nano-confined systems [13,14,,16,22,26]. Recent experimental work in [26,27,28] confirms that mesoscopic modelling yields surprisingly good predictions, even on the nanometre scale. For example, a simplified LdG-type model can reproduce experimental observations in cylindrical nano-channels with radius of order approximately 5 nm [27]. Secondly, our model predicts that the WORS is observable for *R*∼3.28*ξ*^{(0)}_{b} for temperature, *τ*=4. As stated in the Introduction, *ξ*^{(0)}_{b} is typically tens of nanometre [11]. Recent experimental work in a *π*-cell geometry, with conventional optical microscopy techniques, suggests that *ξ*^{(0)}_{b}∼33±9 nm [19]. This would correspond to cell sizes, *R*∼78−138 nm, and such dimensions are widely modelled by LdG approaches in the community. At higher temperatures, the WORS is observable in even larger wells, with dimensions *R*∼240−420 nm (figure 5*a*). Furthermore, polymeric LCs can have correlation lengths which are an order of magnitude larger than that of their conventional nematic counterpart [3]. In such cases, the WORS would be observed even in micrometre-scale wells.

Biaxial patterns have been extensively found in the literature, see [6,9,13,14,16,21], etc., the reported textures are heavily localized near defects or boundaries. The WORS is a new global two-dimensional biaxial texture predicted for sub-micrometre and nano-systems with moderate to strong anchoring conditions, for temperatures below the supercooling temperature. The WORS is qualitatively different to the reported diagonal and rotated solutions in [4,18]. As such, we expect that it will offer different optical properties and new responses to electric fields, etc. Our numerical work can give engineers quantitative information about the scale and shape of prototype devices that can support the WORS. Once such devices are actually engineered, one could experimentally measure the electro-optical responses of the WORS, possibly by analogy with similar experimental work on OR patterns in one-dimensional hybrid cells in [12,17]. While our work sets up a sound foundation for the study of the WORS, there are several directions for further study which include (i) field-induced effects, (ii) dynamic phenomena, (iii) effects of elastic anisotropy, (iv) non-cuboid geometries, and (v) inclusion of interior defects. A comprehensive study on these lines can yield vital information about systems with intrinsic biaxiality and how biaxiality can be exploited for new nano-scale applications.

## Funding statement

A.M. is supported by an EPSRC Career Acceleration Fellowship EP/J001686/1, an Oxford Centre for Collaborative Applied Mathematics (OCCAM) Visiting Fellowship and a Keble Research Grant, Keble College.

## Acknowledgements

A.M. and S.K. thank the Isaac Newton Institute, University of Cambridge for a Visiting Fellowship as part of the ‘Mathematics of Liquid Crystals’ programme. S.K. thanks OCCAM and the University of Bath for supporting collaborative trips in connection to this project. A.M. and S.K. thank all referees for constructive suggestions.

## Footnotes

↵1 The trends would be unchanged for heights

*R*/10≤*h*≤*R*/5.↵2 The tangent conditions may be violated for energetic considerations.

↵3 We could also use weak anchoring or surface energies on

*z*=0 and*z*=*h*to obtain similar results.↵4 We interpret regions with

*β*^{2}<0.01 as being uniaxial for practical purposes. We note that*β*^{2}may not, strictly speaking, vanish at any point inside the well.

- Received April 4, 2014.
- Accepted May 28, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.