Variational and numerical analysis of a $\mathbf{Q}$-tensor model for smectic-A liquid crystals

We analyse an energy minimisation problem recently proposed for modelling smectic-A liquid crystals. The optimality conditions give a coupled nonlinear system of partial differential equations, with a second-order equation for the tensor-valued nematic order parameter $\mathbf{Q}$ and a fourth-order equation for the scalar-valued smectic density variation $u$. Our two main results are a proof of the existence of solutions to the minimisation problem, and the derivation of a priori error estimates for its discretisation of the decoupled case (i.e., $q=0$) using the $\mathcal{C}^0$ interior penalty method. More specifically, optimal rates in the $H^1$ and $L^2$ norms are obtained for $\mathbf{Q}$, while optimal rates in a mesh-dependent norm and $L^2$ norm are obtained for $u$. Numerical experiments confirm the rates of convergence.


INTRODUCTION
Smectic liquid crystals (LC) are layered mesophases that have a periodic modulation of the mass density along one spatial direction. Informally, they can be thought of as one-dimensional solids along the direction of periodicity and two-dimensional fluids along the other two remaining directions. According to different in-layer structures, several types of smectic LC are recognised, e.g., smectic-A, smectic-B and smectic-C (see [26,Figure 9] for illustrations of different smectic phases). In particular, in the smectic-A phase, the molecules tend to align parallel to the normals of the smectic layers. For a broader review of liquid crystals, see [2,16,28].
Several models have been proposed to describe smectic-A liquid crystals. The classical de Gennes model [15] employs a complex-valued order parameter to describe the magnitude and phase of the density variation. This complex-valued parameterisation leads to some key modelling difficulties, which motivated the development by Pevnyi, Selinger & Sluckin (PSS) of a smectic-A model with a real-valued smectic order parameter and a vector-valued nematic order parameter [24]. The use of a vector-valued nematic order parameter limits the kinds of defects the model can permit [2], and hence Ball & Bedford (BB) proposed a version of the PSS model employing a tensor-valued nematic order parameter [3] instead. The model considered in this work is similar to the BB model, but with additional modifications to render it amenable to numerical simulation (see [31] for details). The model was used to numerically simulate several key smectic defect structures, such as oily streaks and focal conic domains, for the first time.
Keywords and phrases: C 0 interior penalty method, a priori error estimates, finite element methods, smectic liquid crystals * The work of JX is supported by the National Natural Science Foundation of China (No. 12201636)  While the numerical modelling of nematic liquid crystals is now mature, relatively little work has considered smectics. Garcìa-Cervera and Joo [19] perform numerical simulations of the classical de Gennes model [15] in the presence of a magnetic field, using a combined Fourier-finite difference approach. Wittmann et al. use density functional theory to investigate the topology of smectic liquid crystals in complex confinement [30]. Monderkamp et al. examine the topology of defects in two-dimensional confined smectics with the help of extensive Monte Carlo simulations [22]. Our goal in this work is to analyse a model for smectic-A liquid crystals, and its finite element discretization, that was recently proposed by Xia, MacLachlan, Atherton and Farrell in [31].
We consider an open, bounded and convex spatial domain Ω ⊂ R d , d ∈ {2, 3} with Lipschitz boundary ∂ Ω. The smectic-A free energy we analyse in this work is given by where I d denotes the d × d identity matrix, D 2 is the Hessian operator, f s (u) is the smectic bulk energy density (2) Here, B > 0 is the nematic-smectic coupling parameter, a 1 , a 2 , a 3 represent smectic bulk constants with a 3 > 0, K > 0 and l > 0 are nematic elastic and bulk constants, respectively, q ≥ 0 is the wave number of the smectic layers, and the trace of a matrix is given by tr(·). The two main contributions of this paper are to prove existence of minimisers to the problem of minimising J over a suitably-defined admissible set, and to derive a priori error estimates for its discretisation using the C 0 interior penalty method [7,9]. We show the existence result of minimising the free energy eq. (1) in section 1. We then derive a priori error estimates for both Q and u in the simplified case q = 0 in section 2. We confirm that the theoretical predictions match numerical experiments in section 3, for both q = 0 and q > 0.
In order to avoid the proliferation of constants throughout this work, we use the notation a b (respectively, b a) to represent the relation a ≤ Cb (respectively, b ≥ Ca) for some generic constant C (possibly not the same constant on each appearance) independent of the mesh.

