{'fraction_non_alphanumeric': 0.03315458389192534, 'fraction_numerical': 0.011781994431267402, 'mean_word_length': 4.281726579520697, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hydrodynamics of domain growth in nematic liquid crystals

We study the growth of aligned domains in nematic liquid crystals. Results are obtained solving the Beris-Edwards equations of motion using the lattice Boltzmann approach. Spatial anisotropy in the domain growth is shown to be a consequence of the flow induced by the changing order parameter field (backflow). The generalization of the results to the growth of a cylindrical domain, which involves the dynamics of a defect ring, is discussed.


I. INTRODUCTION
Liquid crystals [1] are an ideal material for the study of topological defects due to the complex textures they create which are easily visible to the naked eye. As topological defects arise in many situations the observed phenomena in liquid crystals can be used to test theories in other areas of physics from cosmic strings [2] to vortices in superfluid helium [3]. Contrary to the assumption inherent in most previous studies of defect dynamics in liquid crystals [4][5][6], in a recent Letter [8] we found that backflow, the coupling between the order parameter and the velocity fields, has a significant effect on the motion of defects. In particular, the defect speed can depend strongly on the topological strength in two dimensions and on the sense of rotation of the director about the core in three dimensions.
These defects were free, in the sense that they were in an unbounded system. However, it is much easier to study liquid crystals experimentally in a confined system. A very straightforward example is the geometry used for display devices. In such a display the liquid crystal is sandwiched between two plates. As the optical and electrical responses of the liquid crystal are coupled, one can apply an electric field between the two plates and directly observe the behaviour.
The operational state of many such devices, including pi-cells [9], is topologically distinct from its state at zero voltage. Before the device can be used, the operational state must be nucleated and grow to fill the display. (A typical cross section through a domain wall separating an operational state and a zero voltage state is shown in Fig. 2(b).) These interfaces between topologically distinct states can behave differently than nematic-isotropic interfaces studied by other groups [7]. Recent experimental work on pi-cells [9] has shown an unusual anisotropy in the domain growth: one side of a domain grows faster than the opposite side. This would appear to be very similar to the anisotropy observed in our simulations of free defects [8]. However there are important differences in this system related to the wall tilt angle, which can dom-inate the defect-defect interaction energy, and the driving force of the electric field. As such, in order to unambiguously characterize the observed anisotropy we need to directly study the growth of these domain walls-which incorporate topological defects-in a confined geometry.
In additional to the technological applications, similar devices have been proposed as ideal experimental realizations of two-dimensional (2D) Ising models (in the plane defined by the walls of the device). Coarsening of reverse tilt domains in liquid crystal cells with heterogeneous alignment layers has been shown to be consistent with predictions of the random-bond Ising model [10] thus providing experimental confirmation of theoretical predictions for domain growth under conditions of quenched random disorder.
The dynamics of a liquid crystal medium is often modeled by using the Ericksen-Leslie-Parodi equations of motion [1,20]. These equations describe the state of the liquid crystal in terms of a director field n which is related to the orientation of the typically long, thin, rodlike molecules which make up the liquid crystalline material. The Ericksen-Leslie-Parodi equations are restricted to an uniaxial order parameter field of constant magnitude. Thus they cannot model the dynamics of topological defects where in the defect core the magnitude of order has a steep gradient and the order parameter field is biaxial. (However, they provide a good description of the bulk away from the defect core.) In order to describe the hydrodynamics of topological defects correctly, we use the more complex Beris-Edwards formulation of liquid crystal hydrodynamics [11]. The propensity to order, as well as the direction along which the system orders are conveniently described by a tensor order parameter Q [1]. The Beris-Edwards equations allow for variations in the magnitude of the nematic order parameter as well as biaxiality present in defect cores. They model both defect dynamics and the coupling between the velocity field and the motion of the order parameter. We use a recent lattice Boltzmann algorithm [13] which has been shown to successfully model the full Beris-Edwards equations.
Our aim in this paper is to study the growth of a do-main of a nematic liquid crystal at the expense of a second domain with a different director orientation. Defects form at the walls between domains and their dynamics is vital in controlling the rate of growth. We find that a spatial anisotropy in domain growth can result from backflow and discuss how the wall speed varies with the material parameters of the liquid crystal, such as viscosity and elastic constants, the geometry and the surface properties of the confining cell and a external electric field. The Beris-Edwards equations of motion are presented in Section II, and the results of the domain growth are described in Section III and Section IV. In Section V we discuss the relevance of our results to the experiments on pi-cells [9]. An outline of the numerical algorithm is given in the Appendix.

II. THE HYDRODYNAMIC EQUATIONS OF MOTION
We summarize the formulation of liquid crystal hydrodynamics described by Beris and Edwards [11], extended to include an electric field and surface potentials. Similar models have been examined by a number of researchers [12]. The continuum equations of motion are written in terms of a tensor order parameter Q which is related to the direction of individual moleculesm by Q αβ = m αmβ − 1 3 δ αβ where the angular brackets denote a coarse-grained average. (Greek indices will be used to represent Cartesian components of vectors and tensors and the usual summation over repeated indices will be assumed.) Q is a traceless symmetric tensor. Its largest eigenvalue, 2 3 q, 0 < q < 1, describes the magnitude of the order.
We first write down a Landau-de Gennes free energy which describes the equilibrium properties of the liquid crystal [1,14] (1) f bulk describes the bulk free energy For γ = 2.7 there is a first-order transition from the isotropic to the nematic phase. The minimum of f bulk describes a uniaxial nematic with an order parameter of the form Q αβ = q(n α n β − 1 3 δ αβ ) where q is zero in the isotropic phase and has a finite value in the nematic phase, and n is the director field.
f el is the analogue of the Frank elastic free energy density This can be easily mapped to give the Frank elastic constants K 1 , K 2 and K 3 [11]. In particular, the "one elastic constant" approximation, K 1 = K 2 = K 3 corresponds to L 1 > 0 and L 2 = L 3 = 0. For a uniaxial nematic, the dielectric constant is anisotropic measured along or perpendicular to the director. The relation between the electric displacement D and field E is of the form [1] More generally, the dependence of the dielectric constant on the order parameter is described by giving results consistent with Eqn.(4) for the uniaxial nematic. The electric contribution to the thermodynamic potential f f ield is The equation of motion for the nematic order parameter is [11] (∂ t + u · ∇)Q − S(W, Q) = ΓH (9) where Γ is a collective rotational diffusion constant. The first term on the left-hand side of equation (9) is the material derivative describing the usual time dependence of a quantity advected by a fluid with velocity u. This is generalized by a second term S(W, Q) = (ξA + Ω)(Q + I/3) where A = (W + W T )/2 and Ω = (W − W T )/2 are the symmetric part and the anti-symmetric part respectively of the velocity gradient tensor W αβ = ∂ β u α . S(W, Q) appears in the equation of motion because the order parameter distribution can be both rotated and stretched by flow gradients. This is a consequence of the rod-like geometry of the liquid crystal molecules. ξ is a constant which depends on the molecular details of a given liquid crystal.
The term on the right-hand side of Eqn.(9) describes the relaxation of the order parameter towards the minimum of the free energy. The molecular field H which provides the driving motion is related to the derivative of the free energy by where and δ αβ is the Kronecker delta. We work in a twodimensional cross section, assuming that the order parameter does not change in the perpendicular direction (although the director may point out of the simulation plane). In addition, the symmetry and zero trace of Q is exploited for simplification. At the surfaces of the device we assume a pinning potential We typically take Q 0 of the form where q is set to the equilibrium bulk value. This corresponds to specifying a preferred direction n 0 for the director at the surface. There can be other terms in the surface free energy [15] and a complete treatment of surface dynamics can be quite rich [16]. However, in this paper we will be operating in the strong pinning limit (α S large) so that the only effect of the pinning potential is to furnish an almost fixed value of Q αβ at the surface (equal to Q 0 ). In all cases studied here, the results are insensitive to the precise value of α S , so long as it is large enough to be in the strong pinning limit. The fluid momentum obeys the continuity and the Navier-Stokes equation where ρ is the fluid density and η = ρτ f /3 is an isotropic viscosity (which is controlled by the simulation parameter τ f described in the Appendix). The form of this equation is not dissimilar to that for a simple fluid. However the details of the stress tensor reflect the additional complications of liquid crystal hydrodynamics. There is a symmetric contribution and an antisymmetric contribution These additional terms can be mapped onto the Ericksen-Leslie equations to give the Leslie coefficients [13]. The background pressure P 0 is constant in the simulations to a very good approximation (±1%). The differential equations for order parameter field Eqn. (9) and the flow field Eqn. (18) are coupled. The velocity field and its derivatives appear in the equation of motion for the order parameter Eqn. (9). Unless the flow field is zero, u = 0, the dynamics given by Eqn. (9) are not relaxational and hydrodynamics play an important role. Conversely, the order parameter field affects the dynamics of the flow field through the stress tensors (19) and (20) which appear in the Navier-Stokes equation (18) and depend on Q and H. This back-action of the order parameter field on the flow field is usually referred to as backflow. To study these equations we use a lattice Boltzmann algorithm summarized in the Appendix. Other than when explicitly stated, the simulation parameters are those listed in [17].

III. DOMAIN GROWTH
We consider a liquid crystal confined between two planes a distance L x apart. The director field may take topologically distinct states depending on the boundary conditions and applied voltage. In the simulations we set the boundary condition (Q 0 in Eq. (15)) so as to give a tilt angle −θ p between the director and the y axis at x = 0 and +θ p at x = L x . At zero applied voltage these conditions result in a global minimum free energy state with a splayed director configuration, or horizontal (H) state as shown in Figure 1(a). At high voltages, typically on the order of 6V , the H state is no longer the global minimum, and a bend configuration (vertical state) is obtained like the one shown in Figure 1 As the H and V states are topologically distinct, the transition from V to H requires nucleation of H domains and the generation of defects. The problem we will investigate is the growth (or shrinking) of the H state within the V state. In particular, we are interested in how hydrodynamics affects the speed of the domain walls. This is partly motivated by the observation in Ref. [9] that the domain growth in a liquid crystal device can be anisotropic and the speculation that this may be due to hydrodynamics.
We have previously observed that the velocity of defects in unbound systems can be affected by hydrodynamics. In particular the defect speed can depend strongly on the topological strength in two dimensions and on the sense of rotation of the director about the core in three dimensions [8]. The crucial difference between the domain growth problem and the motion of free defects is that in the latter case each defect moves due the director field of the other. In the domain growth problem the defects are not interacting but are dragged by the free-energy-driven movement of the domain walls. The free defects are accelerated as they approach each other while in the domain growth problem the defects move with a constant speed. Due to these differences and the additional geometrical parameters, a separate analysis is needed for the confined system which is also easier to realize experimentally and provides a better control of the parameters influencing the defect speed. In order to study the role of hydrodynamics in the system we will examine the factors affecting the domain growth so that we can clearly identify what causes the wall speed anisotropy. The key parameters are the surface director tilt θ p , the sample thickness L x and material parameters such as coefficients in the bulk free energy (1): A 0 , γ, and elastic constants L 1 , L 2 and L 3 . In addition, the rotational diffusion constant Γ, which appears in the dynamical equation (9) for the order parameter gives an overall (inverse) time scale and is related to the Leslie-Ericksen viscosities [13]. For simplicity, we will first study the undriven case of an H domain growing at the expense of a V state. We will later examine growth under the influence of an electric field. The initial configuration, depicted in Fig.  2(a), is a horizontal (i.e., along the y-direction) domain in an otherwise vertically aligned state. This models a time shortly after the electric field has been switched off when small but macroscopic domains have formed in the device. As the simulation proceeds, the director configuration relaxes rapidly to that shown in Fig. 2(b). During the relaxation defects are formed at the center of each domain wall with strengths + 1 2 and − 1 2 , respectively. Once the two defects have formed the vertical domain begins to grow and the + 1 2 and − 1 2 defects move in opposite directions.
Our simulations correspond to a two-dimensional cross section of the two line defects, assuming that the order parameter does not change in the perpendicular direction (although the director may point out of the simulation plane). The two defects are topologically distinct only in two dimensions, but even in three dimensions they are usually separated by an energy barrier.
A particular advantage of the simulations is that the backflow can easily be switched off by setting σ αβ = −P 0 δ αβ and τ αβ to zero. (Compare this to (19) and (20).) Since there is no flow imposed, there is a zero velocity field throughout the whole simulation. The dynamical equation in this case can be obtained from (9) by setting u to zero. It corresponds to the purely relaxational Ginzburg-Landau model [18] ∂ t Q = ΓH (21) where the molecular field H is given by (11). Comparing the dynamics obtained from the Ginzburg-Landau model and the full hydrodynamic equations, the effect of the backflow can be unambiguously identified. The Ginzburg-Landau equation (21) with a single elastic constant is invariant under a local coordinate transformation mirroring the director on the x axis. This corresponds to the transformation The order parameter fields of the two moving defects with topological charges s = ± 1 2 shown in Fig. 2(b) (even including the deformation due to the boundaries) transform into each other. Thus approaches based on a simple Ginzburg-Landau equation predict that as the defects move they follow symmetric dynamical trajectories.
We can construct a simple model for the domain motion in the absence of hydrodynamic flow. In the bulk regions (away from the domain walls), f bulk is minimized with an uniaxial order parameter of the form Q αβ = q(n α n β − 1 3 δ αβ ). We can then restrict our attention to the elastic free energy f el . With this form of the order parameter and if L 2 = L 3 = 0, the elastic free energy density has the form f el = L 1 q 2 (∇θ) 2 if the director remains in the plane so that n = (cos θ, sin θ, 0). The minima in the bulk regions (away from the domain walls) correspond to the director angle changing linearly along the x coordinate from −θ p to +θ p in the H domain, and from −θ p to (+θ p − 180 deg) in the V domain. The difference of the free energy densities of the two domains can then be written as For θ p < 45 deg, the horizontal domain grows because this decreases the free energy of the system. For θ p > 45 deg the horizontal domain should begin to shrink, and at θ p = 45 deg the two domains have the same free energy and the defects should stop moving.
If the H domain grows by a length of ∆L y then the free energy of the system decreases by ∆f × L x × ∆L y . Simple relaxational arguments then suggest that the speed of domain growth can be described by the formula where η e is an effective viscosity. Surface tilt: We first investigate the effect of the surface tilt θ p on the defect speed. Equations (23) and (24) suggest that as θ p increases, the free energy difference decreases, and the defects should move more slowly. The defect speed v is plotted as a function of surface tilt θ p in Fig. 3. Consider first the diamonds. These correspond to the case with backflow switched off. For this case both defects move at the same speed (but in opposite directions). Notice that the defect velocity is proportional to (θ p − 45 deg). From (23) and (24) this leads to the conclusion that the effective viscosity η e is independent of the tilt angle and its value is found to be 0.138 Pa s for the parameters of the simulation. Back-flow: The triangles and circles in Fig. 3 show the velocity of the defects when backflow is included in the model. The s = + 1 2 defect is considerably speeded up, whereas the s = − 1 2 defect is only slightly accelerated. The defect speed remains proportional to (θ p − 45 deg) within a 2% error. The effective viscosities are η +1/2 = 0.083 Pa s < η −1/2 = 0.123 Pa s. These values are comparable to the rotational viscosity γ 1 = 0.08 Pa s [9]. The speed anisotropy defined as is independent of θ p and its value is α = 36%. The order parameter field affects the flow field through the symmetric and antisymmetric stress tensors, (19) and (20). The total nonviscous stress (i.e. the combination of all the stress terms not related to the velocity gradient tensor) is the sum of three terms Here σ i,αβ = −P 0 δ αβ is the stress due to the isotropic pressure. σ d,αβ = −∂ β Q γν δF δ∂αQγν is the deformation stress. For L 2 = L 3 = 0 the deformation stress is σ d,αβ = −L 1 T r(∂ α Q∂ β Q). The rest of the terms in (19) and (20) give what we will refer to as the molecular field stress, σ H , which is a function of H and Q. σ d and the diagonal σ i do not change under the transformation (22) that transforms the defects of topological charge ±1/2 into each other. Conversely, the off-diagonal elements of σ H have their sign inverted. Thus the stress fields and the resultant backflow are different for the two defects.
The stress field σ d is related to the deformation free energy density, which is the same for both defects. It induces vortices similar to those around a solid cylinder moving in a viscous liquid. The flow points in the direction of defect motion at the defect core. The contribution σ H describes the stress due to the the reorientation of the director. The reorientation is the strongest around the core while molecules near the surfaces reorient much less. The director reorientation induces vortices around the two defects as shown in Fig. 2(b). The direction of these vortices is determined by the gradient of the director angle taken moving around the defect in the positive direction. It is positive (negative) for the + 1 2 (− 1 2 ) defect.  Fig. 2(b). There is a strong vortex pair around the s = + 1 2 defect which, at the defect core, points in the direction of defect movement. The flow at the core of the s = − 1 2 defect is weaker and points opposite to the direction of defect propagation.
The two contributions to the backflow reinforce each other for the s = + 1 2 defect but partially cancel for the s = − 1 2 defect. The resulting flow fields can be seen in Fig. 4. The flow is stronger around the s = + 1 2 defect.
Around the s = − 1 2 defect the flow is much weaker and the flow field points opposite to the defect propagation at the core. However, even in this case backflow accelerates the relaxational dynamics. Sample dimensions: As the sample becomes wider the speed of the defect propagation decreases [19] as can be seen in Fig. 5(a). The dependence of the defect velocity on L x follows from Eqns. (23) and (24) The effective viscosity can be calculated from η e = 1/(2πsL 1 Γ) (∇θ q ) 2 dr, where θ q is the field due just to the defect itself [5,6] and the integral is over the volume of the system. Due to the confining geometry, one expects the integral to be dominated by the near-field contribution (i.e. the field near the defect core) which is the same for both a static or moving defect [5]. The gradient of θ caused by a static defect drops off as 1/r and hence one expects the effective viscosity to go like log(L x /L x0 ) [6]. This dependence is observed in Fig. 5(b). The fit is η e = 0.62 Pa s × log(L x /L x0 ) where L x0 = 0.076µm is comparable to the defect core diameter. The dashed line in Fig. 5(a) shows the fit to Eq. (27) including the L x dependence of the viscosity.
When backflow is considered, the relative speed anisotropy increases with L x and saturates at about 60% as shown in the inset of Fig. 5(a). The increase is probably due to the increasingly larger regions around the core involved in director reorientation. This leads to stronger vortices and hence stronger hydrodynamic effects. Electric Field: When an external electric field is applied, it changes the free energy densities of the H and V domains. Thus it also influences the speed of the domain growth. Fig. 6 shows the speed of the two defects as the function of applied voltage. For low voltages the free energy difference between the H and V domains can be estimated. If we assume that the orientation of the H and V states are unchanged from the zero voltage case (i.e. the director angle remains a linear function of x, as used in Eq. (23)), then substituting this into Equation (8) and integrating over space gives for a sample of length L y . The total free energy difference between the domains is the sum of the elastic and the field contributions, Eqn. (23) and Eqn. (28), respectively. Substituting this into Eqn. (24) results in a parabolic dependence of the velocity on V for low voltages as shown in Fig. 6. At V limit ∼ 0.9V domain growth is reversed since the effect of the surfaces is balanced by the influence of the electric field. At high voltages (V ∼ 0.6V ), the H domain is replaced by H A , an asymmetric horizontal domain shown in Fig.  7. It has a lower free energy than the H domain due to the more favorable alignment with the electric field. As a result of the deformation of the H A domain, the defects move towards the surfaces. For V > 1V the molecules in the bulk are almost completely aligned with the vertical field. The horizontal region is then confined to a thin layer near the surfaces. For this type of configuration, the free energy difference between the domains is proportional to the voltage [9]. This gives the linear slope of the curve in Fig. 6 seen for the higher voltages.
For low voltages the speed anisotropy is 36% and independent of the voltage. As shown in the inset of Fig. 6, at high voltages (V > 1V ) the speed anisotropy decreases. The effect on the anisotropy can be explained by the relative weight of the relaxational dynamics and the backflow. In the V = 1V to 1.6V range the relaxational dynamics are substantially speeded up by the increasing voltage due to the electric field contribution (8) in the free energy. The backflow is induced by the the stress fields given in Eq. (19) and Eq. (20). These stress fields do not depend directly on the electric field, only on the order parameter field, which changes only slightly in this voltage range. Therefore the stress fields do not increase with the increasing voltage, and the hydrodynamics is dominated by the relaxational dynamics at high fields. (In comparison, for the experiment presented in [9] the domain wall speed was ∼ 0.1µm/ms and the anisotropy also decreased with increasing voltage.)

IV. OTHER CONTROL PARAMETERS
The equations governing the dynamics of liquid crystals covered in Section II contain a large set of parameters. In this section we explore some of this parameter space. In particular, we examine the case of multiple elastic constants, important for comparing to any real liquid crystal. In addition we look at the influence of the different viscosities and free energy parameters on the balance of the relaxational dynamics and the backflow. We also examine the case of asymmetric and inhomogeneous surface tilts since they give new insights about the underlying symmetries of the system and are important for practical devices with non-trivial surface anchoring.
Elastic constants: Real liquid crystals have multiple elastic constants. We now consider the effect of non-zero elastic constants L 2 and L 3 . If L 2 = 0 in the expression of the molecular field Eqn. (13), the dynamical equation in the absence of backflow Eq. (21) loses its invariance under the transformation (22). However, the speed anisotropy is very small. The reason for this is that the relaxational dynamics are still invariant under the mirroring transformation for a uniaxial order parameter with a constant magnitude [22]. In our setup these conditions hold except for close to the defect cores. A larger anisotropy in speed is obtained if L 3 = 0. If L 3 > 0 (L 3 < 0) then the s = + 1 2 (s = − 1 2 ) defect is faster. For L 1 =8.73pN and L 3 =15.87pN, we measure a speed anisotropy of α = 3%, in a system without hydrodynamics. This anisotropy due to the unequal elastic constants increases with increasing sample thickness. For L x = 1.25µm (our benchmark system has L x = 0.7µm [17]) the anisotropy due to non-zero L 3 is α = 6%. This may be due to the fact that, for free defects, the elastic anisotropy causes significant deviations from the case of isotropic elasticity in the order parameter field away from the axis determined by the two defect cores [1]. In the thinner sample, the surfaces "cut off" this part of the field, decreasing the anisotropy. If the surface tilt is close to vertical, and the horizontal domain is shrinking, then the L 3 dependence of the speed anisotropy is the opposite. If L 3 < 0 (L 3 > 0) then the s = + 1 2 (s = − 1 2 ) defect is faster. Since in the two cases (growing vs. shrinking domain) the order parameter fields near the axis of the two cores are the same, this should also be attributed to the differences in the order parameter fields far from the defect cores which results from the different tilts.
Viscosities and diffusion: Consider now the effect of the parameters governing the time scales in the equation of motion for the domain growth. τ f is proportional to the viscosity in the Navier-Stokes equation (18). Increasing τ f increases the viscosity, slows the defects, and decreases their velocity anisotropy as shown in Fig. 8. The velocity tends to that measured without backflow, as represented by a dashed line in the figure. This is as expected since backflow will be suppressed by a large viscosity.
Increasing Γ which appears in Eq. (9) increases the velocity of both defects, but decreases the relative speed anisotropy as shown Fig. 9. The defects move faster, due to the fact that Γ governs the speed of the relaxational dynamics. Since the stress fields in Eq. (19) and Eq. (20) do not depend on Γ, they do not increase. As a result, as Γ increases, the relaxational dynamics speed up, but the backflow does not change as significantly, and as a result the anisotropy decreases.
Free energy: γ is the parameter in the free energy (2) which controls the magnitude of order in the bulk of the domains. (The isotropic-nematic phase transition is at γ = 2.7.) When γ decreases, the defect core gets bigger and this results in a smaller effective viscosity [5]. Thus the defects move faster under the relaxational dy-namics. The decrease in the magnitude of order results in a smaller backflow due to the reorientation and α decreases. Changing A 0 in the free energy (2) does not affect the homogeneously aligned bulk H and V domains, only the defect structure. The larger A 0 , the larger the energy cost of any deviation from the magnitude of order corresponding to the minimum of f bulk in (2). Decreasing A 0 increases the size of the defects, and as above, the defects move faster under the relaxational dynamics. Since the defect core size increases, the magnitude of order around the core decreases resulting a smaller backflow due to reorientation and hence α decreases. The effect of increasing L 1 is similar. It increases the defect core size, speeds up the defects, and leads to a smaller velocity anisotropy.
Non-symmetric surface tilt: The director surface tilt at the top and bottom surfaces does not have to be symmetric. In this more general case, tilt angles at the surfaces can be written as a sum of a symmetric and an antisymmetric contribution: The dynamical equations for L 2 = L 3 = 0 and without flow are invariant under the local rotation of all the molecules by the same angle. The dynamics for a nonsymmetric surface tilt can therefore be obtained from the symmetric case by rotating all the molecules by θ p,a . Thus even for non-symmetric tilts the two defects move with the same speed and remain in the center between the two plates. The full dynamical equations however, are not invariant under the rotation of molecules if flow is included (or if L 2 = 0, L 3 = 0). Thus the defects will move with a different speed. Moreover they are no longer constrained by symmetry to lie midway between the two plates. Typical director and velocity fields for an asymmetric surface tilt are shown in Fig. 10.
It is also possible to construct a patterned alignment of liquid crystals on surfaces [23]. If the surface tilt is not homogeneous then, when the defect arrives at a region with a different tilt, it assumes the new velocity corresponding to this tilt. On the boundary of two bend domains with opposite surface tilt the defect can even "change" topological strength by merging with the two defects located at the surface as shown in Fig. 11.

V. GROWTH OF A CYLINDRICAL DOMAIN: THE HYDRODYNAMICS OF A DEFECT RING
In this section we consider the three dimensional analogue of our system, where a liquid crystal is held between parallel plates ∼ µm apart. A domain nucleated at a point will grow to a cylindrical shape with its axis perpendicular to the confining plates. Based on our experience with the two dimensional system we can discuss how backflow can affect the growth of a cylindrical H domain in a V environment.
At the domain boundary there is a defect ring, as shown in Fig. 12(a). The defect configuration of a vertical cross section through the ring (perpendicular to the plates) changes gradually from a − 1 2 to a + 1 2 defect. Although this problem is three-dimensional, a vertical cross section indicated by dashed lines in Fig. 12(a) gives a geometry similar to that considered in Section IV.

V V (a)
Rubbing direction For simplicity, assume that the H domain is a perfect cylinder. In this case the director field of any vertical cross section passing through the middle of this cylindrical domain can be obtained from our two-dimensional simulation plane by rotating the director field locally by a given angle around the vertical axis, x.
Let us now examine the effect of the backflow. σ d and the diagonal σ i do not change during the local rotation of the tensor order parameter field Q around the vertical axis, x, by the same angle. The generated velocity vortices always lie in the plane of the cross sections across the defects ring. These generate flow which is always in such a direction as to expand the defect ring, indicated by straight arrows in Figs. 12(a-c).
The tensor σ H does however change under a local rotation around x. For the cross section which includes a + 1 2 defect the resultant backflow points in the direction of defect motion. For the − 1 2 defect, the flow points opposite to the defect motion, as shown in Fig. 12(b). In both cases the vortices are in the plane of the cross section. Fig. 12(c) shows the cross section indicated by the dashed-dotted lines in Fig. 12(a). Now the directors and the vortices due to the reorientation are in a plane perpendicular to the cross section. Thus the total flow field will vary around the domain wall and will lead to anisotropy in the domain growth.
An experimental setup similar to this was considered by Acosta et al. in their investigation of domain growth and switching in pi-cell liquid crystal devices [9]. The growth of a horizontal domain in a bend (V) or twisted bend environment was studied: such a transition is needed to produce the operational state of the device. (The twisted bend configuration has a lower energy than the bend state for small surface tilts and low voltage, if the Frank elastic constant K 2 is sufficiently small. The twisted bend state is replaced by a bend state for larger (∼ 30 deg) surface tilts.) A cylindrical bend or twisted bend domain was formed in a H environment and the domain wall velocity was measured at four points around the ring, where its cross section corresponds to a + 1 2 and a − 1 2 defect, and at two points halfway between these. It was found that the wall at the + 1 2 defect moved substantially faster than at the other three. It seems very plausible that the essential physics is captured by our model.
Further measurements of defect dynamics in confined geometries have been done very recently [24] and these techniques should allow further testing of the concepts we present here.

VI. SUMMARY
In this paper we explored domain growth in nematic liquid crystals. Defects form at moving domain walls. We find that a wall incorporating a s = + 1 2 defect is substantially speeded up by backflow effects, whereas a wall containing a s = − 1 2 defect is only slightly affected. This was explained in terms of the symmetry properties of the different stresses acting on the defects. These reinforce each other for the s = + 1 2 defect while partially cancelling for the s = − 1 2 defect. The influence of different material and geometrical parameters on the velocity anisotropy was interpreted by comparing the relative weight of the relaxational dynamics and the backflow. By generalizing two-dimensional simulation results, a qualitative picture was proposed for the role of the backflow in three dimensions.
Results were obtained using a lattice Boltzmann algorithm to solve the Beris-Edwards equations of liquid crystal hydrodynamics. Working within the framework of a variable tensor order parameter it was possible to correctly incorporate variations in the magnitude of order and hence the dynamics of domain walls and their associated topological defects.

VII. ACKNOWLEDGMENT
We would like to thank E.J. Acosta, C.M. Care, S. Elston, K. Good, N.J. Mottram, E. Orlandini and T.J. Sluckin for helpful discussions. We acknowledge the support of Sharp Laboratories of Europe at Oxford. C.D. acknowledges funding from NSF Grant No. 0083286 and G.T. from the EPSRC Grant No. M04426 .

APPENDIX A: A LATTICE BOLTZMANN ALGORITHM FOR LIQUID CRYSTAL HYDRODYNAMICS
We now summarize a lattice Boltzmann algorithm which solves the hydrodynamic equations of motion of a liquid crystal (9), (17), and (18). Lattice Boltzmann algorithms are defined in terms of a set of continuous variables, usefully termed partial distribution functions, which move on a lattice in discrete space and time [21].
The simplest lattice Boltzmann algorithm, which describes the Navier-Stokes equations of a simple fluid, is defined in terms of a single set of partial distribution functions which sum on each site to give the density. For liquid crystal hydrodynamics this must be supplemented by a second set, which are tensor variables, and which are related to the tensor order parameter Q [13].
We define two distribution functions, the scalars f i ( x) and the symmetric traceless tensors G i ( x) on each lattice site x. Each f i , G i is associated with a lattice vector e i . We choose a nine-velocity model on a square lattice with velocity vectors e i = (±1, 0), (0, ±1), (±1, ±1), (0, 0). Physical variables are defined as moments of the distribution function The distribution functions evolve in a time step ∆t according to G i ( x + e i ∆t, t + ∆t) − G i ( x, t) = ∆t 2 [C Gi ( x, t, {G i }) + C Gi ( x + e i ∆t, t + ∆t, {G * i })] .
This represents free streaming with velocity e i and a collision step which allows the distribution to relax towards equilibrium. f * i and G * i are first order approximations to f i ( x + e i ∆t, t + ∆t) and G i ( x + e i ∆t, t + ∆t) respectively. They are obtained by using only the collision operator C f i ( x, t, {f i }) on the right of Equation (A2) and a similar substitution in (A3). Discretizing in this way, which is similar to a predictor-corrector scheme, has the advantages that lattice viscosity terms are eliminated to second order and that the stability of the scheme is improved.
The collision operators are taken to have the form of a single relaxation time Boltzmann equation [21], together with a forcing term The form of the equations of motion and thermodynamic equilibrium follow from the choice of the moments of the equilibrium distributions f eq i and G eq i and the driving terms p i and M i . Full details of the algorithm can be found in [13].