Nematic liquid crystal director structures in rectangular regions. Physical Review E, 97 (2).

We consider a shallow rectangular well of nematic liquid crystal subject to weak anchoring on the sides of the well. By considering weak anchoring instead of inﬁnitely strong anchoring, we are able to analyze nematic equilibria in the well without the need to exclude point defects at the corners, as done in previous work in the area. For relatively weak anchoring, we are able to derive analytic expressions for the director alignment angle in terms of an inﬁnite series of modes, involving roots of a transcendental equation. The analytic forms of the director conﬁguration are then used to calculate critical anchoring strengths at which uniform and distorted director structures exchange stability. We also consider the asymptotic behavior of the director structure and energy for very strong anchoring. We show that in both cases—for the transitions from uniform to distorted states and the limit of inﬁnitely strong anchoring—the approximate analytic expansions agree very well with corresponding numerical calculations of the full model.


I. INTRODUCTION
Liquid crystals are liquids in which there is some degree of orientational ordering of the constituent molecules. These molecules, often elongated organic compounds, arrange in such a way that an average orientation of one of the principle molecular axes may be defined [1]. This average molecular orientation is termed the director and is denoted mathematically by the unit vector n(x,t), an orientation that may vary in both space and time [2]. One alternative mathematical description uses both the average molecular orientation and a measure of the order about that average orientation, the general form of which leads to an analysis in terms of a second-order tensor Q, termed the orientation tensor or simply the Q tensor [3]. Liquid crystals may also exhibit positional ordering of the molecules, for instance, forming relatively distinct layered structures (smectics). However, for nematic liquid crystals, the type we consider in this article, only orientational ordering is present.
Interest in the confinement of liquid crystals between solid boundaries originally came about out of necessity since viewing a liquid crystal under a microscope was only possible if the liquid was held in place by at least one solid boundary. Later it became clear that the competition between the orientational influence of a solid boundary and those of internal effects (such as elasticity) or external effects (e.g., an applied electric field) could lead to interesting behavior of both scientific and technological interest. Indeed, in the major application area of liquid crystals, liquid crystal display (LCD) devices, the influence of the bounding surface plays a crucial role in the optical switching effects upon which the device relies [4].
The influence that a solid boundary has on the molecular alignment of a liquid crystal is the subject of a considerable body of research but is often characterized through the specification of a surface energy [5]. This surface energy, often termed an anchoring energy, is a function of a macroscopic variable, such as the director n or the Q tensor, as well as parameters that measure the strength of interaction between the molecular orientation and the boundary. The surface energy is minimized when the molecular orientation is in a preferred orientation with respect to the boundary normal, or a certain direction on the boundary that has been prescribed by a mechanical or chemical treatment. The actual molecular orientation on the boundary may differ from the preferred orientation, possibly due to the influence of elastic effects away from the boundary, leading to an increase in the surface energy. This ability of the director to move away from the preferred direction is termed weak anchoring, while a situation where the director is fixed at the preferred direction is termed infinite or strong anchoring. For weak anchoring it is the balance of the surface energy and other components of the free energy, such as the elastic energy, that lead to distorted equilibrium structures. When there are multiple boundaries in a device, each preferring different orientations, distortion cannot be avoided and it is this situation that we will consider in this article.
Confining a liquid crystal in shallow rectangular wells has been studied by a number of authors in recent years because it offers the possibility of multiple stable director configurations [6]. The relative stability of these different director structures, and the mechanisms through which the system switches between different states, is of particular interest if low-power electro-optic devices are to be developed [7]. This was the motivation behind the work of Tsakonas et al. [8], where a device consisting of an array of shallow rectangular wells was considered experimentally and theoretically. In their theoretical work, where the director was assumed to stay within a single plane, a two-dimensional analysis was effective in modeling the experimentally observed director configurations. The Q-tensor modeling of [8] imposed planar infinite anchoring of the nematic director on the boundaries of the well (where the preferred director orientation is tangential to the boundaries) and predicted multistable configurations of the director exhibiting defects at the corners of the well, all in good agreement with experimental data. Their work has subsequently led a number of researchers to consider similar geometries of confined nematics. Luo et al. [9] extend the analysis of [8] to incorporate surface energies and a degree of weak anchoring, still within the context of Q-tensor, or Landau-de Gennes, theory. They also propose a dynamic model for switching between equilibrium director states based on dielectric effects. Kusumaatmaja and Majumdar [10]a l s o model the device in [8] with a surface energy potential, computing minimum energy pathways between the stable equilibria for variable surface anchoring strength. Robinson et al. [11] compare molecular and Landau-de Gennes models for nematic equilibria in square wells. Landau-de Gennes theory is also employed by Kralj and Majumdar [12], Canevari et al. [13], and Slavinec et al. [14] in studies which allow for biaxial order reconstruction. Other work modeling confined liquid crystal systems include using Monte Carlo techniques [15,16], while Davidson and Mottram [17] derive the director orientation in a variety of geometries via conformal mappings. Studies of confined regions have not been restricted to planar (or near-planar) surface anchoring. For example, Zheng and Hu [18,19] examine models for polydimethylsiloxane microchannels where the liquid crystal molecules exhibit homeotropic ordering on boundary surfaces, so that the preferred director orientation is perpendicular to the boundaries.
In recent work on colloidal nematics in a rectangular geometry, Lewis et al. [20] examine a director model of equilibria in a well when the liquid crystal is subject to infinite planar anchoring on the four sides. However, in that situation difficulties arise because of incompatible conditions on adjacent boundaries, leading to defects at each corner for which a director-based model breaks down. Therefore, in order to calculate the energy of the system, Lewis et al. [20] remove a quarter-disk of small radius around each defect and calculate a regularized free energy over the reduced domain. An asymptotic expansion of the regularized energy is then expressed in terms of the unknown defect core radius.
In this article, we consider a shallow rectangular well of nematic liquid crystal, but instead of imposing infinite planar anchoring on the boundaries of the well, we introduce weak planar anchoring through a Rapini-Papoular [5] surface energy at each boundary. This allows us to carry out an analysis of the director configuration equilibria in the well without the need to exclude point defects at the corners of the rectangle, in other words, avoiding the problems faced in [20]. We are able to derive analytic expressions for the director alignment angle, written as an infinite series involving roots of a transcendental equation, finding a critical anchoring strength at which a uniform director configuration exchanges stability with a distorted structure. Using the analytic form for the director orientation, we are then able to examine the asymptotic behavior of our system both close to the critical anchoring strength and in the limit of infinite anchoring. The latter analysis allows a comparison with the results of previous work in this area. We also show that the asymptotic expansions agree very well with numerical calculations.

