Elastic anisotropy in the reduced Landau-de Gennes model.

We study the effects of elastic anisotropy on Landau-de Gennes critical points, for nematic liquid crystals, on a square domain. The elastic anisotropy is captured by a parameter, L 2 , and the critical points are described by 3 d.f. We analytically construct a symmetric critical point for all admissible values of L 2 , which is necessarily globally stable for small domains, i.e. when the square edge length, λ , is small enough. We perform asymptotic analyses and numerical studies to discover at least five classes of these symmetric critical points-the WORS , Ring ± , C o n s t a n t and p W O R S solutions, of which the WORS , Ring + and C o n s t a n t solutions can be stable. Furthermore, we demonstrate that the novel C o n s t a n t solution is energetically preferable for large λ and large L 2 , and prove associated stability results that corroborate the stabilizing effects of L 2 for reduced Landau-de Gennes critical points. We complement our analysis with numerically computed bifurcation diagrams for different values of L 2 , which illustrate the interplay of elastic anisotropy and geometry for nematic solution landscapes, at low temperatures.

Introduction

Nematic liquid crystals (NLCs) are quintessential examples of partially ordered materials that combine fluidity with the directionality of solids [1]. The nematic molecules are typically asymmetric in shape, e.g. rod- or disc-shaped, and these molecules tend to align along certain locally preferred averaged directions, referred to as nematic directors in the literature. Consequently, NLCs have a long-range orientational order and direction-dependent physical, optical and rheological properties. It is precisely this anisotropy that makes them the working material of choice for a range of electro-optic devices such as the multi-billion dollar liquid crystal display industry [2,3].

There has been substantial recent interest in multistable NLC systems, i.e. NLCs, confined to two-dimensional (2D) or three-dimensional (3D) geometries that can support multiple stable states without any external electric fields [4-10]. Multistable NLC systems offer new prospects for device technologies, materials technologies, self-assembly processes and hydrodynamics. This paper is motivated by a bistable system reported in [11]. Here, the authors experimentally and numerically study NLCs inside periodic arrays of 3D wells, with a square cross-section, such that the well height is typically much smaller than the square cross-sectional length. Furthermore, the authors speculate that the structural characteristics are translationally invariant along the well-height, effectively reducing this to a 2D problem. Hence, the authors restrict attention to the bottom square cross-section with square edge length denoted by , which is typically on the micron scale. The choice of boundary conditions is crucial and in [11], the authors impose tangent boundary conditions (TBCs) on the well surfaces, i.e. the nematic directors, in the plane of the well surfaces, are constrained to be in the plane of the surfaces. However, this necessarily means that the nematic director is tangent to the square edges, creating defects at the vertices, where the director is not defined. The authors observe two classes of stable NLC states: the diagonal states, for which the nematic director aligns along one of the square diagonals and; the rotated states, for which the director rotates by radians between a pair of opposite square edges.

In [6,12], the authors model this square system within the celebrated continuum Landau–de Gennes (LdG) theory for NLCs. The LdG theory describes the nematic state by a macroscopic order parameter—the -tensor-order parameter [1]. From an experimental perspective, the -tensor is measured in terms of NLC responses to external electric or magnetic fields, which are necessarily anisotropic in nature. Mathematically, the -tensor-order parameter is a symmetric, traceless, matrix with 5 d.f. For a square domain with TBCs on the square edges, it suffices to work in a reduced LdG framework where the -tensor only has three degrees of freedom, , , . The degree of nematic order in the plane is captured by and , whereas measures the out-of-plane order, such that positive (negative) implies that the nematic director lies out of the plane (in the plane) of the square, respectively. The TBCs naturally constrain to be negative on the square edges, but could be positive in the interior.

The LdG theory is a variational theory, i.e. experimentally observable states can be modelled by local or global minimizers of an appropriately defined LdG free energy. In the simplest setting, the LdG energy has two contributions—a bulk energy and an elastic energy that penalizes spatial inhomogeneities. In these papers, the authors work with low temperatures that favour an ordered nematic state. The elastic energy is typically a quadratic and convex function of and, in [6,12], the authors work with an isotropic elastic energy—the Dirichlet elastic energy. In a reduced LdG setting, the authors recover the stable and states for large and, in [12], they discover a novel stable Well Order Reconstruction Solution () for small . The is special because it exhibits a pair of mutually orthogonal defect lines, with no planar nematic order along the defect lines as will be described in §3. In [13], the authors generalize this work to arbitrary 2D regular polygons and, in [8], the authors study 3D wells, with an emphasis on novel mixed solutions which interpolate between two distinct solutions.