EXISTENCE OF MINIMISERS
To formulate the minimisation problem for eq. (1), we must first define the admissible space A in which minimisers will be sought. We define A as with the specified Dirichlet boundary data Q b ∈ H 1/2 (∂ Ω, S 0 ) and u b ∈ H 3/2 (∂ Ω, R). Here, S 0 denotes the space of symmetric and traceless d × d real-valued matrices. For simplicity of the analysis, we only consider Dirichlet boundary conditions for Q and u here, but other types of boundary conditions (e.g., a mixture of the Dirichlet and natural boundary conditions for u as illustrated in the implementations in [31]) can be taken. With this admissible space, we consider the minimisation problem min Notice that f n (Q, ∇Q) is the classical Landau de Gennes (LdG) free energy density [16,23] for nematic LC and it is proven by Davis & Gartland [14,Corollary 4.4] that there exists a minimiser of the functional Ω f n (Q, ∇Q) for Q ∈ H 1 (Ω, S 0 ) in three dimensions. Moreover, Bedford [4,Theorem 5.18] gives an existence result for the BB model in three dimensions: with an admissible space where SBV denotes special functions of bounded variation. For simplicity, we have ignored the surface integral here in the energy functional of the BB model. One can observe its resemblance to eq. (1). Motivated by the above results, we prove the existence of minimisers to eq. (1) via the direct method of calculus of variations. Proof. Note that both the smectic density f s and the nematic bulk density f b n are bounded from below as a 3 , l > 0. Thus, J is also bounded from below and we can choose a minimising sequence {(u j , Q j )}, i.e., Here we setQ ∈ H 1 (Ω, S 0 ) (resp.ũ ∈ H 2 (Ω, R)) to be any function with trace Q b (resp. u b ). We tackle the three terms in J (u, Q) separately in the following. First, for the nematic energy term Ω f n (Q, ∇Q), we can follow the proof of [14,Theorem 4.3] to obtain that f n (Q j , ∇Q j ) is coercive in H 1 (Ω, S 0 ), i.e., f n grows unbounded as Q j 1 → ∞, and thus the minimising sequence {Q j } must be bounded in H 1 (Ω, S 0 ). Since H 1 (Ω) is a reflexive Banach Space, we have a subsequence (also denoted as {Q j }) that weakly converges to Q ⋆ ∈ H 1 (Ω, S 0 ) such that Q ⋆ −Q ∈ H 1 0 (Ω, S 0 ), and from the Rellich-Kondrachov theorem it follows that Q j → Q ⋆ in L 2 (Ω) and ∇Q j ⇀ ∇Q ⋆ in L 2 (Ω). The weak lower semi-continuity of the nematic energy density f n in eq. (2) is guaranteed by [14,Lemma 4.2], therefore, For the smectic bulk term Ω f s (u), we can follow the proof in [4,Theorem 5.19]. By eq. (5), we have which implies an upper bound for ∇u j using [4,Ineq. (5.42)]. Hence, {u j } is bounded in H 2 (Ω) and thus there is a subsequence (also denoted as {u j }) such that u j ⇀ u ⋆ in H 2 (Ω) and u ⋆ −ũ ∈ H 2 ∩ H 1 0 (Ω). Moreover, one can readily check that u ⋆ ∞ < ∞ by the Sobolev embedding of H 2 (Ω) into the Hölder spaces C t,κ 0 (Ω) (t + κ 0 = 1 for d = 2 and t + κ 0 = 1/2 for d = 3) and the boundedness of domain Ω. Again, it follows from the Rellich-Kondrachov theorem that u j → u ⋆ in L 2 (Ω) and D 2 u j ⇀ D 2 u ⋆ in L 2 (Ω). Noting f s (u) is bounded from below for all u ∈ H 2 (Ω), we then obtain Now, we consider the nematic-smectic coupling term Ω B D 2 u + q 2 Q + I d d u 2 in J (u, Q). Note that when the wave number q = 0, this term reduces to Ω B|D 2 u| 2 and it is straightforward to obtain the weak lower semi-continuity property. Therefore, we discuss the case of q > 0 in detail as follows. By the H 1 -boundedness property of the minimising sequence {Q j } and the fact that u ⋆ ∞ < ∞, we can deduce Hence, u j Q j → u ⋆ Q ⋆ in L 2 (Ω), and further, Hence, we can conclude that J (u ⋆ , Q ⋆ ) achieves its minimum in the admissible space A by combining eq. (6), eq. (7) and eq. (8).
Remark 1.2. We briefly compare the proof with that of Bedford [4,Theorem 5.18]. In that work, the admissible space A BB included a uniaxial constraint.This assumption is necessary to guarantee the H 1 -boundedness property of the minimising sequence {Q j }, since that work seeks Q ∈ SBV (Ω, S 0 ) instead of H 1 . Enforcing this constraint numerically is difficult [6]; in this work we prefer the choice Q ∈ H 1 , which enables us to remove the uniaxiality assumption.