II. MATHEMATICAL MODEL
We will model a static nematic liquid crystal confined in a rectangular well of depth μd and side lengths d and λd, so that the aspect ratio in the xy plane is λ [see Fig. 1(a)]. We will assume that the depth of the well is considerably smaller than the other two lengths so that μ ≪ 1 and μ ≪ λ, and that planar degenerate anchoring on the faces z = 0 and z = μd forces the director to remain in the xy plane. These conditions lead us to assume that the director lies in the xy plane throughout the well and we may simplify the mathematical model to consider only the director configuration in the cross-sectional area specified as Ŵ ={(x, y) ∈ [0,d] × [0,λd]} [see Fig. 1(b)]. Since the nematic director is now assumed to lie in the xy plane it can be expressed in the form where θ (x, y) is the director angle measured with respect to the positive x direction. The director configuration θ (x, y) will then be determined by a minimization of the total free energy of the system, namely, the sum of the bulk elastic energy and the surface energies at the boundaries. Using the standard oneconstant approximation for the Frank elastic constants [21], we can write the elastic energy density for the nematic as where K is the elastic constant, typically of the order 10 −11 N [2], and ∇θ = (θ x ,θ y ), where subscripts indicate partial derivatives with respect to x and y.
The four boundaries of region Ŵ are denoted by σ i (i = 1to 4) as indicated in Fig. 1(b). As mentioned in the Introduction, the preferred director orientation is parallel to each boundary (i.e., planar alignment). Employing the Rapini-Papoular form for the surface energy density, we can write the surface energy density on each boundary as where ω is the constant, positive anchoring strength and ν represents the outward normal for the boundary. Although it is possible to proceed with the general situation in which the anchoring strengths on each side of the region are all different (e.g., ω i for i = 1,2,3,4), the analysis is cumbersome and little is gained in terms of general insights into this problem. Therefore, we restrict our attention here to a uniform anchoring strength ω.
Combining the integration of the elastic energy density (1) in the bulk of the well with the integration of the surface energy density (2) on each boundary, we can write the total free energy as where ds i represents integration along the surface σ i in the positive direction of the corresponding Cartesian coordinate. We will nondimensionalize the Cartesian coordinates using (x,ŷ) = (x/d ,y/d) so that the cross-sectional region is noŵ Ŵ ={(x,ŷ) ∈ [0, 1] × [0,λ]}. We now need only consider the case λ 1 since the transformation λ → 1/λ with (x,ŷ) → (ŷ,x) will provide solutions for the case λ<1.
With the nondimensionalizations mentioned above, we obtain the dimensionless total free energŷ where∇θ = (θx,θŷ) and the parameter τ = ωd/K represents a dimensionless anchoring strength, a ratio of anchoring and elastic effects or, equivalently, it can be thought of as the ratio of the well dimension d and the surface extrapolation length K/ω. Our aim is now to minimize the dimensionless energŷ W in Eq. (4) with respect to the possible director angle configurations θ (x,ŷ) for a given anchoring parameter τ and aspect ratio λ. By a standard application of the calculus of variations (see, for example, Stewart [2]), it is straightforward to show that θ (x,ŷ) must satisfy Laplace's equation in the bulk of the nematic cell, subject to a nonlinear Robin condition on each boundary, Equations (5) and (6) are solved by the trivial solutions θ (x,ŷ) ≡ 0 and θ (x,ŷ) ≡ π/2 [or the equivalent solutions θ (x,ŷ) ≡ nπ and θ (x,ŷ) ≡ (n + 1/2)π for n ∈ Z], and in the next section we will first consider solutions that bifurcate from these undistorted states as the anchoring parameter increases. We find that there are four important states, two that bifurcate from θ (x,ŷ) ≡ 0 and two from θ (x,ŷ) ≡ π/2. Importantly, the ordering of the bifurcations to the distorted states can change for different values of λ. Therefore, as anchoring increases (or, equivalently, as elasticity decreases) the first nontrivial state to appear will depend on the value of the aspect ratio λ.
We will also consider the asymptotic limit for large anchoring parameter, i.e., τ →∞, which approximates the infinite anchoring limit considered by Lewis et al. [20]. For the model with infinite anchoring at the boundaries presented in [20], there exists a discontinuity in the director angle θ at each of the four corners. Because of this discontinuity, which can be thought of as a line defect along the z direction, the free energy diverges logarithmically on approach to these points [22]. To remove this singularity, Lewis et al. [20] regularized the free energy by removing from the region a quarter-disk of radius ǫ at each corner of the well, effectively a defect core region. With infinite anchoring and this form of regularizing the problem, it is then possible to find an infinite number of equilibrium solutions. As mentioned by Lewis et al. [20], in the limit of infinite anchoring the lowest energy solutions correspond, up to symmetry, to one of three basic types of equilibrium. With our present approach we also find these three states and can link them to three of the states that bifurcate from the trivial solutions at lower values of the anchoring parameters. We also consider a fourth state found to bifurcate from a trivial solution and determine the energy of this state in the large anchoring parameter limit.
In comparison to the work of Lewis et al. [20], the model presented here in Eqs. (5) and (6) has no director discontinuities and there is no need to remove any regions to regularize the energy. Indeed, at each corner the boundary conditions (6) for two adjacent walls do not contradict each other and can both be satisfied independently. In the limit as τ →∞ in Eq. (6), we can also reach the situation modeled by [20] without their approximation and are able to determine the global energy minimizer, as well as analyze the director angle over a range of anchoring parameters and aspect ratios. While it is relatively straightforward to implement a numerical scheme to solve Eqs. (5) and (6), the nonlinear nature of our boundary conditions makes any type of analysis difficult. However, in the following sections we will show how linearization can lead to very effective results in certain limits of the anchoring strength. In the sections that follow, we omit theˆfrom all nondimensional quantities (specifically,x,ŷ,∇,Ŵ,σ i ,Ŵ ) with the understanding that, henceforth, all quantities are nondimensional.