In this paper, we study the same problem of NLCs on a square domain with TBCs, with an anisotropic elastic energy as opposed to the isotropic energy studied in [6,12]. Notably, we take the elastic energy density to be , where is the anisotropy parameter. Physically speaking, positive implies that splay and bend deformations of the nematic director are energetically expensive compared to out-of-plane twist deformations, i.e. we expect the physically observable states to have positive in the square interior, as increases. Therefore, there are competing effects of the TBCs on the square edges, which prefer in-plane director orientation, and the preferred out-of-plane director orientation in the square interior, for larger values of . We construct a symmetric critical point of the LdG energy, for any . This symmetric critical point is globally stable for small enough. The is a special case of this symmetric critical point with on the square domain, for . For , the WORS does not survive with the perfect cross symmetry along the square diagonals. We perform an asymptotic analysis for small and small , about the . The anisotropy has a first-order effect on , i.e. is perturbed linearly by , and increases at the square centre for positive , relative to its value for , corroborating the trend of increasing with increasing . The globally stable symmetric critical point for small and small , labelled as the solution, effectively smoothens out the and exhibits a stable central -degree point defect. We perform formal calculations to show that as , energy minimizers (and consequently the symmetric critical point described above for small ) approach the state with , away from the square edges and exhibits four boundary layers near the edges. Thus, there are three different classes of the symmetric critical point discussed above: the , which only exists for ; the solution, which only can be stable for moderate values of and non-zero and; the solution, which exists for large enough and is always stable according to our heuristics and numerical calculations. Additionally, we also find two unstable classes of this symmetric critical point, both of which exist for moderate values of and . These are the solution which exhibits a central -degree point defect, and the novel which exhibits an oscillating sequence of nematic point defects along the square diagonals. We provide asymptotic approximations for the novel solution branch.

While most of our work is restricted to the small -limit, we also touch on energy minimizers in the limit. The competitors in the large -limit are the familiar and states, and the solution. Using Gamma-convergence arguments, we show that the solution has lower energy than the and states, for large enough . We complement our analysis with numerical computations of bifurcation diagrams for five different values of . To summarize, our notable findings on the response of the NLC solution landscape for this model problem, to the elastic anisotropy parameter, , are (i) novel stable states ( and ) for small , and (ii) enhanced multistability in the large -limit due to the competing and states, for large . As increases, we expect that there are further, not necessarily energy-minimizing, LdG critical points with positive , or out-of-plane nematic directors in the square interior. Furthermore, has a stabilizing effect with respect to certain classes of planar perturbations and out-of-plane perturbations, and so we expect enhanced multistability as increases, for all values of .

A lot of open questions remain with regards to the interplay between , and temperature on NLC solution landscapes, but our work is an informative forward step in this direction. Our paper is organized as follows. We provide all the mathematical preliminaries in §2. We construct the symmetric critical points described above in §3. In §4, we perform separate asymptotic studies in the small and small limit; large limit; large -limit. In §5, we present bifurcation diagrams for five different values of , accompanied by some rigorous stability results. We conclude with some perspectives in §6.

Preliminaries

In this section, we review the LdG theory of NLCs. Within this framework, the nematic state is described by a macroscopic LdG order parameter—the -tensor-order parameter. The -tensor is a symmetric, traceless, matrix, which is a macroscopic measure of the nematic anisotropy. The eigenvectors of represent the preferred material directions, the corresponding eigenvalues measure the degree of order about these directions. The nematic director is often identified with the eigenvector that has the largest positive eigenvalue. The -tensor is said to be: (i) isotropic if ; (ii) uniaxial if has a pair of degenerate non-zero eigenvalues; and (iii) biaxial if has three distinct eigenvalues. A uniaxial -tensor can be written as , where is the identity matrix, is the distinguished eigenvector with the non-degenerate eigenvalue, and is a scalar order parameter. The LdG theory is a variational theory with an associated free energy, and the basic modelling hypothesis is that the physically observable configurations correspond to global or local energy minimizers subject to imposed boundary conditions. We work with 2D domains, , in the context of modelling thin 3D systems. In the absence of a surface anchoring energy and external fields, the LdG free energy is given by where and are the elastic and thermotropic bulk energy densities, respectively. We consider a two-term elastic energy density: where is an elastic constant, and is the ‘elastic anisotropy’ parameter. The elastic energy density penalizes spatial inhomogeneities, typically quadratic in . In terms of notation, we use and , , where the Einstein summation convention is assumed throughout this manuscript. We work with a 2D confining geometry in this paper and hence, for all . We work with the simplest form of , that allows for a first-order isotropic–nematic transition: Here, , and , for . We take to be the rescaled temperature and are material-dependent constants. In this regime, is the absolute temperature, and is the characteristic nematic supercooling temperature. The rescaled temperature, , has three physically relevant values: (i) , below which the isotropic state, , loses stability; (ii) the nematic super-heating temperature , above which is the unique critical point of and (iii) the nematic-isotropic phase transition temperature , at which is minimized by the isotropic phase and a continuum of uniaxial states. We work with low temperatures, , for which is minimized on the set of uniaxial -tensors defined by where is the space of traceless symmetric matrices and We non-dimensionalize the system using, , where is a characteristic geometrical length-scale, e.g. half edge length of a 2D regular polygon. The rescaled LdG energy functional (up to a multiplicative constant) is given by where is the rescaled domain in , and is the rescaled area element. We drop the ‘bars’ but all computations should be interpreted in terms of the rescaled variables.