A PRIORI ERROR ESTIMATES
We now consider the discretisation of the minimisation problem eq. (4). For simplicity, we analyse the decoupled case with q = 0, where eq. (4) splits into two independent problems: one for the smectic density variation u: and one for the tensor field Q, Here, One can derive the following strong forms of their equilibrium equations using integration by parts. The optimality condition for u yields a fourth-order partial differential equation (PDE) where ν denotes the outward unit normal. Note that the second boundary condition of u is a natural boundary condition on the second derivative of u. For simplicity, we only consider the case of a cubic nonlinearity (i.e., a 2 = 0) here as the quadratic term can be tackled similarly. Therefore, we consider the following governing equations Meanwhile, the optimality condition for Q yields a second-order PDE We now consider these two problems (P1) and (P2) in turn.
Remark 2.1. The uniqueness of solutions is not expected. It is well-known that eq. (10) can support multiple solutions [25], while numerical experiments in [31] indicate the existence of multiple solutions to the optimality conditions for eq. (4).

A priori error estimates for (P1)
Since the PDE eq. (9) for the density variation u is a fourth-order problem, a conforming discretisation requires a finite dimensional subspace of the Sobolev space H 2 (Ω), which necessitates the use of C 1 -continuous elements. The construction of these elements is quite involved, particularly in three dimensions; without a special mesh structure, the lowest-degree conforming elements are the Argyris [1] and Zhang [32] elements, of degree 5 and 9 in two and three dimensions respectively. One approach to avoid such complexity is to use mixed formulations by solving two secondorder systems, and we refer to [12,27] for instance. However, this substantially increases the size of the linear systems to be solved. Alternatively, one can directly tackle the fourth-order problem with non-conforming elements, that do not satisfy the C 1 -requirement. For instance, the so-called continuous/discontinuous Galerkin methods and C 0 interior penalty methods (C 0 -IP) are analysed in [9,17], combining concepts from the theory of continuous and discontinuous Galerkin methods. Essentially, these methods use C 0 -conforming elements and penalise inter-element jumps in first derivatives to weakly enforce C 1 -continuity. This has the advantages of both convenience and efficiency: the weak form is simple, with only minor modifications from a conforming method, and fewer degrees of freedom are used than with a fully discontinuous Galerkin method.
We thus adopt the idea of C 0 -IP methods to solve the nonlinear fourth-order problem (P1) and derive a priori error estimates regarding u. We adapt the techniques of [21] to prove convergence rates with the use of familiar continuous Lagrange elements for the problem (P1). The weak form of eq. (9) is defined as: where for t, w ∈ H 2 (Ω), and for µ, ζ , η, ξ ∈ H 2 (Ω), Since eq. (11) is nonlinear, we derive its linearisation: where ·, · H 2 represents the dual pairing between H 2 (Ω) ∩ H 1 0 (Ω) ⋆ and H 2 (Ω) ∩ H 1 0 (Ω). It is straightforward to derive the coercivity and boundedness of the bilinear operator A s (·, ·) with the semi-norm | · | 2 (in fact, this is indeed a norm in H 2 (Ω) ∩ H 1 0 (Ω)).
Let T h be a mesh of Ω with T denoting an element, and let E I (resp. E B ) denote the set of all interior (resp. boundary) edges/faces e of the mesh, and E := E I ∪ E B . Define the broken Sobolev space We take a nonconforming but continuous approximation u h for the solution u of eq. (11), that is to say, with Q deg denoting piecewise continuous polynomials of degree deg on a mesh of quadrilaterals or hexahedra. Following the derivation of the C 0 -IP formulation given in [7, Section 3], we introduce the discrete nonlinear weak form: find where for all u h ,t h ∈ W h , Here, ε is the penalty parameter (to be specified in section 3), h e indicates the size of the edge/face e and the average of the second derivatives of u along the normal direction across the edge/facet e is defined as For any interior edge e ∈ E I shared by two cells T − and T + , we define the jump v by v = v − · ν − + v + · ν + with ν − , ν + representing the restriction of outward normals in T − , T + respectively. On the boundary edge/face e ∈ E B , we define v = v · ν. The operator P s h penalises the first derivatives across the interior edge/facet since the function in H 1 (Ω) is not necessarily continuously differentiable.
The linearised version is to seek v h ∈ W h,0 such that where We also define the mesh-dependent H 2 -like semi-norm for v ∈ W h , Note that |||·||| h is indeed a norm on W h,0 . This norm will be used in the well-posedness and convergence analysis below. We first give an immediate result about the consistency of the discrete form eq. (13).
. The solution u of the continuous weak form eq. (11) solves the discrete weak problem eq. (13).
Proof. Multiplying the fourth-order term 2B∇ · (∇ · (D 2 u)) in eq. (9) with a test function t ∈ W h,0 and using piecewise integration by parts with the boundary condition specified in eq. (9) for u, one can obtain Since u ∈ H 4 (Ω) implies ∇u is continuous on the whole domain Ω, the jump term ∇u then becomes zero and we can thus symmetrise and penalise the form eq. (18). This leads to the presence of A s h (u,t) + P s h (u,t). The remaining terms involving B s and C s are straightforward as one takes the test function t ∈ W h,0 . Therefore, u satisfies eq. (13).
By noting that D 2 : it is natural to extend the classical elliptic regularity result [5] for the biharmonic operator ∆ 2 to the case of the bi-Hessian operator D 2 : D 2 . If the domain Ω is smooth, the weak solutions of the biharmonic problem (e.g., [7, Example 2]) belong to H 4 (Ω) by classical elliptic regularity results and thus we make this assumption henceforth.
Moreover, to facilitate the analysis, we further assume that u is an isolated solution, i.e., there is only one solution u satisfying eq. (9) within a sufficiently small ball {v ∈ H 2 (Ω) ∩ H 1 0 (Ω) : |v − u| 2 ≤ r b } with radius r b . These assumptions then imply that the linearised operator DN s (u)·, · H 2 satisfies the following inf-sup condition: as the edge/facet size h e < 1. With the estimate eq. (20) at hand, we then apply the Cauchy-Schwarz inequality and use the definition eq. (17) of |||·||| h to obtain the boundedness of A s h (·, ·) and P s h (·, ·). We omit the details of the proofs here and only illustrate the boundedness result for B s (·, ·, ·, ·) and C s (·, ·) below. Lemma 2.4. (Boundedness of B s (·, ·, ·, ·) and C s (·, ·).) For u, v, w, p ∈ W h,0 , we have Proof. By Hölder's inequality, the Sobolev embedding H 1 (Ω) ֒→ L 4 (Ω), and the fact that the It then follows from a Poincaré inequality [11, Eq. (5.7)] for piecewise H 2 functions that The boundedness of C s (·, ·) follows similarly by the Cauchy-Schwarz inequality, the Sobolev embedding H 1 (Ω) ֒→ L 2 (Ω) and the use of eq. (23). The proof of eq. (22) is analogous to that of eq. (21) with a use of the embedding result H 2 (Ω) ֒→ L ∞ (Ω) and the Cauchy-Schwarz inequality.

