A Moving Mesh Method for Modelling Defects in Nematic Liquid Crystals

The properties of liquid crystals can be modelled using an order parameter which describes the variability of the local orientation of rod-like molecules. Defects in the director field can arise due to external factors such as applied electric or magnetic fields, or the constraining geometry of the cell containing the liquid crystal material. Understanding the formation and dynamics of defects is important in the design and control of liquid crystal devices, and poses significant challenges for numerical modelling. In this paper we consider the numerical solution of a $\bf{Q}$-tensor model of a nematic liquid crystal, where defects arise through rapid changes in the $\bf{Q}$-tensor over a very small physical region in relation to the dimensions of the liquid crystal device. The efficient solution of the resulting six coupled partial differential equations is achieved using a finite element based adaptive moving mesh approach, where an unstructured triangular mesh is adapted towards high activity regions, including those around defects. Spatial convergence studies are presented using a stationary defect as a model test case, and the adaptive method is shown to be optimally convergent using quadratic triangular finite elements. The full effectiveness of the method is then demonstrated using a challenging two-dimensional dynamic Pi-cell problem involving the creation, movement, and annihilation of defects.


Introduction
The orientational properties of liquid crystal materials can be manipulated by applying an electric or magnetic field, leading to particular characteristics of the reflection and transmission of light waves. These effects make liquid crystals key materials in the construction of a broad range of commonly-used display devices, such as the Twisted Nematic Device (TND) [33], the Pi-cell [8] and the Zenith Bistable-Device (ZBD) [10,30]. More recently, there has been growing interest in liquid crystals in a wider context. Examples include active liquid crystals [26] (which are relevant to natural applications such as modelling cytoskeletal structure in cell biology or animal flocking as well as in synthetic manufacture of colloids and granular matter), liquid crystal shells and drops [23], and materials design and self-assembly of ordered fluids [37]. Because of their importance in these and other technological applications, there is a great deal of interest in modelling the properties of liquid crystals mathematically.
The most commonly-used continuum models utilise one or more unit vector fields as state variables. For the uniaxial nematic phase, which is the simplest and most common liquid crystal phase, the orientation of the molecules is represented by a unit vector denoting the direction in which their main axis points. This is known as the liquid crystal director and is traditionally denoted by n. More generally, taking n and −n to be equivalent, the average molecular orientation can be represented by an order tensor, usually denoted by Q. This tensor can be written as where S and T are scalar order parameters, {l, m, n} is a set of orthonormal directors, and I is the identity (see, for example, [36]). Note that the uniaxial case can be recovered by setting T = 0. In this paper, we propose an efficient numerical method for computing the orientational state of a nematic liquid crystal based on a Q-tensor model. In particular, we focus on tracking the movement of defects in the material, that is, local regions (of point, line or wall type) where the symmetry of the ordered material is broken. The switching behaviour of liquid crystal material between two equilibrium states (by means of an applied field), which is the basis of most liquid crystal devices, is strongly influenced by the existence of such defects, so it is important to be able to model these features accurately. Our use of the Q-tensor (as opposed to a director-based) model in this paper is driven by the fact that in this formulation, topological defects do not appear as mathematical singularities.
A Q-tensor theory of nematic liquid crystals, which allows for changes in the scalar order parameters, has been developed from the theory of Landau by de Gennes [13]. Minimisation of the total free energy in the case of a nematic liquid crystal coupled with an applied electric field leads to a set of six coupled partial differential equations (PDEs) for the five degrees of freedom of the order parameter tensor Q and the electric potential U , which poses a challenge for numerical solvers. Furthermore, additional physical features such as flow and temperature change require the Q-tensor equations are coupled to the Navier-Stokes and energy equations. Even in the absence such additional complications, the Q-tensor equations are difficult to solve numerically due to their highly non-linear nature. Also, the defects mentioned above induce distortion of the director over very small length scales as compared to the size of the cell. It can therefore be difficult to accurately represent their nature and behaviour with a standard numerical model. The large discrepancies in length and time scales which occur mean that numerical difficulties are even more acute for models of dynamic problems involving the movement of defects, such as the Pi-cell problem studied in §5.2.
For identifying static equilibrium states, relatively straightforward numerical methods are often good enough (see, for example, [12,18,22,29,34,35]). There have also been several studies using more sophisticated adaptive techniques. These include the h (grid parameter) and p (degree of basis function) adaptive finite element methods presented in [11,16,17,21]. Additional methods have been proposed based on moving meshes [1,2,24,25,31,32], which move existing mesh points so as to cluster them in areas of large solution error whilst maintaining the same mesh connectivity. These techniques are particularly appropriate for resolving localised solution singularities such as defects, as maintaining a fixed connectivity is very efficient in terms of computing time (as opposed to adding or removing grid points in areas of interest). Also, for transient problems, it is sometimes possible to use bigger time steps if the solution remains almost stationary relative to the moving mesh frame of reference. This motivates our use of adaptive moving mesh techniques to capture defect structure and track defect movement within the cell.
In [25], we proposed a robust and efficient numerical scheme for solving the system of six coupled partial differential equations which arises when using Qtensor theory to model the behaviour of a nematic liquid crystal cell under the influence of an applied electric field in one space dimension. The numerical method uses a moving mesh partial differential equation (MMPDE) approach to generate an adaptive mesh which accurately resolves important solution features. In this paper, we extend this adaptive moving mesh strategy to solve liquid crystal problems in two dimensions. This involves addressing a number of significant new challenges, including the choice of appropriate adaptivity criteria for problems with moving singularities, the efficient solution of the large systems of highly non-linear algebraic equations arising after discretisation, and how to deal with the creation and annihilation of defects in a realistic model.
The remainder of the paper is structured as follows. In §2, we give a brief overview of the derivation of the physical PDEs arising from the Q-tensor framework coupled with an applied electric field, along with some details of their finite element discretisation when an adaptive moving mesh is utilised. In §3.1, the details of the two-dimensional moving mesh PDE are given. We consider a number of different mesh adaptivity criteria through the use of monitor functions, and present a series of numerical experiments which indicate that monitor functions based on a local measure of biaxiality perform well. We then apply the biaxiality-based monitor function to a problem first presented by Bos [40]: a two dimensional Pi-cell problem with a sinusoidal perturbation across the centre of the cell. This is a dynamic two-dimensional version of the test problem described in [25, §1.2].
2 Derivation and discretisation of physical PDEs