Next, we define the working domain and Dirichlet boundary conditions, although we believe that our methods can be generalized to arbitrary 2D domains. We focus on square domains, building on the substantial work in [7,14,15]. We impose almost uniaxial Dirichlet TBCs on the square edges, which require the nematic director to be tangent to the edges, necessarily creating a mismatch at the square vertices. This is consistent with the experimentally and numerically investigated TBCs, [6,11,12]. To avoid the discontinuities at the vertices, we take to be a truncated square whose edges are parallel to the coordinate axes: Provided , the truncation does not change the qualitative properties of the LdG energy minimizers away from the square vertices. The boundary, , has four ‘long’ edges parallel to the coordinate axes, defined in a clockwise fashion as , where lies parallel to the -axis at . The truncation creates four additional ‘short’ edges, of length , parallel to the lines and , labelled as in a clockwise fashion, starting at the top-left corner of the domain. The domain is illustrated in figure 1. In particular, we fix the uniaxial director to be on the edges, and , and on and . From a physical standpoint, this constitutes strong (infinite) anchoring on the long edges. One could also model weak (finite) anchoring condition with an additional surface energy in the LdG free energy [16], but that would make the analysis more complicated. We set where where and are unit vectors in the - and -directions, respectively. In particular, on . On the short edges, , we prescribe a continuous interpolation between the boundary conditions on the long edges (2.8) given by where is a unit vector in the -direction, and is a smoothing function, i.e. Although the boundary conditions (2.9) do not minimize on , and do not respect TBCs, they are short by construction and are chosen purely for mathematical convenience. Given the Dirichlet boundary conditions (2.8) and (2.9), the admissible space is The energy minimizers, or indeed any critical point of the LdG energy (2.5), are solutions of the associated Euler–Lagrange equations: which comprise a system of five nonlinear coupled partial differential equations. The terms and are Lagrange multipliers associated with the tracelessness constraint.

Figure 1

The truncated square domain .

Finally, we comment on the physical relevance of the 2D domain, . Consider a 3D well, where , and is a characteristic length scale of . In this limit, one can assume (at least for modelling purposes) that physically relevant -tensors are independent of the -coordinate, i.e. the profiles are invariant across the height of the well, and that is a fixed eigenvector (see [17,18] for some rigorous analysis and justification). This implies that we can restrict ourselves to -tensors with three degrees of freedom: subject to the boundary conditions where , on ; , on ; , on ; , on ; , on ; , on , and The conditions (2.13) and (2.14) are equivalent to Dirichlet conditions in (2.7).

Qualitative properties of equilibrium configurations

In [12], the authors numerically compute critical points of (2.5) with , satisfying the Dirichlet boundary conditions (2.7), on the square cross-section . For the edge length, small enough, the authors report a new . The has a constant set of eigenvectors, and , which are the coordinate unit vectors. The is further distinguished by a uniaxial cross, with negative scalar order parameter, along the square diagonals. Physically, this implies that there is a planar defect cross along the square diagonals, and the nematic molecules are disordered along the square diagonals. In [7], the authors analyse this system at a fixed temperature with , and show that the is a classical solution of the associated Euler–Lagrange equations (2.11) of the form: The single degree of freedom, , is a solution of the Allen–Cahn equation with the following symmetry properties: Notably, everywhere for the (refer to (2.12)), which is equivalent to having a set of constant eigenvectors. They prove that the is globally stable for small enough, and unstable for large enough, demonstrating a pitchfork bifurcation in a scalar setting. Their analysis is restricted to the specific temperature and, in [8], the authors extend the analysis to all , with . In this section, we analyse the equilibrium configurations with , including their symmetry properties in the small limit. Notably, we show that the cross structure of the does not survive with .