We give the coercivity result for the bilinear form
Proof. By eq. (20) and the inequality of geometric and arithmetic means, we deduce for v ∈ W h , provided the penalty parameter ε is sufficiently large with the generic constant C from eq. (20).
An important question about the well-posedness is the coercivity of the bilinear operator DN s h (u h )·, · . Due to the presence of B s and C s terms in DN s h (u h )·, · , it is not trivial to derive its coercivity. We first discuss the weak coercivity of the bilinear form DN s h (u)·, · defined as To this end, we will employ the enrichment operator E h : W h → W C ⊂ H 2 (Ω) with W C being the Hsieh-Clough-Tocher macro finite element space [7]. The following lemma is adapted to our notation and definition of |||·||| h from [7, Lemma 1].
With this, we obtain the following discrete inf-sup condition for the discrete bilinear operator DN s h (u)·, · . Theorem 2.7. (Weak coercivity of DN s h (u)·, · ) Let u be a regular isolated solution of the nonlinear continuous weak form eq. (13). For a sufficiently large ε and a sufficiently small mesh size h, the following discrete inf-sup condition holds with a positive constant β c > 0: Proof.
For v ∈ H 2 (Ω) ∩ H 1 0 (Ω), it follows from the boundedness result of B s ,C s that B s (u, u, v, ·) and C s (v, ·) ∈ (L 2 (Ω)) ⋆ . Furthermore, since A s (·, ·) is bounded and coercive as given by lemma 2.2, for a given v h ∈ W h with |||v h ||| h = 1, there exists ξ and η ∈ H 4 (Ω) ∩ H 1 0 (Ω) that solve the linear systems: (27b) We require the H 4 -regularity for the derivation of eq. (34) below. It then follows from the standard elliptic regularity result that η 4 C BC with a constant C BC depending on u 2 .
Subtracting eq. (27a) from eq. (27b), then taking w = η − ξ and using the coercivity of A s (·, ·), boundedness of B s (·, ·, ·, ·), C s (·, ·) and enrichment estimate lemma 2.6, we get Since u is a regular isolated solution of eq. (11), it yields by eq. (19), eq. (27b), lemma 2.2, the fact that ∇(E h v h + η) = 0 and the triangle inequality, that there exists Note that ∇ξ = 0 on E I since ξ ∈ H 4 (Ω). We can thus calculate using lemma 2.6 that Further, by the triangle inequality, we get Since v h + I h ξ ∈ W h , it follows from the coercivity result eq. (24) that there exists w h ∈ W h with |||w||| h = 1 such that where in the last equality we have used the fact that because of eq. (27a) and ∇ξ = ∇E h w h = 0.
Using the boundedness result lemma 2.4 and the enrichment estimate lemma 2.6, we obtain By the boundedness of the bilinear form A s h + P s h and standard interpolation estimate, we have Moreover, by the fact that ∇ξ = ∇(E h w h ) = 0, the enrichment estimate lemma 2.6 and eq. (18), there holds Combine eqs. (32) to (34) in eq. (31) to obtain Substituting eq. (35) into eq. (30) and using standard interpolation estimates yield that A use of eqs. (35) and (36), standard interpolation estimates and eq. (28) in eq. (29) leads to Then, by the triangle inequality, we have Therefore, for the mesh size h satisfying h 2 + h min{deg −1,2} < 1 2C t , the discrete inf-sup condition eq. (26) holds for β c = 1 2C t . We can now obtain the discrete inf-sup condition for the perturbed bilinear form DN s h (I h u)·, · given by Here, we employ the interpolation operator I h : Let u be a regular isolated solution of the nonlinear continuous weak form eq. (13) and I h u the interpolation of u. For a sufficiently large ε and a sufficiently small mesh size h, the following discrete inf-sup condition holds: It follows from the definition of B s and its boundedness result lemma 2.4 that where C 1 is the generic constant arising in the boundedness result lemma 2.4 for B s (·, ·, ·, ·). Therefore, we obtain that Now using the inf-sup condition theorem 2.7 for the bilinear form DN s h (u)·, · , boundedness result lemma 2.4 and interpolation estimates, we get for a sufficiently small mesh size h such that h min{deg −1,k u −2} < β c 2C 2 . Here, C 2 depends on C 1 and u k u and k u ≥ 4 gives the regularity of u, i.e., u ∈ H k u (Ω). Therefore, the inf-sup condition eq. (38) holds.