Derivation of physical PDEs
To characterise the molecular alignment of a nematic liquid crystal, we define a uniaxial Q-tensor using a local ensemble average of the molecular axes as (see, e.g., [14, §2.1.2]). The unit vectors u lie along the molecular axes and the angle brackets denote the ensemble averaging: the factor 3/2 is included for convenience so that, for a uniaxial state with director n and scalar order parameter S, tr(Q 2 ) = S 2 . The tensor (2) has five degrees of freedom and is symmetric and traceless, so it can be represented in matrix form as where each of the five quantities q i , i = 1, . . . , 5, is a function of the spatial coordinates and time. Note that the orthonormal eigenvectors of this matrix are the vectors {l, m, n} used in the representation of the Q tensor given in (1). The globally stable state of a nematic liquid crystal under the influence of an applied electric field corresponds to a minimum point of the free energy. Using Landau-de Gennes theory, in which the free energy density is assumed to depend on Q and its gradient, the free energy may be written as where F t , F e , F u and F s represent the thermotropic, elastic and electrostatic terms, respectively. Note that, as here we only consider problems with fixed (strong anchoring) boundary conditions, we omit any (constant) surface energy terms. Expressions for the individual terms in the integrand of (4) can be derived in a variety of different ways: here we expand the thermotropic energy, F t , up to fourth order in Q and the elastic energy, F e , up to second order in the gradient of Q. Details of the resulting expressions can be found in [25, §2], along with a description of the contribution from the applied electric field, F u , which includes flexoelectricity. As in [25], values for material constants throughout this paper are those used in [3], which are commensurate with a liquid crystal cell of the 5CB compound 4-cyano-4 -n-pentylbiphenyl.
To derive time-dependent PDEs for the quantities q i in (3), we use a dissipation principle with viscosity coefficient ν and dissipation function where the dot represents differentiation with respect to time (see, e.g., [36, eq. (4.23)]). For a physical domain with spatial coordinates {x 1 , x 2 , x 3 }, this produces a system of equations which can be written as involving the bulk energy F b = F t + F e + F u , where the vectorΓ i has entries andf i is given byf To add the coupling with an electric field, E say, we introduce an additional unknown in the form of a scalar electric potential U such that E = −∇U .
Assuming that there are no free charges, the electric field within the cell can then be found by solving the Maxwell equation where the specific form of the electric displacement D can be found in [25, eqn (2.5)]. Minimisation of the total free energy (4) therefore leads to a set of six coupled non-linear PDEs for the five degrees of freedom of Q and the electric potential U . Specifically, combining (5) and (6) and using some algebraic manipulation for notational convenience, we obtain the equations where and Note that we non-dimensionalise the equations (7) for computational purposes: details of the precise scalings applied in terms of lengths and energies are given in [25, §2].
The governing physical PDEs in (7) now have to be adapted to account for the movement of the finite element mesh. To do this, we introduce a family of bijective mappings such that, at a given time t in time domain T ⊆ R + , the point ξ = (ξ, η) of a two-dimensional computational reference domain Ω c is mapped to the point x = (x, y) of the original physical domain Ω. The temporal derivative of a mapping g : Ω → R (from the physical domain) in the computational frame can then be defined as and applying the chain rule for differentiation (with appropriate smoothness assumptions on g) gives which includes an additional convection-like term due to the mesh movement. Incorporating these changes into (7a) gives the final set of six coupled PDEs