Proposition 3.1.

There exists at least one solution to the Euler–Lagrange equations () of the form () in , given the Dirichlet boundary conditions () and (). For this solution, the functions  satisfy the following systems of PDEs: and the boundary conditions (2.13) and (2.14).

Proof.

Our proof is analogous to theorem 2.2 in [18]. Substituting the -tensor ansatz (2.12) into the general form of the LdG energy (2.5), let where and are the elastic and thermotropic bulk energy densities, respectively. We prove the existence of minimizers of in the admissible class which will also be solutions of (2.11) in the admissible space, . Since the boundary conditions (2.13) and (2.14) are piece-wise of class , is non-empty. If , is in the form of (3.7). If , the elastic energy density can be rewritten as a function of in the following way: The difference between the expressions for in (3.7) and (3.10), is a null Lagrangian, and hence can be ignored with the Dirichlet boundary condition. Since we assume that , the elastic energy density is the sum of non-negative terms for any and, more specifically, Furthermore, also satisfies for some constant, , depending only on and . Hence is coercive in . Finally, we note that is weakly lower semi-continuous on , which follows immediately from the fact that is quadratic and convex in . Thus, the direct method in the calculus of variations yields the existence of a global minimizer of , among the finite energy triplets , satisfying the boundary conditions (2.13) and (2.14) [19]. One can verify that the semilinear elliptic system (3.3)–(3.5) are the Euler–Lagrange equations associated with , and the minimizers for are solutions of (3.3)–(3.5). The corresponding -tensor (2.12) is an exact solution of the LdG Euler–Lagrange equations (2.11).

Proposition 3.2.

There exists a critical point  of the energy functional (3.6) in , for all , such that  on the square diagonals  and , and  on  and .

We follow the approach in [7]. Consider of a square located in the positive quadrant of : The following boundary conditions on are consistent with the boundary conditions (2.13) and (2.14) on the whole of : where represents the outward normal derivative. We minimize the associated LdG energy in , given by on the admissible space As the boundary conditions on are continuous and piecewise of class , is non-empty. Furthermore, we have shown that is coercive on and convex in the gradient . Thus, by the direct method in the calculus of variations, there exists a minimizer . We define a function by odd reflection of about the square diagonals, and even reflection about - and -axes. We do the same for defined by even reflections of about the square diagonals, and odd reflection about - and -axes and lastly, for the function defined by even reflections of about the square diagonals and - and -axes. By repeating arguments in [20], the new triple, , is a weak solution of the Euler–Lagrange equations on . One can verify that is a critical point of on with the desired properties.

Proposition 3.3.

For  and , the critical point constructed in proposition 3.2, denoted by , has non-constant  on , for all .

We proceed by contradiction. Assume that is constant on . Recalling the boundary conditions (2.14), we necessarily have that in . Substituting into (3.4), we obtain for arbitrary real-valued functions , with on . Therefore, in . Substituting and into (3.3) and (3.5) yields and where From the reduced PDEs for , (3.16) and (3.17), one can calculate Furthermore, from the symmetry properties of the constructed solution in proposition 3.2, we have Substituting (3.22) into (3.21), we obtain If , then equation (3.23b) at is non-zero, which leads to a contradiction. If , then and (3.17) reduces to and hence, for arbitrary real-valued functions .

From proposition (3.2), we know that for any , satisfies the symmetry property and hence, Subtracting on both sides of the equality (3.24), we get for some constant . The function may now be rewritten as This formulation cannot be extended continuously on the boundary since, for , and , we have which again leads to the required contradiction.

Proposition 3.4.

There exists a critical edge length  such that, for any , the critical point, , in proposition 3.2 is the unique critical point of the LdG energy ().