Convergence analysis
We proceed to the error analysis for the discrete nonlinear problem eq. (13).
for v h , w h ∈ W h,0 . Due to the weak coercivity property in theorem 2.8, the nonlinear map µ h is well-defined. The existence and local uniqueness of the solution u h to the discrete nonlinear problem eq. (13) will be proven via an application of Brouwer's fixed point theorem, which necessitates the use of two auxiliary lemmas illustrating that (i) µ h maps from a ball to itself; and (ii) the map µ h is contracting. Lemma 2.9. (Mapping from a ball to itself) Let u be a regular isolated solution of the continuous nonlinear weak problem eq. (11). For a sufficiently large ε and a sufficiently small mesh size h, there exists a positive constant R(h) > 0 such that: Proof. Note that the solution u ∈ H 2 (Ω) ∩ H 1 0 (Ω) of eq. (11) satisfies the discrete weak formulation eq. (13) due to the consistency result theorem 2.3, that is to say, there holds that By the linearity of DN s h (I h u)·, · H 2 , the definition eq. (39) of the nonlinear map µ h and eq. (40), we calculate We rearrange terms in N 3 and use the boundedness result lemma 2.4 and the interpolation result [9, Eq. (5. 3)] to obtain We use the definition of B s (·, ·, ·, ·) and its boundedness to deduce that By the inf-sup condition eq. (38) for the perturbed bilinear form, we further deduce that there exists a w h ∈ W h with Note that there are other cases when k u > 4 and we only focus on the case of k u ≤ 4 here for brevity. The idea of the remainder of the proof is to choose an appropriate R(h) so that |||I h u − µ h (v h )||| h ≤ R(h). For simplicity of the calculation, we illustrate the case when 2 ≤ deg ≤ 3, k u ≤ 4. To this end, we take R(h) = 4C u h min{deg−1,k u −2} and choose h satisfying This completes the proof.