Finite element discretisation
With a space of time-independent finite element test functionsv ∈ H 1 0 (Ω c ), mesh mapping (8) defines the test space We denote the approximation spaces with essential boundary conditions on q i and U by H Eq and H E U , respectively. In an analogous way to the onedimensional case described in [25, §3.1], Reynolds's transport formula can be used to derive the following conservative weak formulation of (9): find q i ∈ H Eq (Ω), i = 1, . . . , 5, and U ∈ H E U (Ω) such that ∀v ∈ H 0 (Ω) To discretise (10), we assume that the reference domain Ω c is covered by a fixed triangulation T h,c with straight edges, so that Ω c = ∪ K∈T h,c K, and introduce the Lagrangian finite element spaces where P k (K) is the space of polynomials on K of degree less than or equal to k. Using a piecewise linear discretisation of the mesh mapping (8) to produce a discrete mapping A h,t ∈ L 1 (Ω c ), finite element spaces on the physical domain Ω can be defined as (again, analogously to the one-dimensional setting studied in [25, §3.2]). Letting H h,Eq ⊂ L k (Ω) and H h,E U ⊂ L k (Ω) be the finite dimensional approximation spaces satisfying the Dirichlet boundary conditions for the q i 's and U , respectively, the finite element spatial discretisation of the conservative weak formulation (10) is therefore: Finally, introducing vectors q i (t) and u(t) which contain the degrees of freedom defining q ih and U h , respectively, (11) can be rewritten to obtain the highly nonlinear differential algebraic system where M (t) is the (time-dependent) finite element mass matrix.
3 Derivation and discretisation of moving mesh PDEs