We adapt the uniqueness argument in lemma 8.2 of [21]. Let be a global minimizer of in (3.6), for . Let be such that a.e. . Defining , where is constant, we have for some constant, , depending only on and . Thus, we restrict ourselves to the following admissible space of -tensors: The second derivatives of are quadratic polynomials in . By an application of the relevant embedding theorem in [22] (theorem 9.16 which implies that for a bounded domain with Lipschitz boundary, for any , , with constant depending only on ), there exist some constant , depending only on and , such that We apply the Hölder inequality to get, for any , Therefore, for any , we have where . Using (3.11), an application of the Poincaré inequality, and repeating the same arguments as above, we have for some constant, , depending only on and the sign of . Using both (3.33a,) and (3.34a,), we have for . Thus, is strictly convex for the finite energy triplets , and has a unique critical point for . Hence, the critical point constructed in proposition 3.2 is the unique minimizer of and, in fact, the unique global LdG energy minimizer (when we consider -tensors with the full 5 d.f.), for sufficiently small .

Lemma 3.5.

Let  be the unique global minimizer of the energy (3.6), for  given by proposition 3.4. Then for any , the function  vanishes along the square diagonals  and , and the function  vanishes along  and .

This is an immediate consequence of proposition 3.2, but we present an alternative short proof based on symmetry. Suppose that is a global minimizer of , in for a given . Then is a solution of the Euler–Lagrange system (3.3)–(3.5), subject to the boundary conditions (2.13) and (2.14). So are the triples that are compatible with the imposed boundary conditions. We combine this symmetry result with the uniqueness result in proposition 3.4 to get the desired conclusion. For example, use with to deduce that . Also, with yields that . Furthermore, we use the relation with to deduce , and similarly, with to obtain .

As in [8,15], we refer to the following dimensionless parameter in our numerical simulations: In figure 2, we plot the unique stable solution of (3.3)–(3.5) with , for . In this figure, and all subsequent figures, we fix with and . When , the solution is the defined by (3.1). When , and 10, and are non-constant as proven above. One can check that vanishes along the square diagonals and , and the function vanishes along and , as proven in lemma 3.5. When and 10, we observe a central -point defect in the profile of , and we label this as the solution. We then perform a parameter sweep of , from 5 to 500, and find one of the symmetric solution branches in proposition 3.2, which is a continuation of the branch. The solutions with are plotted in figure 3. When , we find the for all . When , the solution exhibits a -defect at the square centre, continued from the branch and hence, we refer to it as the solution. When is positive and moderate in value, we again recover the solution branch and the corresponding at the square centre for negative , but for positive . When is large enough, we recover a symmetric solution which is approximately constant, , away from the square edges, shown in the fourth column of figure 3 for . We refer to this solution as the Constant solution throughout this manuscript.

Figure 2

The unique stable solution of the Euler–Lagrange equations (3.3)–(3.5), with , and (from the first to fourth row) , , and , respectively. In the first column, we plot the scalar order parameter by colour from blue to red, and the director profile in terms of white lines. The and profiles are plotted in the second to fourth columns, respectively. The same convention is used throughout the paper. (Online version in colour.)

Figure 3

A solution branch for the system (3.3)–(3.5) with , and , , and from left to right. This solution branch is a symmetric solution branch, as described in proposition 3.2, is unstable for , and , and stable for . We plot and in the first row. and in the second row. (Online version in colour.)

Asymptotic studies

In the following, we work on a square domain without truncated vertices. For a truncated domain, we keep fixed, with short edges of length , where is sufficiently small.

The small and small anisotropy, limit

We work at the special temperature to facilitate comparison with the results in [23], where the authors investigate solution landscapes with . Notably, for and , reduced LdG solutions have for our choice of TBCs on 2D polygons, and it is natural to investigate the effects of the anisotropy parameter, , in these 2D frameworks. At , the governing Euler–Lagrange equations are given by the following system of partial differential equations: satisfying , and on . We take a regular perturbation expansion of these functions in the limit. The leading-order approximation is given by the , , where is a solution of the Allen–Cahn equation, as in [7]: We may assume that can be expanded in powers of as follows: for some functions which vanish on the boundary. For small enough, one can show that there exists a unique solution with and the symmetry property , i.e. on diagonals. Hence, for small enough, the cross structure of the is lost mainly because of effects of on the component , as we discuss below.

From [24], the solutions of (4.3) with , are a good approximation to the solutions of (4.3) for sufficiently small . When , where The analytical solution of (4.5) is given by [25]: The formula will also hold on a truncated square, with Dirichlet conditions on the truncated edges extracted from the explicit formula for , i.e. we can choose boundary conditions on the short truncated edges that are compatible with , once we compute on the full square domain without truncations. When , the unique solution of in (4.4) is , and where with on . See figure 4 for a numerical comparison between the asymptotic solution and relevant solutions of the Euler–Lagrange equations. The asymptotic solution in (4.3) is a good approximation of the solution of the Euler–Lagrange equations in (3.3)–(3.5) when is small enough (see figure 4).