III. BIFURCATIONS FROM UNDISTORTED STATES
As mentioned in the previous section, regardless of the anchoring parameter τ and aspect ratio λ, the system in Eqs. (5) and (6) always admits undistorted solutions θ ≡ 0 and θ ≡ π/2. Using Eq. (4), we see that these two constant equilibria correspond to total energies W 0 = 2λτ and W π/2 = 2τ , respectively, so that when λ = 1 the undistorted states have equal energies, and when λ ≶ 1weseeW 0 ≶ W π/2 . As we will examine below, distorted states occur at bifurcations from the two constant equilibria at nonzero τ , with different behavior depending on the value of λ. Each of these bifurcations is a supercritical pitchfork from the two undistorted states and we now proceed by linearizing about these states.
Linearizing the boundary conditions (6) around θ ≡ 0, we derive The solution to Laplace's equation in region Ŵ subject to boundary conditions (7) is then for any constant A 0 , and where p and τ must satisfy the simultaneous equations Equation (9) can be solved to find two solutions for τ in terms of p, which, upon substitution into Eq. (10), lead to two corresponding transcendental equations for p, A similar approach can be adopted when we linearize about θ ≡ π/2, with the resulting solution being with the solution for τ being and q satisfies the corresponding transcendental equation The solutions to the transcendental equations (12) and (15) are an infinite set of values, corresponding to mode numbers p i and q i , respectively. Each mode has associated with it a corresponding value τ given by the appropriate equation, (11) or (14). As we will see later, these values are critical anchoring strengths at which each mode appears in the system in order to reduce the free energy. Corresponding to each mode is a set of amplitudes A 0,i and A π/2,i so that the general solutions, for any value of τ , are then for the linearization about θ ≡ 0, and for the linearization about θ ≡ π/2. It is worth noting that, depending on the choice of ± solutions in Eqs. (12) and (16), In other words, one solution for θ (x, y)i nE q .( 16) will be symmetric with respect to x = 1/2, and the other antisymmetric. [Symmetry about x = 1/2 is equivalent to θ (x, y) = θ (1 − x, y) and antisymmetry is θ (x, y) =−θ (1 − x, y).] The same can be said of Eq. (17) with solutions symmetric or antisymmetric with respect to y = λ/2 through a combination of hyperbolic terms. In Fig. 2(a), we plot the smallest, positive zeros of f ± (p; λ) and g ± (q; λ) as the aspect ratio λ varies. These values have then been used to calculate the corresponding critical anchoring strengths τ in Fig. 2(b). These four lowest critical values of τ correspond to four director distortion modes, given by the appropriate solutions in Eqs. (8) and (13). Plots of these four director structures are illustrated in Fig. 3 for an aspect ratio λ = 1.5. For three of these states, D, U 1 ,U 2 ,wehaveusedthe same notation as in [20], indicating the diagonal or U-shaped nature of the distortion. We denote the fourth state by DD to recognize that it is essentially a double D state with symmetric diagonal distortions in 0 <y<λ/2 and λ/2 <y<λ.Weare able to classify the four different branches, and associate them with solutions in [20], because the states exhibit particular (anti)symmetries with respect to x = 1/2 and/or y = λ/2. For example, for the D state, θ (x, y) = θ (1 − x, λ − y). Figure 2(b) indicates that for the two states which bifurcate from the θ ≡ π/2 trivial solution (the D and U 2 states), it is always the D state which bifurcates at the lower critical τ value. In fact, it is simple to show analytically that the first nonzero solution to g + (q; λ) = 0 tends to q = π/2 from above as λ →∞, and the first nonzero solution to g − (q; λ) = 0 tends to q = π/2 from below as λ →∞. The asymptotic behavior for the corresponding critical values of the anchoring parameter is then τ → π/2 + for the U 2 state and τ → π/2 − for the D state. For the two states bifurcating from the θ ≡ 0 trivial solution (the DD and U 1 states) the situation is slightly more complicated. While the solutions of f ± (p; λ) = 0 both tend to p = 0 from above as λ →∞, the critical value of the anchoring parameter behaves as τ → 2 + for the U 1 state and τ → 0 + for the DD state. The values of the critical anchoring parameters at λ = 1 are found from Eq. (12)tobeτ ≈ 2.55 for the U 1 state and τ ≈ 4.61 for the DD state. There is, therefore, a value of λ at which the critical values of τ for the U 1 and DD states cross, as can be seen in Fig. 2(b). This value is found numerically to be λ c ≈ 1.75, at which τ c ≈ 2.24.
The order of the bifurcations is illustrated in Fig. 4, where we have considered equilibrium states obtained both numerically from the full nonlinear problem in Eqs. (5) and (6)  . We see that, for λ<λ c ,t h eU 1 state bifurcates at a lower value of τ than the DD state while, for λ>λ c , the ordering exchanges and the DD state bifurcates at the lower value of τ . In fact, for much larger values of λ,w e would expect a DDD state (or even states we might term D n which are similar to n repeated D states) to emerge and eventually to be the state that bifurcates from the trivial θ ≡ 0 state at the lowest value of τ . For bifurcations from the trivial state θ ≡ π/2, however, the D state always bifurcates at the lower value of anchoring parameter τ . As would be expected with a perturbation method, the analytical energy calculation is only in good agreement with the full numerical solution close to the respective bifurcation points. It is worth noting that relaxation of the one-constant approximation for the Frank elastic constants may lead to a change in the behavior of the system since the relative energy cost of splay and bend distortions would be altered. With differing amounts of splay and bend distortions in each of the states, we might expect the values and order of the bifurcations from the trivial states to be changed. However, previous work in similar systems [24] suggests that elastic anisotropy may not significantly affect the stability of states so that the qualitative behavior would remain the same.