Equations governing mesh movement
Having formulated equations to represent the physical PDEs, we now establish a mechanism for moving the mesh: this will be done by constructing so-called moving mesh PDEs. To avoid potential mesh crossings or foldings, we derive a suitable evolution equation for the inverse mapping (8) (see, for example, the discussion in [15]). A mesh T h,t on Ω can then be generated as the pre-image of a fixed mesh T h,c on Ω c . As introduced in [19], we choose the mapping ξ(x) corresponding to a fixed value of t in order to minimise the functional where G is a 2 × 2 symmetric positive definite matrix referred to as a monitor matrix, and ∇ is the gradient operator with respect to x. Rather than directly attempt to minimise (13), a more robust procedure is to evolve the mapping according to the modified gradient flow equations Here, τ > 0 is a user-specified temporal smoothing parameter which affects the temporal scale over which the mesh adapts, and P is a positive function of (x, t), chosen such that the mesh movement has a spatially uniform time scale [20]. The selection of an appropriate monitor matrix is crucial to the success of mesh adaptation. In this paper, we will consider the monitor matrix proposed by Winslow [39] where w(x, t) is a positive weight function called a monitor function. The choice of the monitor function should ideally be based on a local a posteriori error estimate but if no such estimate exists then the monitor function can be any smooth function designed to adapt the mesh towards important solution features. Suitable choices for w for applications to the Q-tensor equations are discussed below. In practice, we interchange the roles of the dependent and independent variables in (14), since it's the location of the physical mesh points {x i (t)} N i=1 that defines the mapping A t . With a Winslow-type monitor matrix (15) the resulting MMPDEs take the form where a = 1 w and J = x ξ y η − x η y ξ is the Jacobian of A t . To complete the specification of the coordinate transformation, the MMPDE (16) must be supplemented by suitable boundary conditions g(ξ, t), ξ ∈ ∂Ω c ; these are obtained using a onedimensional moving mesh approach. The numerical solution of (16) requires spatial and temporal discretisation. We discretise in space using standard linear Galerkin finite elements. For time discretisation, we use a backward Euler integration scheme to update the solution at t = t n+1 and, to avoid solving nonlinear algebraic systems, we evaluate the coefficients a, c, . . . , e at the time t = t n . We therefore seek x n+1 The resulting linear systems are solved using the iterative method BiCGSTAB [38] with an incomplete LU (ILU) factorization [28] as a preconditioner. An analysis of the performance of this iterative solver for the discretised MMPDE equations can be found in [7].

Details of the monitor functions
An appropriate choice of monitor function w(x, t) in (15) is essential to the success of any adaptive moving mesh method. In this paper we consider twodimensional analogues of the monitor functions which were shown in [25] to be appropriate for one-dimensional Q-tensor models. We first describe these assuming that we have an input function T (x, t) which represents a physical quantity associated with the particular problem under consideration: a discussion of appropriate input functions for our problem involving finite element approximation of the Q-tensor matrix follows in §3.3.
We consider three different forms of monitor function: • AL. This is based on a measure of the arc-length of T : • BM1. This is a generalisation of the Beckett-Mackenzie monitor function introduce in [4,5], based on first-order partial derivatives of T : where scaling parameters α and m are discussed below.
• BM2. This is a second variant of the Beckett-Mackenzie monitor function which takes into account information about second-order partial derivatives of T : again involving scaling parameters α and m.
The value of the parameter α in (19) and (20) is determined a posteriori from the numerical approximation itself as for BM2. Its purpose is to avoid mesh starvation external to layers in the solution as, without it, the resulting mesh would have almost all mesh points clustered inside the layers due to the monitor function being very small elsewhere. The lower bound on α in (3.2) removes unwanted oscillations in the mesh trajectories caused by amplification of errors in approximating T (x, t) which could otherwise cause the mesh to adapt incorrectly. The parameter m in (19) and (20) also helps to regulate mesh clustering: when m > 1, any large variations in T (x, t) are smoothed so that mesh points are more evenly distributed over the domain. In [5] it is shown that, for a function in one dimension with a boundary layer, the optimal rate of approximation order using polynomial elements of degree p can be obtained by ensuring that the parameter m ≥ p + 1. With no specific guidance on the choice of m in higher-dimensional settings, we follow our work in [25] and choose m = 3. In general, the monitor function often has large spatial and temporal variations, so to improve the robustness of the moving mesh method we employ both spatial and temporal smoothing procedures. This results in an MMPDE that is easier to integrate forward in time and a smoother mesh, which can improve spatial solution accuracy. Temporal smoothing is done by under-relaxing the monitor function so that the monitor function at the current time level n is given by where 0 < ω < 1 is an under-relaxation parameter (in this paper, we set ω = 0.8). Following [20], spatial smoothing of the monitor function is done by taking a local average of the monitor function across cells that have a common vertex. That is, the smoothed monitor functionw is defined as where x m ∈ Ω is a mesh point in the physical domain, ξ m ∈ Ω c is the corresponding a mesh point in the computational domain, and C(ξ m ) ⊂ Ω c represents all neighbouring cells of vertex ξ m . If required, spatial smoothing can be repeated in an iterative fashion to further smooth the monitor function.