Figure 4

The difference between the solution of (3.3)–(3.5) with , and asymptotic solution . (Online version in colour.)

Proposition 4.1.

(Proof in electronic supplementary material) The analytical solution of (4.7) is given by where  is positive.

The limit

Consider a regular perturbation expansion, in powers of , of the solutions, , of the Euler–Lagrange system (3.3)–(3.5), subject to (2.13) and (2.14). Let be the leading-order approximations of , respectively, in the limit. Then we have:

Proposition 4.2.

The leading-order system in the  limit, given by (4.9)–(4.11), is not an elliptic PDE system.

Proof.

The system of equations (4.9)–(4.11) can be written as where and The system is said to be elliptic, in the sense of I.G. Petrovsky [26], if the determinant for any real numbers . We can check that for any real numbers . Hence, the limiting problem (4.9)–(4.11) is not an elliptic problem.

Proposition 4.3.

There is no classical solution of the limiting problem (4.9)–(4.11), with the boundary conditions (2.13) (in the  limit) and (2.14), where  is the short edge length of the truncated square.

As , the minimizers of in (3.6), with as in (3.7), are constrained to satisfy subject to the Dirichlet TBCs (2.13) and (2.14). Up to , this corresponds to the following PDEs for the leading-order approximations : and almost everywhere, subject to the same TBCs, , , on . As , the boundary conditions for are piecewise constant, and hence the tangential derivatives of and vanish on the long square edges. On , the tangential derivative , hence we obtain in (4.12). Similarly, we have on . This implies that on , where is the outward pointing normal derivative, and we view equation (4.10) to be of the form By the Hopf lemma, when on the boundary, we have . Following the same arguments as in proposition 3.3, this requires that . Substituting into equations (4.12) and (4.13), we obtain , contradicting the boundary condition (2.13). Hence, there are no classical solutions of the system (4.9)–(4.11).

Although there is no classical solution of (4.9)–(4.11) subject to the imposed boundary conditions, we can use finite difference methods to calculate a numerical solution, see figure 5. We label this solution on as the Constant solution, where and are discontinuous on . We now give a heuristic argument to explain the emergence of the Constant solution, as . With constant in the interior of , we have in up to . The choice of constant value is determined by the boundary conditions to minimize the elastic energy . Therefore, , is the unique stable solution of (4.9)–(4.10), except for zero measure sets and we label as the physically relevant Constant solution in the limit. This is consistent with the numerical results in figure 5.

Figure 5

Solutions, , of the leading-order system (4.9)–(4.11) in the limit. (Online version in colour.)

The set of minimizers of, , in the -plane can be written as , where The limit is equivalent to the vanishing elastic constant limit, and converges uniformly to its minimum value in this limit [27].

Proposition 4.4.

Let  be a simply connected bounded open set with smooth boundary. Let  be a global minimizer of  in the admissible class  in (3.9), when . Then there exists a sequence  such that  strongly in  where . If , i.e. then  is a minimizer of in the admissible class . The boundary condition, , is compatible with  on  by the relation . Otherwise, , i.e.

Our proof is analogous to lemma 3 of [27]. See electronic supplementary material.

This proposition cannot be applied to square domains directly, because cannot be defined at the square vertices. The TBCs necessarily mean that is constant on each square edge, and hence, discontinuous at the vertices. However, we can still use the Proposition above to understand the qualitative properties of energy minimizers in the limit, by smoothening the boundary near the vertices and by defining appropriately. Firstly, the TBCs imply that must be a multiple of on the horizontal edges, and an odd multiple of on the vertical square edges. Secondly, we prescribe so that the degree of is zero on the square boundary. For example, experiments suggest that there are two classes of stable equilibria, which are almost in the set , for large —the diagonal and rotated states. The diagonal states, , are such that the nematic director (in the plane) is aligned along one of the square diagonals. The rotated states, labelled as states, are such that the director rotates by radians between a pair of opposite square edges. There are two rotationally equivalent states, and four rotationally equivalent states. The corresponding boundary conditions in terms of are given by or , respectively, where These conditions can be translated to Dirichlet conditions for on the truncated square as well; for example, we can solve on a square domain, subject to these boundary conditions and use this solution to prescribe the Dirichlet conditions on the short edges of the truncated square.