IV. LARGE ANCHORING STRENGTH ANALYSIS
We now consider the situation of a large anchoring strength τ corresponding to, for example, a well dimension d that is much larger than the surface extrapolation length K/ω.I n this case, the director is anchored relatively strongly at each well boundary and we may assume that on σ i ,awayfromthe corners, the director angle θ is close to the constant angle the director takes in the infinite anchoring limit. We will denote these infinite anchoring limiting angles at each boundary by i . This allows us to linearize each of the boundary conditions (6) about θ = i , replacing them with The four bifurcation states considered in the previous section are characterized by the preferred boundary orientations with the DD state having a change in orientation along the walls at x = 0 and x = 1sothat We can associate these linearized boundary conditions with quadratic forms of the surface energy densities in Eq. (2) via (ν · n) 2 ≈ (θ − i ) 2 (i = 1,2,3,4), up to an additive constant that will play no role in the minimization of the total energy. As mentioned above, this linear approximation of the nonlinear boundary conditions will be valid everywhere except at the corners of the region and also the points (0,λ/2), (1,λ/2) for the DD state. However, as we will show later, up to first order, when τ is large, our analysis suggests that the differences between the energy associated with solutions of the linear problem and the nonlinear problem are parameter-independent constants. Therefore, after these parameter-independent constants are determined [for each of the states (19)- (22)], we obtain the τ →∞asymptotic expressions for the energies of each state.
The solution of Laplace's equation in region Ŵ subject to boundary conditions (18) is found by separation of variables. For the D, U 1 , and U 2 states, since the system is now linear, one need only derive the solution in the case 1 = 0, 2 = 3 = 4 = 0, and then employ the principle of superposition together with appropriate rescaling and rotation, as is done in [20]. The complication in our analysis is the presence of Robin boundary conditions which lead to eigenvalues that are solutions of a transcendental equation [25].Detailsofthe calculations for all four distorted states are given in Appendix A [including an explanation for why we need consider only odd j in Eq. (24) and subsequent analysis]. The solution for the D, U 1 ,orU 2 states is most compactly written as where with the T dependence entering through the eigenvalues P j (T ), solutions of the transcendental equation T − P j tan(P j /2) = 0,P j ∈ (jπ − π/2,j π )( j = 1,3,5,...).
We see immediately from Eq. (25) that j is symmetric with respect to U = 1/2 through a single trigonometric term, cos[P j (U − 1/2)]. This will lead to symmetry in the x or y direction about the center of the rectangle for each particular state, depending on the combination of terms in solution (24) and the appropriate i (i = 1,2,3,4). It is less obvious from the nature of the V -dependent term in j , but when combined with the different choices of i in Eq. (24), symmetry or antisymmetry is also introduced for the other xy coordinate through the addition of the hyperbolic terms. The special form of the series solution θ (x, y)f o rt h eDD state is given in Appendix A. Figure 5 shows the director configuration for the series solutions θ (x, y) for the four different sets of preferred directions i given by Eqs. Having found the equilibrium solution (24), we can calculate the total free energy W in Eq. (4) associated with the director structure, albeit using the quadratic forms of the surface energy densities (23). Although W is now quadratic in θ (x, y), through both the elastic and surface energy terms, it is possible to simplify the expression using Green's first identity [26] and the boundary condition (18): The free energy can now be calculated using the series solution (24) integrated along the four boundaries.
In order to derive a compact expression for the W D energy, we first consider the E j term in Eq. (27). Given that we are considering large anchoring strengths, we can simplify and solve the transcendental equation (28) to obtain This allows us to approximate E j as It is possible to express the sum of the approximation for E j over odd j 1 in terms of the digamma function (z); the full expression is given in Appendix B and we will denote it by E sum (λ,τ). Hence we may now write W D as Asymptotic expansions for E j F j and E sum are then possible as τ →∞, leading to where γ ≈ 0.57721 is the Euler-Mascheroni constant [27] and we have introduced the functions With analysis similar to that above, we can also derive the asymptotic expansion of the energies of the U 1 and U 2 states. The analysis for the DD state is more complicated, but we recognize that, in the limit of infinite anchoring, the DD state in a confined region of aspect ratio λ is effectively the same as two adjacent D states of aspect ratio λ/2. The approximate energies for the four states are then Following Bruckman [28], we could write the two infinite sums s 1 (·) and s 2 (·) in the alternative forms where m is related to λ by λ = K(1 − m)/K(m) and K is the complete elliptic integral of the first kind. The energies (32)-(35) now reduce to The differences between the approximate energies in Eqs. (32)-(35) and an asymptotic result for the full nonlinear problem using Rapini-Papoular boundary conditions (6)a r e the errors due to the linear approximation of the boundary conditions close to the corners of the region. However, comparison of Eqs. (37)-(40) with numerical calculations for the nonlinear system (discussed in the following section), restricting attention to large values of τ , indicates that this error is neither a function of λ nor τ . This suggests that, in the leading order, O(ln(τ )), and first order, O(1), terms, these errors need only be obtained numerically once, and are the same for all instances of a nematic confined in a rectangular region being independent of any geometric or material properties. It should be noted that the parameter independence of these constants, ǫ DU and ǫ DD in the expressions below, is numerically obtained and we have not proved this result analytically. However, after extensive calculations for a wide range of physically relevant parameters we can have a high level of confidence in this assertion. The asymptotic results for the four full nonlinear energies are then where ǫ DU =−0.06 and ǫ DD = 0.22. From Eqs. (41)-(44), we see that the asymptotic expressions for W U 1 and W U 2 coincide for the special case of a square domain, i.e., λ = 1. (This is to be expected from the symmetry of a square nematic well.) When λ is small, for which m is close to 1 in Eqs. (37)-(39), it follows from Eqs. (37) and (38) that the asymptotic behavior of W D is very similar to that of W U 1 . The same can also be said about W D and W U 2 when λ is large, corresponding to m close to zero in Eqs. (37)-(39). Given that W DD (λ) ≈ 2W DD (λ/2), we can find similar relationships between W DD and W U 1 or W U 2 , except that the constant error ǫ DD = ǫ DU due to the presence of high distortion regions at the midpoints of the side walls x = 0 and x = 1. The asymptotic expansions for the U 1 and U 2 states differ from the D state through the logarithmic terms in Eqs. (38) and (39), respectively. Therefore, since m ∈ (0, 1), it follows that, in the limit τ →∞,W D is the state with lowest energy. As mentioned in the previous section, anisotropy in the elastic constants could influence the nature of the minimum energy state.
Comparing Eqs. (37)-(40) to equivalent expressions in Lewis et al. [20], we see that we have agreement in a number of terms [for example, those involving ln(τ ), ln(2λ/π) and the s 1 (·),s 2 (·) functions], but we have additional terms at O(1) which are missing from [20]. Presumably this is due to the need in [20] to remove parts of the region (at the corners) to produce an analytically tractable problem.