Details of the monitor input functions
Having specified monitor functions, it remains to decide on an appropriate input function T (x, t). We will consider two variants: a. Order parameter. We recall from §2.1 that, for a uniaxial state with scalar order parameter S, we have S 2 = tr(Q 2 ). This has led to the function being used to generate monitor functions in previous studies [1,2,31,32]. This quantity is known to vary rapidly in regions where order reconstruction occurs, and was shown in [24] to be ideal for certain one-dimensional uniaxial problems.
b. Biaxiality. For biaxial problems, an alternative input function (based on the direct invariant measure of biaxiality used in [3]) is This takes values ranging from 0 (for a uniaxial state) to 1 (for a wholly biaxial state).
In the numerical experiments in §5, we compare the performance of the AL, BM1 and BM2 monitor functions with various input functions. Details of the specific combinations highlighted in the results presented are summarised in Table 1

Solution algorithm
The numerical algorithm for solving the full problem involves an iterative solution strategy as originally proposed in [7,6]. Full details of how it can be used in a Q-tensor setting are given in [25, §5], so are not reproduced here. Instead, we simply highlight the main features and point out any modifications needed for application in two dimensions. The iterative solution algorithm involves completely decoupling the solution of the physical PDE system (12) from the solution of the MMPDE (16). This has a key advantage in that different convergence criteria can be used for the two systems: it is well understood that the computational mesh rarely needs to be resolved to the same degree of accuracy as the solution of the physical PDEs. The system for the Q-tensor components (12a) is integrated forward in time using a second-order singly diagonally implicit Runge-Kutta (SDIRK2) method, with the electric potential values being updated by solving (12b) at each step of the Newton iteration used to generate intermediate stages (see [25, §5.2]). We also use an adaptive time-stepping procedure when integrating forward in time. This is based on both the computed solutions of (12a) for q i and on the solution of the MMPDE (16). The specific two-dimensional error indicator used here is where N E is the number of mesh elements and ∆ n j denotes the j th element of the mesh at time level n. The individual error terms

Numerical results
In this section, we illustrate the performance of our novel time-adaptive method with a moving two-dimensional finite element mesh using two test problems involving Q-tensor models of liquid crystal cells. Note that we use quadratic finite elements on triangualr meshes throughout.