In figure 6, we study the effect of increasing on a state with . When , we see that , almost everywhere. In [13], the authors show that the limiting profiles described in proposition 4.4 are a good approximation to the solutions of (3.3)–(3.5), for large . The differences between the limiting profiles and the numerically computed solutions concentrate around the vertices, for large . As increases, deviates significantly from the limiting value , near the square vertices. From an optical perspective, we expect to observe larger defects near the square vertices for more anisotropic materials with , on large square domains.

Figure 6

The solution with , with and , respectively. In the first row, we plot and . In the second row, we plot the -profiles. (Online version in colour.)

In [25], the authors compute the limiting energy, , of and solutions on a unit square to be: and respectively, where Since is positive, we have . The numerical values of and are approximately zero, so and are approximately for small , and the limiting energies are linear in .

The solution has transition layers on the boundary from to or . The limiting energy in (4.16) is the same for the and cases. Therefore, there are no additional complexities from the term. Analogous to section 4 of [15], using classical arguments in the theory of -convergence, the limiting energy of the solution is the sum of four transition costs: , where is the geodesic distance between and associated with the Riemannian metric , where . The numerical value of is in [15]. The limiting energy is independent of . Hence, there is a critical value such that for , the limiting solution is energetically preferable to the and solutions, i.e. .

The Novel pWORS solutions

For all and , the is a solution of (3.3)–(3.5) given by , where satisfies (4.3). In §4a, we study the Euler–Lagrange equations, in the small and small limit, up to ; see (4.4). However, is a solution for all and we assume that the solution, , of (3.3)–(3.5), can be expanded as follows: In figure 7, we plot a branch of the solutions. As increases, we observe an increasing number of zeroes on the square diagonals, where . For any , we can use the initial condition to numerically find a new branch of unstable solutions, referred to as configurations in figure 7. In the plane, the has a constant set of eigenvectors away from the diagonals, and has multiple -point defects on the two diagonals, so that the is similar to the away from the square diagonals. As increases, the number of alternating and point defects on the square diagonals increases, for the numerically computed . This is mirrored by the function that encodes the second-order effect of on the .

Figure 7

Left three figures: plots of in (4.22), with , 100 and 500 from the left to right. Right two figures: the configurations of the numerically computed with , and 1000. (Online version in colour.)

Bifurcation diagrams

We use the open-source package FEniCS [28] to perform all the finite-element simulations, numerical integration and stability checks in this paper [29,30]. We apply the finite-element method on a triangular mesh with mesh-size , for the discretization of a square domain. The non-linear equations (3.3)–(3.5) are solved by Newton’s methods with a linear LU solver at each iteration. The tolerance for convergence is set to . We check the stability of the numerical solution by numerically calculating the smallest real eigenvalue of the Hessian matrix of the energy functional (3.6), using the LOBPCG (locally optimal block preconditioned conjugate gradient) method [31]. If the smallest real eigenvalue is negative, the solution is unstable, and stable otherwise. In what follows, we compute bifurcation diagrams for the solution landscapes, as a function of , with fixed temperature , for five different values of . The and are fixed material-dependent constants, so is proportional to and we will use these diagrams to infer qualitative solution trends. in terms of the edge length, .

For small enough, there is a unique solution for any value of ; see the results in §4. For , the unique stable solution, for small enough , is the . The unique solution deforms to the , with a central point defect, for relatively small such as . For , the unique solution is the Constant solution, on the grounds that this solution approaches the constant state, , in the square interior as . The and solution branches coexist for some values of ( for , for ). When is large enough, the solution has lower energy than the solution.

We distinguish between the distinct solution branches by defining two measures and . In addition to the , , Constant solutions, there also exist the unstable and unstable solution branches with the same symmetries in proposition 3.2, which are indistinguishable by these measures. Hence, they appear on the same line in bifurcation diagram in figure 8 for all . The difference between the , , , and can be spotted from the associated -profiles. If on and , the corresponding solution is the solution. If on and , the corresponding solution is the solution. The and solutions also exist for . If , the solution is either the or solutions. If has isolated zero points on the square diagonals, the corresponding solution is identified to be the solution branch.

Figure 8

Bifurcation diagrams with , and from top to bottom. Left: plot of , verses . Right: plot of the energy verses . (Online version in colour.)