V. COMPARISON OF ANALYTICAL AND NUMERICAL ENERGIES
We can now compare various aspects of the approximate asymptotic expansions found in Eqs. (41)-(44) to the energies derived from the numerical solution of the full nonlinear problem in Eqs. (5) and (6). Unlike the previous section, here we consider a range of anchoring strengths, including those close to the critical bifurcation values derived in the analysis of very weak anchoring. We solve the full nonlinear system using COMSOL with a nonuniform mesh, refined at the corners and sides of the region, ensuring that the mesh is fine enough so that further refinement does not alter the energy calculation by more than 1%.
In Fig. 6 we plot the four asymptotic energy expansions from Eqs. (41)-(44), up to O(ln(τ )/τ ), against λ for τ = 10 and τ = 100 (dashed lines). In addition, we also plot the numerical energies for the nonlinear system (solid lines). As expected, there is a significant difference between the asymptotic forms and numerical results for a relatively small value τ = 10, whereas for a large value τ = 100, the leading terms in the asymptotic results show good accuracy over a range of aspect ratios λ. them to the energies obtained numerically (solid lines) for λ = 1.5. The solutions in Eqs. (16) and (17) provide approximate energies for the weakly anchored system close to bifurcations from the trivial states (small τ , dash-dot lines), while the asymptotic behavior as τ →∞is given by the leading terms in Eqs. (41)-(44)(largeτ , dashed lines). The energies W 0 and W π/2 are also indicated. Note that it is difficult to distinguish the graphs of U 2 and D states in Fig. 7(b), as was suggested by the form of the asymptotic expansions in Eqs. (41), (43), and mentioned earlier. The analytical approximate energies agree very well with the equivalent numerical graphs at anchoring strengths that are close to the bifurcation from the trivial states and also at large values of τ , with less than 1% discrepancy for τ 10 2 .