Test problem 1: resolving the core structure of a stationary defect
We begin by considering the resolution of a stationary liquid crystal defect, as this problem is ideal for assessing the ability of particular monitor functions to resolve the core defect structure. Specifically, we assume we have a disclination line in the z-direction and that far from the defect core the director n lies in the x-y plane. Such a defect can be simulated by imposing the initial condition on the director where d i is referred to as the disclination index. We note that on travelling a closed path around the disclination line, the director rotates through the angle 2πd i . Here we restrict our attention to the case d i = 1/2: a plot of the steady state director field is shown in Figure 1(a), where we observe that the director rotates through an angle of π radians as the centre of the defect is circled. Using numerical simulations of a Q-tensor model, Schopohl & Sluckin  [34] found that the defect core structure was contained in a circular region of radius approximately 5ζ, where ζ is the nematic coherence length. For this example, we solve the Q-tensor equations on the square region [−10ζ, 10ζ] 2 , with ζ ≈ 4.06 nm (as commensurate with our use of the physical parameters in [3]). On the domain boundaries, we impose Dirichlet conditions corresponding to the director being given by (25) and the order parameter S taking a value associated with the equilibrium (nematic) phase (S = S eq 0.65 with our parameters). An analytical solution of the Q-tensor equations does not exist for this problem, so a reference solution was obtained using a fine adaptive mesh using the BM2b monitor function with N = 5334 quadratic triangular elements, using a time step ∆t = 10 −8 until a steady state solution was obtained. We consider this reference solution to be independent of the specific choice of the monitor function as calculations using reference solutions based on the other three monitor functions gave very similar results. Figure 1(b) shows the three (numerically computed) eigenvalues of the Qtensor along the line y = 0 at the centre of the cell. We observe that an exchange of eigenvalues takes place at the centre of the core region, as the material passes through a biaxial transition (corresponding to the switch from horizontally to vertically aligned directors along y = 0 in Figure 1(a). Contour plots of the computed order parameter S = tr(Q 2 ) 1/2 (cf. (23)) and the biaxiality (as measured by (24)) are shown in Figure 2. We can see that the order parameter takes  its equilibrium value S ≈ 0.65 outside a central circular region of diameter 10ζ. Furthermore, within this region the order parameter varies significantly within a core of diameter 2ζ. We can also see that, outside a core of diameter of approximately 4ζ, the biaxiality is zero, and inside it has a volcano-like structure with a rim of value 1 representing the purely biaxial state, and a base with value 0 representing the uniaxial state: this sudden variation takes place over a core only a few nanometres in diameter, so is very difficult to capture accurately using a numerical method.

Estimated rate of spatial convergence
It is important to check that approximations obtained on a sequence of adaptive meshes are convergent as we increase the number of mesh elements. In [24], we presented convergence results for a scalar model of a one-dimensional uniaxial problem, with a much more complicated one-dimensional Pi-cell order reconstruction problem, modelled by a full Q-tensor model, being considered in [25].
In a similar vein, we now consider the rate of spatial convergence of the moving mesh finite element approximation for this fully two-dimensional defect model problem.
In what follows, we use q i * (x, t) to denote the reference approximation to the exact solution q i (x, t), and q iN to denote the finite element approximation calculated on a grid with N quadratic elements. We assume throughout that The error in the approximation q iN is denoted by e N qi . Since the approximate solution grid points will not in general coincide with the reference grid points, the reference solution at a fixed time t = t * , q i * (x jk , t * ), is interpolated (using the MATLAB function TriScatteredInterp [27]) onto the approximate solution grid. The spatial error in the l ∞ norm is then estimated using the maximum error computed at the grid nodes, that is, We fix time t * = 0.2 ms as by this time the solution has entered a steady state. All approximate solutions are obtained using the BM2b monitor function. The error norm (26) for the non-zero components of Q (components q 2 and q 5 are exactly zero for this problem) for the +1/2 defect problem are plotted in Figure  3(a) for various values of N . We observe that e N qi l∞ appears to converge at  the rate O(N −3 ) which is the optimal rate expected using quadratic triangular elements.

Resolving the defect core
An example of an adapted mesh using N = 1388 elements and the BM2b monitor function is shown in Figure 3(b). Although it is clear that the mesh has been adapted isotropically towards the core of the defect, at this scale it is difficult to observe any detail of exactly how the adaptivity resolves the defect core structure. To give some insight into the resolution obtained using the different monitor functions, Figures 4 and 5 show cross-sections (taken along the line y = 0) of the order parameter and biaxiality, respectively (plotted as solid blue lines). The location of grid nodes for the monitor and input function  Table 1 are also plotted (as red circles) in each case to help visualise how each method copes with adapting to the nano-structure of the defect core. All of the meshes are clearly adapting to resolve the core structure of both the order parameter and the biaxiality. However, we note that the BM2b monitor function in particular has placed a significant number of nodes exactly at the defect core, right inside the volcano structure coming from the biaxiality. We know from our previous experience with the one-dimensional Picell problem [25] that it is particularly difficult to resolve this structure, but the BM2b monitor function is still doing a good job here for the full two-dimensional problem.
In addition to the accuracy of approximations produced, we must take into account the computational cost using each monitor function. The plots in Figure 6 show the l ∞ error (26) for the three non-zero components of the Q-tensor (q 1 , q 3 , q 4 ) plotted against the total CPU time in seconds required for each method. Also included here for comparison are the results with uniform meshes with the same number of elements. The first important observation from these results is that, regardless of which monitor function is used, the MMPDE-based adaptive methods always outperform a standard uniform mesh in terms of this measure of efficiency. Furthermore, the results also show clearly that, as more accurate solutions are sought, the BM2b monitor function proves to be the most cost-effective choice, as in the cases other methods are cheaper, the error is unattractively large. Although the more traditional arc-length based monitor function (AL) comes closest to matching the accuracy of BM2b, it does so at a far greater cost. Hence, overall, we conclude that the BM2b combination of monitor and input functions is the method of choice. Figure 6: The l ∞ error in q 1 , q 2 , q 3 plotted against the total CPU time in seconds for each method, measured at time t * = 0.2 ms. The data points correspond to grids using 122, 162, 218, 286, and 342 quadratic triangular elements.

Test problem 2: two-dimensional Pi-cell problem
We now consider a fully two-dimensional time-dependent problem, involving defects which move in time through a liquid crystal cell, eventually annihilating each other to leave an unperturbed state. The geometry is that of a Pi-cell [9] of width two microns and thickness one microns, with liquid crystal parameters again taken from [3]. At both boundaries, the cell surface is treated so as to induce alignments uniformly tilted by a specified tilt angle, θ t , but oppositely directed. If a sufficiently high voltage is applied across the Pi-cell for long enough, then a transition from the splay state (which has mostly horizontal alignment of the director with a slight splay, as depicted in Figure 7(a)) to the bend state (which has mostly vertical alignment of the director with a bend of almost π radians) can be achieved. Based on experimental results, two different physical mechanisms for this transition have been proposed: a homogeneous transition via the material melting in the central plane of the Pi-cell, or an inhomogeneous transition mediated by the nucleation of defect pairs which move and eventually annihilate each other. The homogeneous transition problem is essentially one-dimensional and has previously been modelled by several authors using moving mesh techniques [1,2,24,32]. For a more challenging test of our two-dimensional adaptive moving mesh approach, we will concentrate here on the simulation of the inhomogeneous transition type with moving defects. This problem is still in theory relatively unchallenging: at t = 0, if no perturbation is applied, the director angle simply varies linearly between the tilt angles, as in Figure 7(a), with a director angle across the middle of the cell of θ = 0 • . In practice, however, it is unrealistic for this to be achieved exactly in a physical cell due to small variations in the pretilt angles or thermal fluctuations. We therefore follow [40] and modify the initial director angle across the middle of the cell so that it follows the sinusoidal function sin(2πx/p), where x is the spatial  coordinate in the horizontal plane, and p is the cell width. This perturbation is fixed only at t = 0 for one time step, but introduces solution gradients in two dimensions, as portrayed in Figure 7(b), which provide a bigger challenge for our numerical method.
We consider a cell of width 2µm and thickness 1µm, with a pre-tilt angle of θ = ±6 • . An electric field of strength 18V µm −1 is applied parallel to the cell thickness at time t = 0. Based on the evidence from the static defect problem in §5.1, we present results from the BM2b monitor function only.
Initially, immediately before the application of the electric field, the cell is in an equilibrium state where the order parameter and biaxiality take constant values of 0.65 and 0, respectively. The mesh at this stage is quasi-uniform as no adaptivity has yet taken place. As time evolves, the combined effect of the perturbation and the applied electric field become apparent. Figure 8 shows the cell state 12µs after the application of the electric field; at this time there is a region of concentrated splay distortion at the centre of the cell. Within this area, the order parameter and biaxiality are no longer at their constant equilibrium values: the mesh, as expected, has started to adapt as depicted in Figure 8(c). After 15.5µs the distortion at the centre of the cell has become more pronounced, and we can clearly observe pairs of +1/2 and -1/2 defects within this area,   as shown in Figure 9. Outwith the distorted area, the cell is largely in an equilibrium state, with the order parameter and biaxiality still at their constant equilibrium values. However, the cores of the defects are now completely biaxial, and the measure of biaxiality approaches its maximal value of 1. Recalling that the BM2b monitor function uses biaxiality as its input, it is not surprising that the mesh has now adapted significantly from its quasi-uniform initial state, and has started to adapt well to resolve the defects (see Figure 10). As time evolves further, the oppositely signed defects are attracted to each other, moving ever closer until they ultimately meet and annihilate each other. Figures 11 and 12   show snapshots of the order parameter and biaxiality respectively, measured after 15.5µs, 16µs and 17µs, calculated on the meshes shown in Figure 10. In  Figure 11, after 16µs, the mesh is still well adapted to the sinusodial shape of the initial perturbation, consistent with the presence of large variations in the biaxiality throughout the central area of the cell. However, after 17µs, by which point the defects have almost coalesced, the mesh has relaxed in areas where the biaxiality is now back to its equilibrium value, and instead is completely focused on resolving the defects. After the defects meet and annihilate, the biaxiality and order parameter again relax towards their equilibrium value everywhere in the cell, and the mesh also relaxes back to a quasi-uniform state. Overall, throughout the simulation, the adaptive moving mesh method does an excellent job of tracking the development, movement, and annihilation of the defects in the liquid crystal cell. In particular, the method is able to cope well with the small-scale structure of the defect core, and the short timescales associated with the establishment and annihilation of defects.

Summary
The focus of this paper is on the description and application of a new efficient moving mesh method for Q-tensor models of liquid crystal cells. Although some of the ideas contained here are described in a one-dimensional setting in [25], extending the method to tackle the more physically realistic fully two-dimensional problems presented here required us to address a number of significant new challenges. As with all moving mesh methods, the choice of an appropriate adaptivity criterion is crucial: here we have established that using a monitor function based on second-order partial derivatives of a local measure of the biaxiality of the liquid crystal material (BM2b in Table1) is extremely successful in this regard. Using a test problem based on a static +1/2 defect, we demonstrated in §5.1 computed solutions from all of our proposed methods converge optimally in space. However, a comparison of efficiency demonstrated that the BM2b monitor function is clearly the method of choice in terms of computationally efficiency when a reasonable level of accuracy is required. Furthermore, when applied to the more realistic but more challenging fully two-dimensional Pi-cell problem described in §5.2, the adaptive MMPDE method based on the BM2b monitor function proved to be very effective for resolving the movement and core details of defects, including the creation and annihilation of these moving singularities. This is particularly impressive given the very short length and time scales involved in these aspects of the material's behaviour.
Of course, some challenges still remain. Particularly useful in practice would be the extension of our method to multi-dimensional problems with irregular geometries. This would pose a further challenge to the adaptive moving mesh method as it would potentially have to resolve defects present around the areas where the cell geometry is most complex. A prime candidate would be the Zenith Bistable-Device (ZBD) described in [10,30], where the liquid crystal cell has an alignment layer on the upper surface and a periodic grating structure on the lower surface.