Lemma 2.10. (Contraction result) For a sufficiently large ε, a sufficiently small mesh size h and any v
Proof. For w h ∈ W h , we use the definition eq. (39) of the nonlinear map µ h , the definition eq. (37) and linearity of DN s h (I h u)·, · to calculate We make some elementary manipulations and use the boundedness of B s and the inequality of geometric and arithmetic means to yield By the inf-sup condition eq. (38), we know that there exists w h ∈ W h with |||w h ||| h = 1 such that The existence and local uniqueness of the discrete solution u h can now be obtained via the application of Brouwer's fixed point theorem [20].
It follows from theorem 2.11 that optimal convergence rates are achieved in the mesh-dependent norm |||·||| h . This will be numerically verified in section 3.

Estimates in the L 2 -norm
We derive an L 2 error estimate using a duality argument in this subsection. To this end, we consider the following linear dual problem to the primal nonlinear problem eq. (9): for f dual ∈ L 2 (Ω). For smooth domains Ω, it can be deduced by a classical elliptic regularity result that χ ∈ H 4 (Ω). The corresponding weak form is: find that is to say, Remark 2.12. The first equality in eq. (44) holds since u ∈ H 2 (Ω), χ ∈ H 2 (Ω) and v ∈ H 2 (Ω).
We give two auxiliary results in the following.
We are ready to derive the L 2 a priori error estimates. (48) Proof.
We bound each U i separately using the boundedness of A s h , P s h , B s and C s , theorem 2.11 and standard interpolation estimates. This leads to Setting e 3 = u h − u and estimating U 4 as in R 4 of lemma 2.9 with the use of theorem 2.11 and |||I h χ||| h Combining the above estimates for U i (i = 1, 2, 3, 4) and using the regularity estimate eq. (45) and χ 4 I h u − u h 0 , we obtain Hence, eq. (48) follows from the triangle inequality and standard interpolation estimates. theorem 2.15 implies that for quadratic approximations to the sufficiently regular solution of eq. (9), there is a suboptimal convergence rate in the L 2 -norm while for higher order (≥ 3) approximations, we expect optimal L 2 error rates. We shall see numerical verifications of this in the subsequent sections.

The inconsistent discrete form
The above analysis considers the consistent weak formulation eq. (13). In practice, Xia et al. [31] adopted the inconsistent discrete weak form in the implementation due to its cheaper assembly cost: find u h ∈ W h,b such that ComparingÃ s h and A s h , the missing terms are the interior facet integrals arising from piecewise integration by parts and symmetrisation. Due to the absence of these terms inÃ s h , one can immediately notice that the discrete weak formulation eq. (49) is inconsistent in the sense that the solution u of the strong form eq. (9) does not satisfy the weak form eq. (49), as opposed to the result of theorem 2.3.
Despite this inconsistency, in practice this also leads to a convergent numerical scheme with similar convergence rates, as illustrated in section 3. This is not surprising; a similar idea has also been applied and introduced as weakly overpenalised symmetric interior penalty (WOPSIP) methods in [10] for second-order elliptic PDEs and in [8] for biharmonic equations.