VI. CONCLUSIONS
We have found approximate analytical solutions for the nematic director angle configuration in a weakly anchored rectangular region. The use of a standard weak anchoring energy (the Rapini-Papoular surface energy) means we are able to find solutions and calculate energies, without the need to extract the core of defects to allow the system to be tractable analytically. In two important limits, for anchoring strengths very close to the bifurcation from a trivial state (i.e., close to the point at which the distorted state comes into existence), and in the high anchoring strength or weak elasticity limit, we find good agreement with a numerical solution of the full nonlinear problem.
Given typical values of the Frank elastic constants K ≈ 10 −11 N and anchoring strength ω ≈ 10 −4 Nm −1 [2], the high-τ energy expressions in Eqs. (41)-(44) will be good approximations for wells of side length d 10 μm, so that τ = ωd/K 100. For larger anchoring strengths of ω ≈ 10 −3 Nm −1 , the energy expressions are accurate for a wider range of well dimensions, with d 1 μm. Since the accuracy of most common forms of construction of such wells (i.e., photolithography) is around the length scale of microns, it is clear that the asymptotic energies are most likely to be valid for all but the weakest of anchoring strengths. However, in this high-τ limit it will always be the D state that is the global energy minimizer. Therefore, if bistability is required, with the possibility of switching between stable states, it may be useful to consider anchoring strengths closer to those which occur at the critical anchoring parameters for bifurcation from trivial states. It is at these anchoring strengths that it will be easiest to switch between states since the energy barriers are smaller.