We numerically solve the Euler–Lagrange equations (3.3)–(3.5) with by using Newton’s method to obtain the unique stable solutions for the different values of . The initial condition is not important here, since the solution is unique and the nonlinear term is small for . We perform an increasing sweep for the , and solution branches and a decreasing sweep for the diagonal , and rotated solution branches (as described in §4c). The stable branch for is obtained by taking the stable branch, with as the initial condition. The unstable and are tracked by continuing the stable and stable branches. If the branch is given by for a fixed , then the initial condition for the unstable -solution is given by the corresponding solution, for any . The initial condition for the unstable branch is given by , where are perturbations in (4.4) and perturbation in (4.22), for any (figure 7).

Consider the case . For , there is the unique . For , the stable bifurcates into an unstable , and two stable solutions. When , the unstable bifurcates into two unstable , which are featured by isotropic lines or defect lines localized near a pair of opposite square edges. When , unstable solutions appear simultaneously. When , the and solution have the same energy. Each unstable further bifurcates into two unstable solutions. As increases, the unstable solutions gain stability. The has the highest energy amongst the numerically computed solutions for , for large . For , the ceases to exist and the unique solution is the stable solution. For , the solution is stable for and the unstable and appear for large . At the first bifurcation point , the bifurcates into two stable solutions. At the second bifurcation point, , it further bifurcates into two unstable solutions and for , the unstable and unstable solution branches appear. The and are always unstable and the solution has slightly lower energy than the . The unstable has higher energy than the unstable solutions when is large. For larger , the unique stable solution, for small , is the solution, which remains stable for . We can clearly see that the solution approaches as gets large. For , the and states disappear, and the solution does not bifurcate to any known states. The and solution branches are now disconnected from the stable solution branch. As we perform a decreasing sweep for the or solution branches, we cannot find a or solution for or , for small and . The solution has lower energy than the and solutions for large , as suggested by the estimates in §4c. For much larger values of , we only numerically observe the solution branch.

To summarize, the primary effect of the anisotropy parameter, , is on the unique stable solution for small . The elastic anisotropy destroys the cross structure of the , and also enhances the stability of the and solutions. A further interesting feature for large , is the disconnectedness of the and solution branches from the parent solution branch. This indicates novel hidden solutions for large , which may have different structural profiles to the discussed solution branches, and this will be investigated in greater detail, in future work.

In the next proposition (proof in electronic supplementary material), we prove a stability result which gives partial insight into the stabilizing effects of positive . Let be an arbitrary critical point of the energy functional (3.6). As is standard in the calculus of variations, we say that a critical point is locally stable if the associated second variation of the energy (3.6) is positive for all admissible perturbations, and is unstable if there exists an admissible perturbation for which the second variation is negative.

Proposition 5.1.

For , where  is some constant depending only on  and , the critical points of the energy functional () in the restricted admissible space are locally stable with respect to the perturbations and

Conclusion and discussions

We study the effects of elastic anisotropy on stable nematic equilibria on a square domain, with TBCs, primarily focusing on the interplay between the square edge length, , and the elastic anisotropy, . We study LdG critical points with three degrees of freedom, . We use symmetry arguments on an th of the square domain, to construct symmetric LdG critical points for which vanishes on the square diagonals, and vanishes on the coordinate axes. The is a special symmetric critical point for with . In particular, cannot be identically zero for . There are different classes of these symmetric critical points, and we perform asymptotic studies in the small and small limit, and large limits, to provide good asymptotic approximations for the novel and solutions, both of which can be stable in physically relevant regimes. The large -picture for is qualitatively similar to the case, with the stable diagonal, and rotated, solutions. The notable difference is the emergence of the competing stable solution, which is energetically preferable to the and -solutions, for large and large . This suggests that for highly anisotropic materials with large , the experimentally observable state is the solution with in the square interior. In other words, the state is almost perfectly uniaxial, with uniaxial symmetry along the -direction, and will offer highly contrasting optical properties compared to the and solutions. This offers novel prospects for multistability for highly anisotropic materials.

Another noteworthy feature is the stabilizing effect of , as discussed in §5. The solution has a central point defect in the square interior and is unstable for . However, it gains stability for moderate values of , as increases, and ceases to exist for very large positive values of . We note some similarity with recent work on ferronematics [32], where the coupling between the nematic director and an induced spontaneous magnetization stabilizes interior nematic point defects, with . It remains an open question as to whether elastic anisotropy or coupling energies (perhaps with certain symmetry and invariance properties) can stabilize interior nematic defects for tailor-made applications.