A priori error estimates for (P2)
Problem (P2) is a special form of the classical LdG model of nematic LC. Finite element analysis for a more general form using conforming discretisations has been studied in [13,14]. More specifically, Davis and Gartland [14] gave an abstract nonlinear finite element convergence analysis where an optimal H 1 error bound is proved on convex domains with piecewise linear polynomial approximations, but do not derive an error bound in the L 2 norm. Recently, Maity, Majumdar and Nataraj [21] analysed the discontinuous Galerkin finite element method for a two-dimensional reduced LdG free energy, where optimal a priori error estimates in the L 2 -norm with exact solutions in H 2 are achieved for a piecewise linear discretisation. Both works only focus on piecewise linear approximations. In this section, we will follow similar steps to section 2.1 to prove the H 1 -and L 2 -convergence rates for the problem (P2) with the use of common continuous Lagrange elements of arbitrary positive degree. Since the approach is similar to the previous subsections, we omit some details for brevity.
The continuous weak formulation of (P2) in two dimensions (the three-dimensional case can be tackled similarly) is given by: find Q ∈ H 1 b (Ω, S 0 ) such that N n (Q)P := A n (Q, P) + B n (Q, Q, Q, P) + C n (Q, P) = 0 ∀P ∈ H 1 0 (Ω), where the bilinear forms are A n (Q, P) := K Ω ∇Q . . . ∇P, C n (Q, P) := −2l Ω Q : P, and the nonlinear operator is given by Since eq. (50) is nonlinear, we need to approximate the solution of its linearised version, i.e., find Θ ∈ H 1 0 (Ω) such that where ·, · represents the dual pairing between H −1 (Ω) and H 1 0 (Ω). Suppose Q h ∈ V h approximates the solution of eq. (50) with the conforming finite element method on a finite dimensional space V h := {P ∈ H 1 (Ω) : P ∈ Q deg (T ), deg ≥ 1, ∀T ∈ T h }. Throughout this subsection we take deg ≥ 1. Furthermore, we denote V h,0 := {P ∈ V h : P = 0 on ∂ Ω} and V h,b := {P ∈ V h : P = Q b on ∂ Ω}. We assume that the minimiser Q to be approximated is isolated, i.e., the linearised operator DN n (Q)·, · is nonsingular. This is equivalent to the following continuous inf-sup condition [21, Eq. (2.8)]: With this inf-sup condition for DN n (Q)·, · , we can obtain a stability result for the perturbed bilinear form DN n (I h Q)·, · by following similar steps as in the proof of theorem 2.8.
We give some auxiliary results about the operators A n (·, ·), B n (·, ·, ·, ·) and C n (·, ·) that can be verified via the Cauchy-Schwarz inequality, the Poincaré inequality, and Sobolev embeddings. Lemma 2.19. (Boundedness of B n (·, ·, ·, ·), C n (·, ·)) For Ψ, Φ, Θ, Ξ ∈ H 1 (Ω), there holds and for Ψ, Φ ∈ H k (Ω), k ≥ 2, Θ, Ξ ∈ H 1 (Ω), To proceed to error estimates for the nonlinear problem eq. (50), we define the nonlinear map ψ : We again employ an Aubin-Nitsche duality argument to derive L 2 error estimates. To this end, we consider the following linear dual problem to the primal nonlinear problem eq. (10): find N ∈ H 1 0 (Ω) such that for a given G ∈ L 2 (Ω) (we will make a particular choice for G in the proof of theorem 2.27). The weak form of eq. (58) is to find N ∈ H 1 0 (Ω) such that To derive the L 2 a priori error estimates, we need two more auxiliary results.
Finally, we are ready to deduce an optimal L 2 error estimate.
Theorem 2.27. (L 2 error estimate) Let Q be a regular solution of the nonlinear weak problem eq. (50) and Q h be the approximate solution to the discrete problem (having the same weak formulation as eq. (50)). Then We will verify these results in the next section.

NUMERICAL EXPERIMENTS
The proceeding section presents some a priori error estimates for both Q and u in the decoupled case q = 0. We now test the convergence rate of the finite element approximations by the method of manufactured solutions (MMS) and experimentally investigate the coupled case q = 0 in two dimensions. To this end, we choose a nontrivial solution for each state variable and add an appropriate source term to the equilibrium equations, thus modifying the energy accordingly. We can then compute the numerical convergence order.