ACKNOWLEDGMENTS
This work was supported by the United Kingdom Engineering and Physical Sciences Research Council (EPSRC) via Grant No. EP/M508159/1.

APPENDIX A
A standard separation of variables approach to find the solution of Laplace's equation in region Ŵ subject to Robin boundary conditions (18) leads to series solutions for the D, U 1 , and U 2 states in the form of Eq. (24), only with the sums taken over all j ∈ Z + rather than restricted to odd values. In their most general forms, j (U,V, ,T ) = M j × N j (j = 1,2,3,...), where M j (U,T ) = √ 2[cos(P j U )P j + T sin(P j U )] P j 2 + T 2 + 2T are orthonormal with respect to U ∈ (0, 1) and N j (V, ,T) = √ 2T 2 P j 2 + T 2 cos(P j ) + P j 2 − T 2 [cosh(P j V )P j + T sinh(P j V )] P j P 2 j − T 2 P j 2 + T 2 + 2T P j 2 + T 2 sinh(P j ) + 2 P j T cosh(P j ) .
However, upon substitution of T we subsequently find N j (V, ,T) = 0forj even. Therefore, the only contributions to θ (x, y) come from eigenvalues P j (T ) lying in the second quadrant satisfying the transcendental equation (26). Furthermore, if we also replace T in j with P j tan(P j /2) (odd j 1), we eventually obtain the simplified form given by (25). This expression is then used to construct the solution θ (x, y)forthe D, U 1 , and U 2 states. The DD solution can also be found by separation of variables, though the derivation is slightly different due to the piecewise nature of boundary condition (22). Following some analysis, we can show that the director angle θ (x, y)f o rt h e DD state can be expressed as the series where the eigenvalues Q j are the positive solutions of the transcendental equation λτ tan(Q j ) + 2Q j = 0.
One obvious feature of this director profile is its symmetry in the horizontal direction with respect to x = 1/2 and antisymmetry about the center in the vertical direction.