Test 1: on the unit square
In this test, the numerical runs are performed on the unit square Ω = (0, 1) 2 and we take the following exact expressions for each state variable, Then, in conducting the MMS, we are to solve the following governing equations subject to Dirichlet boundary conditions for both u and Q and a natural boundary condition for u. Here, source terms s 1 , s 2 and s 3 are derived by substituting eq. (62) to the left hand sides, and t 1 and t 2 are given by t 1 := (Q 11 + 1/2)∂ 2 x u + (−Q 11 + 1/2)∂ 2 y u + 2Q 12 ∂ x ∂ y u, t 2 := ∂ 2 x (u (Q 11 + 1/2)) + ∂ 2 y (u(−Q 11 + 1/2)) + 2∂ x ∂ y (uQ 12 ).
We partition the domain Ω = (0, 1) 2 into N × N squares with uniform mesh size h = 1 N (N = 6, 12, 24, 48) and denote the numerical solutions by u h , Q 11,h and Q 12,h . The numerical errors of u and Q in the · 0 -, · 1 -and |||·||| h -norms are defined as The convergence order is then calculated from the formula log 2  Regarding the density variation u, we first present the convergence behaviour of the consistent discrete formulation eq. (13) with penalty parameter ε = 1, since we have proven the optimal error rate in the mesh-dependent norm |||·||| h . Optimal rates are observed in the |||·||| h -norm. Furthermore, optimal orders of convergence in the L 2 -norm are shown for approximating polynomials of degree greater than 2, while a suboptimal rate in the L 2 -norm is given for piecewise quadratic polynomials, exactly as expected. Sub-optimal convergence rates for quadratic polynomials were also illustrated in the numerical results of [29]. We also tested the convergence with the penalty parameter ε = 5 × 10 4 and found that the discrete norms are very similar to We next give the error rates for the inconsistent discrete formulation eq. (49). We illustrate the discrete norms and the computed convergence rates in table 3 with the penalty parameter ε = 1. It can be observed that only first order convergence is obtained in the H 2 -like norm |||·||| h even with different approximating polynomials. Moreover, we notice by comparing tables 2 and 3 that the convergence rate deteriorates slightly for polynomials of degree 3 (although not for degree 4). This, however, can be improved by choosing a larger penalty parameter, as shown in table 4 with ε = 5 × 10 4 , where optimal rates are shown for the discrete norms |||·||| h , · 1 and · 0 for all polynomial degrees (except only sub-optimal in · 0 when a piecewise quadratic polynomial is used as the approximation). The inconsistent discrete formulation appears to be a reasonable choice when a sufficiently large penalty parameter is used. We next investigate the numerical convergence behaviour in the coupled case, i.e., q = 0, in this subsection. Its analysis remains future work, but since it is the coupled case that is solved in practice it is important to assure ourselves that the discretisation is sensible. For brevity, we fix the model parameter q = 30.
We examine the inconsistent discretisation for u with penalty parameter ε = 5 × 10 4 . In unreported preliminary experiments, we observed that the error in Q is governed by the lower of the degrees of the polynomials used for Q and u. We thus give the convergence rates for u and Q separately in tables 5 and 6, with the other degree fixed appropriately. It can be seen that Q retains optimal rates in both the H 1 and L 2 norms, and though there are some fluctuations of the order for u, it still possesses very similar convergence rates when compared with the decoupled case described in Convergence rates for Q with q = 30 and penalty parameter ε = 5 × 10 4 with the inconsistent discretisation eq. (49) for u, fixing the approximation for u to be with the Q 3 element.

Remark 3.2.
We also tested the convergence with the consistent weak formulation for u under the same numerical settings as in tables 5 and 6. We found that in both cases they present very similar convergence behaviour and thus we omit the details here.

Test 2: on the unit disc
For the second set of experiments, we provide numerical results for an exact solution u e with only H 3 -regularity, instead of C ∞ as in eq. (62). Our goal is to investigate whether the H 4 -regularity assumption on u can be relaxed. To this end, we consider a triangular mesh of Ω = {(x, y) | x 2 + y 2 < 1} and choose u e to be u e = (x 2 + y 2 ) 3/2 , and choose the same exact solution for Q e 11 and Q e 12 as in eq. (62). The exact solution given by eq. (63) is in H 3 (Ω) but not in H 4 (Ω), hence violating the regularity assumption of the analysis in section 2.1.
The resulting convergence rates are reported in tables 7 and 8. Table 7 shows that optimal H 1 and L 2 rates are achieved for Q with three different choices of finite elements [P 1 ] 2 , [P 2 ] 2 , [P 3 ] 2 . Table 8 shows the convergence behaviour with penalty parameter ε = 1 when using the inconsistent discrete formulation eq. (49). In contrast to table 3, only first order convergence is obtained for the discrete norm |||·||| h and second-order convergence for · 0 and · 1 , with both P 3 and P 4 . Interestingly, table 8 indicates no convergence when using P 2 elements. It appears that the assumption u ∈ H 4 (Ω) is necessary for our analysis, and that a different analysis should be carried out when this assumption no longer holds.