Numerical Method of Lines for the Relaxational Dynamics of Nematic Liquid Crystals
A. K. Bhattacharjee, Gautam I. Menon, and R. Adhikari The Institute of Mathematical Sciences,
arXiv:0804.4758v1 [cond-mat.soft] 30 Apr 2008
C.I.T. Campus, Taramani, Chennai 600013, India
(Dated: April 30, 2008)
We propose an e?cient numerical scheme, based on the method of lines, for solving the Landaude Gennes equations describing the relaxational dynamics of nematic liquid crystals. Our method is computationally easy to implement, balancing requirements of e?ciency and accuracy. We benchmark our method through the study of the following problems: the isotropic-nematic interface, growth of nematic droplets in the isotropic phase and the kinetics of coarsening following a quench into the nematic phase. Our results, obtained through solutions of the full coarse-grained equations of motion with no approximations, provide a stringent test of the de Gennes ansatz for the isotropic - nematic interface, illustrate the anisotropic character of droplets in the nucleation regime and validate dynamical scaling in the coarsening regime.
PACS numbers: 64.60.Cn, 61.30.Jf, 82.20Mj, 05.70.Fh
Liquid crystalline phases of matter are rich in examples of subtle order parameters, complex phase behaviour, and exotic defect structures . In the nematic phase of liquid crystals, broken rotational symmetry produces elastic and hydrodynamic modes and topological defects of integer and half-integer charge . Understanding the statics and dynamics of nematics is di?cult because the order parameter is neither a scalar or a vector, but a more complicated tensorial quantity constrained by symmetry and normalisation . For example, the simplest relaxational dynamics which follows from a time-dependent Ginzburg-Landau equation describing a nematic close to thermal equilibrium, entails the solution of a set of 5 coupled non-linear parabolic partial di?erential equations for the independent components of the order parameter tensor Q . For all but the simplest situations, these equations do not have analytical solutions and thus need to be solved numerically. The e?cient numerical computation of the solutions of the relaxational dynamics of nematics - nematodynamics is the problem addressed here. In this paper, we propose an e?cient numerical scheme, based on the method of lines [4, 5], for solving the Landau-de Gennes equations for nematodynamics. The method of lines (MOL) is a powerful technique for discretizing initial-value partial di?erential equations (PDEs). The essence of the MOL is semi-discretisation, where a PDE in spatial and temporal variables is discretised in the spatial variable only. This reduces the PDE to a system of ODEs in the temporal variable. The great advantage of the MOL is that there are powerful numerical methods implemented in general-purpose numerical libraries for the solution of systems of ODEs. The MOL also provides a great degree of freedom in discretising space, allowing one to chose from ?nite-di?erence, ?nite-volume, ?nite-element, spectral collocation, or any other suitable spatial discretisation. The strategy of semi-discretisation, where the spatial discretisation and temporal integration are treated as separate steps, allows for considerable ?exibility, as both these operations can be optimised according to the demands of the problem to ensure maximum accuracy at minimum computational cost. Here we combine a ?nite-di?erence spatial discretisation with a Runge-Kutta temporal integration to solve the relaxational equations of nematodynamics. We benchmark our method through the study of the following three problems: (a) the isotropic - nematic interface and tests of the de Gennes ansatz, (b) the kinetics of the growth of single nematic 2
droplets in the isotropic phase and (c), coarsening kinetics following a quench from the isotropic into the nematic phase. To enable a comparison with previously published results our numerical results are in two dimensions, though the method and our code are applicable in one, two and three dimensions. Our results are summarised as follows. Our numerical scheme permits a stringent test of the de Gennes ansatz for the uniaxial nematic isotropic interface . We show that biaxiality is absent in the interfacial region and that the hyperbolic tangent pro?le proposed by de Gennes provides an accurate ?t to the numerical data. We show next that single nematic droplets placed in an isotropic solvent at isotropic-nematic coexistence coarsen anisotropically , with initially circular droplets deforming into more ellipsoidal structures as time evolves. In our study of coarsening upon quenches from the isotropic phase, we observe the characteristic annihilation of integer and half integer charge defects (disclinations) as revealed by schlieren plots. For the coarsening problem, we ?nd that order parameter correlation functions at di?erent times can be rescaled to lie on a master curve, and that the . Our numerical results coarsening exponent associated with the growing length scales is 1 2 are in excellent agreement with analytical results where available and consistent with numerical results where such result exists. In addition, we extend previous results for each of the problems we study in signi?cant ways. In Section II, we summarise the order parameter description of nematics in terms of the symmetric, traceless order parameter tensor Qαβ . We discuss the Landau - de Gennes free energy functional describing the free energy associated with order parameter con?gurations, the phase diagram for the free energy and the equations of relaxational dynamics. In section III we present the MOL scheme for partial di?erential equations and show how it may be applied to nematodynamics. In section IV we discuss each of the applications described above in some detail, presenting our numerical results. We conclude by describing extensions of the present method to other problems in the dynamics of nematic liquid crystals and compare our method with other methods available in the literature.
Orientational order in the nematic phase is quanti?ed through a second-rank, symmetric traceless tensor, Qαβ . The principal axes of this tensor, obtained by diagonalizing 3
Q, specify the direction of ordering. The principal values represent the strength of ordering. Locally, nematic order is de?ned through Q(x, t) = duf (x, u, t) uu ≡ uu where
f (x, u, t) is the molecular orientational distribution function at the position x at time t which counts the number of molecules with long axis oriented in the direction u, X de?nes the symmetric traceless part of an arbitrary real second rank tensor X with X = 1 1 X + X T ? T r (X ), 2 d (1)
where X T denotes the transpose of the tensor i.e. (X T )αβ = Xβα and d denotes the dimension. The average < · >≡ duf (x, u, t)(·). The three principal axes specify the director n, the codirector l and the joint normal to these, m. Given the two principal values S and T , the order parameter can be written as 3 1 1 Qαβ = S (nα nβ ? δαβ ) + T (lα lβ ? mα mβ ), 2 3 2 where α, β ≡ x, y, z and, 1 S = P2 (cos θ) = (cos2 θ ? ) , 3 2 T = sin θ cos 2φ ,
1 ≤ S ≤ are measures of alignment with ? 3 2 , 3
0 ≤ T < 3S . S =
corresponds to the
nematic phase whereas S = 0 to the isotropic phase. The (polar) angle between n and u is θ, while the (azimuthal) angle between l and u is φ. In the uniaxial nematic phase, T = 0 and the tensor (2) takes the form Q = (3/2)S nn . In the biaxial phase Q = (3/2)S nn + (1/2)T (ll ? mm). The tensor Q is diagonal in a coordinate system aligned with the principal axes, ? ? ?(S + T )/2 0 0 ? ? ? ? Q=? 0 ?(S ? T )/2 0 ? . ? ? 0 0 S
i in a set of basis tensors Tαβ  which are symmetric, traceless and orthonormal, 5
In a frame of reference which is not aligned with the principal axes, Q must be expanded
Qαβ (x, t) =
i ai (x, t)Tαβ .
Explicitly, these are T1 = √ ?? z. T5 = 2 y
? z? z , T2 =
1 (x ? 2
x ??y ?y ? ), T 3 =
2 x ?y ? , T4 =
2 x ?? z
The Landau-Ginzburg free energy functional F  is obtained from a local expansion in powers of rotationally invariant combinations of the order parameter Q(x, t), 1 1 1 Fh [Q] = AT r Q2 + BT r Q3 + C (T r Q2 )2 + E ′ (T r Q3 )2 . 2 3 4 order parameter must be added : 1 1 Fel [? Q] = L1 (?α Qβγ )(?α Qβγ ) + L2 (?α Qαβ )(?γ Qβγ ), 2 2 (8) (7)
To this, non-local terms arising from rotationally invariant combinations of gradients of the
where α, β, γ denote the Cartesian directions in the local frame, and L1 and L2 are elastic constants. Surface terms of the same order in gradients may also be added, but we omit them in the present work. (T r Q2)3 ≥ (T r Q3 )2 , higher powers of cooling transition temperature. From the inequality 1 6 T r Q3 can be excluded for the description of the uniaxial phase. Thus the uniaxial case is described by E ′ = 0 whereas E ′ = 0 for the biaxial phase. For nematic rod-like molecules B < 0 whereas for disc-like molecules, B > 0. The quantities C and E ′ are always taken to be positive to ensure stability and boundedness of the free energy in both the isotropic and nematic phases. The total free energy including local and gradient terms is 1 1 1 1 1 d3 x[ AT r Q2 + BT r Q3 + C (T r Q2 )2 +E ′ (T r Q3 )2 + L1 (?α Qβγ )(?α Qβγ )+ L2 (?α Qαβ )(?γ Qβγ )]. 2 3 4 2 2 (9) In the local free energy density Eq.(7), A = A0 (1 ? T /T ? ), where T ? denotes the super-
The ?rst order isotropic to uniaxial nematic transition at the critical value S = Sc is obtained from the equations, 9 3 2 4 CSc + E ′ Sc , 4 4 9 3 B = ? CSc ? 9E ′ Sc . 2 A = 9 ′2 , A = ? CSc 2 45 ′3 . B = ? E ′ Sc 4 5 (10) (11)
′ The second order uniaxial to biaxial transition at the critical value S = Sc is obtained from,
In the phase diagram of Fig. (1), there are lines separating the nematic phase with S > 0 from the discotic phase in addition to the isotropic to nematic transition lines. There are also two uniaxial to biaxial nematic second order lines as well as the two ?rst order isotropicuniaxial nematic lines. One of these transition lines marks the transition to the nematic phase with S > 0 and the other for the discotic phase with S < 0. The four phases meet at the bicritical Landau point A = B = 0. At the bicritical point, the isotropic to nematic transition is second order rather than ?rst order. In the absence of thermal ?uctuations and hydrodynamic ?ow, the equation of motion of the order parameter can be taken to be relaxational , ?t Qαβ (x, t) = ?Γαβ?ν with 2 Γαβ?ν = Γ[δα? δβν + δαν δβ? ? δαβ δ?ν ]. 3 lessnes of the order parameter. In the above, Γ is a constant. With the choice of the Landau-de Gennes free energy, the equation of motion takes the form ?t Qαβ (x, t) = ?Γ [(A + CT rQ2 )Qαβ (x, t) + (B + 6E ′ T rQ3 ) Q2 αβ (x, t) ?L1 ?2 Qαβ (x, t) ? L2 ?α (?γ Qβγ (x, t)) ]. (16) (15) δF , δQ?ν (14)
The kinetic coe?cient Γαβ?ν ensures that the dynamics preserves the symmetry and trace-
i The equations can be projected onto the tensorial basis Tαβ to give a set of equations for
the independent components ai ,
j i 2 i ?t ai = ?Γ [(A + CT rQ2 )ai + (B + 6E ′ T rQ3 )Tαβ Q2 αβ ? L1 ? ai ? L2 Tαβ Tβγ ?α ?γ aj ]. (17) j i Using the orthonormality of the basis tensors Tαβ Tβα = δij , this can be expanded to the set
of ?ve equations,
a2 1 a2 5 4 2 2 + )? ?t a1 = ?Γ [(A + CT rQ2 )a1 + (B + 6E ′ T rQ3 ) √ (a2 ? a ? a + 1 2 3 2 2 6 1 1 2 1 2 2 L1 ?2 a1 ? L2 ( ?2 a1 + ?z a1 + √ ((?y ? ?x )a2 + ?x ?z a4 + ?y ?z a5 ) ? 6 2 12 1 √ ?x ?y a3 )], 3 2 a2 a2 a1 a2 + √4 ? √5 ) ? ?t a2 = ?Γ [(A + CT rQ2 )a2 + (B + 6E ′ T rQ3 )(? 3 8 8 1 1 2 2 2 2 ? ?x )a1 + ((?x L1 ?2 a2 ? L2 ( √ (?y + ?y )a2 + ?x ?z a4 ? ?y ?z a5 ))], 2 12 2a1 a3 a4 a5 ?t a3 = ?Γ [(A + CT rQ2 )a3 + (B + 6E ′ T rQ3 )(? √ + √ ) ? 6 2 1 1 2 2 L1 ?2 a3 ? L2 (? √ ?x ?y a1 + ((?x + ?y )a3 + ?y ?z a4 + ?x ?z a5 ))], 2 3 a1 a4 a2 a4 a3 a5 ?t a4 = ?Γ [(A + CT rQ2 )a4 + (B + 6E ′ T rQ3 )( √ + √ + √ ) ? 6 2 2 1 1 2 2 L1 ?2 a4 ? L2 ( √ ?x ?z a1 + ((?x + ?z )a4 + ?x ?z a2 + ?y ?z a3 + 2 12 ?x ?y a5 ))], a1 a5 a3 a4 a2 a5 ?t a5 = ?Γ [(A + CT rQ2 )a5 + (B + 6E ′ T rQ3 )( √ + √ ? √ ) ? 6 2 2 1 1 2 2 L1 ?2 a5 ? L2 ( √ ?y ?z a1 + ((?y + ?z )a5 ? ?y ?z a2 + ?x ?z a3 + 2 12 ?x ?y a4 ))].
Note that these equations are parabolic, non-linear, contain anisotropic gradient terms and are coupled strongly to each other. The e?cient computation of the solutions to these equations in a variety of physical situations is described in succeeding sections of this paper.
METHOD OF LINES DISCRETIZATION
The equations for the relaxational dynamics of the orientational order parameter presented in the previous section are a set of ?ve coupled non - linear parabolic di?erential equations. No analytical solutions are available for these equations in general and e?cient and accurate numerical methods must therefore be sought. Below we present such a method, based on the strategy of semi-discretisation, whereby an initial value PDE is discretised only in the spatial variables to yield a set of coupled ODEs. These ODE’s can be solved using powerful general purpose solvers. 7
In the literature, this semi-discretisation strategy goes by the name of the ‘method of lines’ (MOL) [4, 5], since the solutions are obtained along ?xed lines in the space - time plane. To illustrate this, consider a parabolic initial-value problem in a single spatial variable x, the di?usion equation for a scalar ψ (x, t): ?t ψ (x, t) = D ?2 ψ (x, t). (23)
The MOL discretisation proceeds by restricting the ?eld ψ (x, t) to a set of N discrete ‘collocation’ points, xn = n?x, n = 1, 2, . . . , N , spaced ?x apart, and then uses a discrete approximation for the Laplacian based on these points. The simplest of these approximations is based on local polynomial interpolation resulting in ?nite di?erence schemes . However, the MOL is not restricted to ?nite di?erences. When high accuracy is needed, global interpolation based on trigonometric or Chebyshev polynomials can be used to generate spectral approximations of derivatives . When conservation laws need to be respected, ?nite volume approximations to the derivatives can be used . In complicated geometries, a ?nite-element discretisation may be the most appropriate. In this example, the simplest approximation is based on nearest-neighbour ?nite di?erences, for which (?2 ψ )(xn ) = 1 (ψ (xn+1 ) ? 2ψ (xn ) + ψ (xn?1 )) + O ((?x)2 ). (?x)2 (24)
Inserting this into the di?usion equation, we obtain a set of coupled ordinary di?erential equations ?t ψ (xn , t) = 1 (ψ (xn+1 ) ? 2ψ (xn ) + ψ (xn?1 )). (?x)2 (25)
This coupled set of ordinary di?erential equations, together with initial and boundary conditions, can now be integrated with a suitable numerical integration scheme which ensures accuracy, e?ciency, and stability. It is at this stage that the ?exibility of the MOL is most apparent, since any number of numerical integration schemes can be implemented without a?ecting the spatial discretisation. Depending on the nature of the system of ODEs, the optimal choice may be either an explicit scheme, an implicit scheme, or one which is designed to handle sti?ness. In the above example, for instance, it is well known that an explicit Euler integration scheme leads to an instability unless the time step ?t is constrained by the Courant-Friedrichs-Lewy condition ?t < (?x)2 /2D . In contrast, an implicit integration scheme based on the trapezoidal rule gives the stable Crank-Nicolson update . In general,
semi-discretisation followed by numerical quadrature provides an elegant way of deriving many of the well-known schemes for parabolic PDEs. The practical implementation of PDE solvers using the MOL discretisation is simple since the spatially discretised ODE system can be passed directly to a general purpose ODE solver. The complexity, both algorithmic and computational, in the numerical integration can thereby be transferred directly to the ODE solver library. We have followed the MOL discretisation to construct a solver for PDEs describing the dynamics of the orientational order parameter presented in the previous section. We have used standard nearest-neighbour second-order accurate ?nite di?erence formulae for ?rst and second derivatives in the spatial discretisation . More accurate di?erence approximations to derivatives can be obtained using Fornberg’s general formula . For the temporal integration, we have used the ODE solver routines in the GNU Scienti?c Library. For the problems we study here, we ?nd that an explicit multi-stage integrator gives a good compromise between accuracy, stability and computation expense. Finally, it should be mentioned that the MOL discretisation is not restricted to parabolic initial value problems, but also applies to hyperbolic and mixed parabolic - hyperbolic PDEs . This allows the method described in this paper to be extended to situations where e?ects of advection from hydrodynamic ?ow need to be accounted for.
In this section we report benchmarking results using the methodology described in the previous section for the following experimentally relevant situations: the properties of the isotropic-nematic interface at coexistence, the kinetics of droplet evolution between the binodal and spinodal lines and spinodal decomposition into the nematic phase. In each case, we compare our numerical results with available analytical and numerical results in the literature. A suitable non-dimensionalisation of the equations of motion before they are discretised is essential for controlling artefacts and error introduced by the discretisation. For our problem, we non-dimensionalise the governing equations (16) using the values of dimensionful quantities at coexistence (T = Tc ). The resulting dimensionless quantities are superscripted with an asterisk, while the dimensionful quantities are subscripted with c. 9
We non-dimensionalise the order parameter Q? αβ = Qαβ /Sc by the value of the strength then S ? = S/Sc = 3(1 + of ordering at coexistence Sc = ?2B/9C . The non-dimensionalised strength of ordering is 1 ? 24AC/B 2 )/4. The gradient terms in the free energy are non54C (L1 + 2L2 /3)/B 2, which is related to the width
dimensionalised by the length lc =
of the isotropic-nematic interface at coexistence. The free energy is non-dimensionalised as
4 F ? = F/Fc , where Fc = 9CSc /16 is proportional to the free energy barrier between the
isotropic and nematic minima at coexistence. Finally, time is non-dimensionalised by the characteristic relaxation time τ = (ΓFc )?1 . To solve the equations of motion Eq.(17), we choose the discretisation length ?x = ?y = 1 and the integration time-step ?t = 1. These de?ne simulation units of length and time. To ensure that discretisation errors and artefacts are kept to a minimum, the discretisation length must be much smaller than the characteristic length ?x = ?y ? lc . The integration time step must be much smaller than the characteristic time scale ?t ? τ . 1 and ?t? = ?t/τ ? 1. Our simulations are performed maintaining the above conditions, The dimensionless discretisation scales must then satisfy ?x? = ?x/lc ? 1, ?y ? = ?y/lc ? on grid sizes ranging from 32 × 32 to 256 × 256, with periodic boundary conditions, using the 4th-order Runge - Kutta method as implemented in the GNU Scienti?c Library for the ODE solver.
The isotropic-nematic interface at coexistence
The isotropic-nematic interface was studied in a seminal paper  by de Gennes within the framework of a Landau-Ginzburg description. To render the problem analytically tractable, de Gennes made a speci?c assumption regarding the variation of the components of the order parameter across the interface. For an in?nitely extended interface where, by homogeneity, variations perpendicular to the interface alone are allowed, de Gennes assumed that the only quantity which changed across the interface was the strength of ordering S . In the de Gennes ansatz, there is no biaxiality and no variation of the director across the interface. This reduces a problem with ?ve degrees of freedom to a more manageable problem involving only a single degree of freedom. The variation of the ordering strength S along the coordinate z normal to the interface located at z0 can then be obtained analytically as, S (z ) = z ? z0 Sc (1 ? tanh ), 2 w 10 (26)
where w =
2/S ? is the non-dimensional interfacial width.
We have veri?ed this remarkable ansatz with a direct numerical solution. In our numerical calculations, a strip of nematic interface was sandwiched between two isotropic domains with periodic boundary conditions. The system was then allowed to relax to the minimum of the free energy. The parameters were chosen such that the width w ? ?z = 1, ensuring that discretisation errors were kept to a minimum. The resulting pro?les for the variation of S and T are shown in Fig. (2). The values obtained for T are consistent with de Gennes’ assumption of vanishing biaxiality. The variation of S at each of the two isotropic-nematic interfaces were ?tted, using the leastsquares method, to the analytical pro?le, with saturation value of the order Sc , the location of the interface z0 and the interface width w as ?tting parameters. As shown in the inset to Fig. (2), ?tted values of w agree remarkably well with the analytic result for a range of parameter values. The agreement is accurate to within a fraction of a percent. Similar results were obtained for the saturation value of the order parameter. This benchmark clearly demonstrates the accuracy of the MOL scheme in reproducing the equilibrium limit of Eqn.(14). The de Gennes ansatz also predicts that the energies of planar and homeotropic anchoring are identical when elasticity is isotropic (L2 = 0). We compared the value of the free energies of the interface for both planar and homeotropic anchoring of the director. To machine precision, these answers are identical. Our results for this problem represent the ?rst direct veri?cation of the de Gennes ansatz retaining all degrees of freedom of the orientational tensor.
Nematic droplet in an isotropic background
Since the isotropic-nematic transition is ?rst order, there exists a regime of parameters where the kinetics proceeds by nucleation. In the phase diagram of Fig. (1) this regime is bounded by the binodal and spinodal lines. A droplet of nematic phase in an isotropic background shrinks if it is smaller than a critical size Rc . Droplets which are larger than Rc grow in size till they expand to the size of the system. A droplet of size Rc is thus metastable. The development in time of a ?uctuation-induced droplet is a key factor in determining 11
the nucleation rate . Unlike conventional isotropic ?uids, the nucleus in a nematic is not expected to be spherical . Allowing the shape to deviate from perfect sphericity reduces the total energy, which is a sum of the elastic energy associated with the director deformation in the bulk and a surface energy associated with the anchoring condition at the interface of the droplet. Thus, determining the droplet shape of least energy involves a minimisation over the strength of ordering as well as the director degrees of freedom [13, 14, 15]. In our simulations we have studied the time evolution of droplet shapes. For studying conformational changes, our initial condition was chosen to be S = Sc , with the director anchored at an angle of π/4 with the nematic - isotropic interface inside the droplet. We took S = 0 in the isotropic region. We relaxed the system at a temperature intermediate between the binodal and spinodal temperatures. For droplets larger than the critical radius, the nematic region was observed to grow in size, ?nally evolving into the full nematic state. For the parameters used, the alignment of the director did not change signi?cantly during the evolution of the droplet. Figs. (3) shows three stages in the evolution of a single, initially circular droplet, illustrating how the droplet shape evolves into an elongated ellipsoidal form. In the absence of elastic anisotropy L2 = 0, the droplet remains circular. In the presence of elastic anisotropy, the droplet orients along the direction of nematic order for planar anchoring (L2 > 0) and perpendicular to it for homeotropic anchoring (L2 < 0). In our numerics, we quantify this geometrical change through measurements of the change in the aspect ratio. The aspect ratio, de?ned as the ratio of the major to minor axis of the ellipse, are indicated in Figs. (3). This is calculated by extracting a contour at a ?xed value of S (say Sc /2) and ?tting it with an ellipse. To the best of our knowledge, there are no analytical expressions for order parameter variations for nematic droplets. Thus, our results cannot be compared directly against analytical theory. However, the results obtained are in qualitative agreement with previous work on the shape of nematic droplets as a function of anchoring strength [13, 14, 15, 16].
In the previous section, we studied the dynamics of droplets when the quench was to a temperature between the binodal and the spinodal. For a quench below the spinodal 12
temperature, the isotropic phase becomes locally unstable to a nematic perturbation and the system proceeds to the nematic phase by spinodal decomposition. Coherent regions of local nematic order develop in time, with a distinct axis of order in each of these domains. Topological defects form at the intersections of these di?erently ordered domains. Coarsening proceeds through the annihilation of topological defects, increasing the correlation length of local orientational order. To study the coarsening kinetics, we start from a random initial con?guration where we draw the order parameters S and T randomly from a Gaussian random variable with variance proportional to Fc , ensuring that 0 < T < S . We obtain cos θ from an uniform distribution between -1 and 1, and choose φ similarly between 0 and 2π to generate the director, codirector and the joint normal. We then relax the system from this initial condition at a temperature below the supercooling spinodal temperature. The data presented below is averaged over 100 di?erent initial conditions for a 256 × 256 system with periodic boundary conditions. From the coarsening simulations we obtain the strength of ordering, the biaxiality and the director. The director is used to construct the schlieren plots shown in Figs. (4). These plots are constructed by ?rst projecting the director into the x ? y plane, ?nding the angle χ made by this projection with an arbitrary axis (say the x- axis) and then computing sin2 (2χ).
The presence of both integer and half-integer defects is clearly visible in these plots as the meeting points of four and two dark brushes, respectively. In the corresponding plots for the strength of ordering, the defects are clearly visible as localised regions where S rapidly decreases. This is the core region of the topological defect, shown in Figs. (5). We con?rm the surprising ?nding that there is strong biaxial ordering inside the defect core . These results are in perfect qualitative agreement with both theoretical predictions and previous numerical results . To make a quantitative comparison with previous work, we compare results for the timedevelopment of correlation functions during coarsening. We calculate the real-space correlation function, C (r, t) = and its Fourier transform, S (k, t) = Qαβ (k, t)Qβα (?k, t) . d3 k Qαβ (k, t)Qβα (?k, t) 13 (28) d3 x Qαβ (x, t)Qβα (x + r, t) , d3 x Qαβ (x, t)Qβα (x, t)
Theoretical predictions and analytical work have veri?ed that these correlation functions have a scaling form C (r, t) = F [r/L(t)], and S (k, t) = Ld (t)G [kL(t)] . Here, C (r, t) =
| r | =r
C (r, t) and S (k, t) =
| k| = k
S (k, t) are angular averages of the correlation functions
in real and Fourier space respectively. The length L(t) is extracted from the real-space correlation function using the implicit condition C (r = L(t), t) = 1/2. In Fig. (6) we con?rm that the real-space correlation function does indeed scale as expected. Our numerical data for the scaling function is in close agreement with an analytical calculation for the n = 2 model due to Bray and Puri . In Fig. (7) we show the corresponding scaling of the Fourier space correlation function. The wavenumber k is the root of the second moment of the S(k, t) de?ned by k
1 = L(t)2
k 2 S (k, t)/
S (k, t).
The inset shows the growth of the length scale as a function of time. Theoretically, this is expected to grow as a power L(t) ? tα . Our estimate for this exponent is α = 0.5 ± 0.005. Our results are consistent with both analytical predictions and an earlier numerical simulation. The Fourier space correlation function is expected to exhibit a short-wavelength scaling S (k, t) ? k ?4 known as a generalised Porod’s law . We see a clear range of wavenumbers where the Porod scaling is obtained. At very short wavelengths, corresponding to the size of the defect core, the Porod scaling breaks down. We see evidence for this as well, where the very highest wavenumbers in Fig. (7) show deviations from the Porod scaling. Our numerical results for spinodal decomposition, then, agree both qualitatively and quantitatively with theoretical results and previous numerical work . There are numerical results reported in the literature are con?icting [21, 22]. We indicate possible reasons for this discrepancy in the following section.
For the nematic-isotropic interface, elastic anisotropy (L2 = 0) induces biaxiality near the interface. This problem has been studied using various approximations by several authors. Sen and Sullivan  introduced a symmetry-adapted parametrisation of the order parameter which reduces the independent degrees of freedom to three. Using this parametrisation PopaNita and Sluckin [24, 25] obtained numerical solutions for the variation of S and T across 14
the interface. The parametrisation used in these calculations is only applicable to the cases of planar or homeotropic anchoring at the interface. For a general anchoring condition this parametrisation breaks down, as is obvious from the fact that the order parameter has, in general, ?ve independent components. With the present method, the same problem can be studied without any approximation, retaining all degrees of freedom of the orientation tensor. This allows us to obtain the variation of the order parameter in situations where the symmetry conditions of Sen and Sullivan do not apply. These include curved interfaces and the e?ect of higher order elasticity. With regard to the problem of nematic droplets a detailed study becomes possible, especially in three dimensions. Indeed the only numerical work we know of is the Monte Carlo simulation of Cuetos and Dijkstra . The advantage of the Landau-de Gennes approach for this problem is that the free energy of anchoring need not be postulated, (as is done when only the director degrees of freedom are retained and a Rapini-Papoular  free energy used), but is e?ectively present in the Landau-de Gennes free energy itself. This is the calculational strategy used by Prinsen and van der Schoot . Non-trivial director con?gurations have been proposed for di?erent strengths of bulk-to-surface coupling. Such predictions can be veri?ed cleanly with the present method. For nematic coarsening in three dimensions, simulations within the Landau-de Gennes framework have not been performed so far. The extended nature of topological defects, which are lines rather than points in three-dimensional nematics, makes the dynamics of the coarsening problem quite di?erent from two dimensions. The strength of the Landaude Gennes approach is particularly clear here. Having access to the strength of ordering makes it easy to locate the defect lines as tubes of zeros of the ordering strength. This is considerably easier than using a Burgers-like circuit integral of the director ?eld to locate defects as is done in the work of Zapotocky et al . The dynamics of defect lines has implications in cosmology, where the Kibble mechanism  predicts a scaling relation for the number density of defects. Experiments studying the Kibble mechanism in liquid crystals have been performed . However, we know of no numerical study of the Kibble mechanism along the lines proposed here for nematic liquid crystals. This provides further impetus for three-dimensional simulations. There are two alternatives to the MOL discretisation for tensor order parameter descriptions of the nematic phase. These are the cell dynamical scheme of Oono and Puri  as 15
implemented by Zapotocky et al  and Dutta and Roy , and the lattice Boltzmann method of Denniston et al . Our method di?ers in an important way from both these approaches, in that it provides a direct discretisation of the governing equations of motions. In cell-dynamical simulations, a local map is constructed whose ?xed points coincide with the extrema of the local part of the free energy. This local map is phenomenological and does not follow from a Landau expansion. The advantage of this method is that it produces sharp interfaces, ensuring a rapid scale separation between the microscopic interfacial length and the macroscopic domain size. In coarsening simulations, this reduces the non-universal o?set time at which the system enters the scaling regime. In comparison, the Landau-de Gennes free energy leads to a less sharp interface. Consequently, the structure of the core region of topological defects is di?erent in the two methods, being smaller in cell dynamics and larger in the present method. We see evidence of this in the dynamic structure factor, Fig. (7), where there is a violation of Porod scaling at the highest wavenumbers. As shown in Figs. (5) the core region of our defect spans several lattice spacings in each direction. Thus, as regards scale separation between what are microscopic and macroscopic lengths, the cell dynamical scheme fares better. On the other hand, when physics at the scale of the interface is important, as in the problem of the planar interface and nematic droplets, the MOL discretisation o?ers a clear advantage. Due to the special nature of the cell dynamical Laplacian, it is not clear how to generalise the cell dynamical method to include elastic anisotropy (corresponding to the gradient term with coe?cient L2 ) or higher order terms involving antisymmetric contractions. Due to the phenomenological nature of the cell dynamical map, it is di?cult to make direct contact with experiment, a procedure which is straightforward in the MOL when the non-dimensionalisation we outlined above is applied. In the lattice Boltzmann method, the governing parabolic equations are replaced by a hyperbolic superset  through the introduction of a distribution function for the tensor order parameter. As with cell dynamics, the Landau-de Gennes equations are not solved directly. Rather, lattice Boltzmann relies on a temporal scale separation which allows the hyperbolic equations to mimic parabolic behaviour. The lattice Boltzmann method needs to be carefully formulated  to avoid the presence of micro-currents and this procedure appears to be non-universal . We ?nd no evidence of spurious micro-currents in the MOL discretisation presented here. The lattice Boltzmann method includes the coupling of 16
order parameter and ?ow. It is entirely straightforward to extend the MOL discretisation to include order parameter advection. To include coupling to the ?uid momentum, we advocate a hybrid strategy, where the order parameter dynamics including advection is solved by a MOL discretisation, but the dynamics of the ?uid momentum including order parameter stresses is solved by the lattice Boltzmann method. In terms of algorithmic complexity, cell dynamics, lattice Boltzmann and the MOL with ?nite-di?erences appear to be matched, since both use fairly local information to calculate derivatives and involve explicit temporal updates. For N degres of freedom the algorithmic complexity of all three algorithms is O (N ). On the other hand, the storage requirements of the lattice Boltzmann formulation are larger by a factor of 6 to 9 in two-dimensions, and about 15 to 19 in three dimensions. On an Intel(R) Core(TM)2 Duo CPU with speed 2.66GHz, our code takes 0.33 seconds for one time step once a lattice size of 256 × 256.
This paper has proposed an e?cient numerical scheme based on the method of lines for the solving the Landau-de Gennes equations of nematodynamics. The numerical results obtained in previous sections are in excellent agreement with analytical results where available and consistent with previous numerical data. We expect that the method presented here will ?nd broad application in exploring the rich physics of the nematic phase of liquid crystals.
We thank the Indo-French Centre for the Promotion of Advanced Research and the DST, India for support.
 P. G. de Gennes and J. Prost, The Physics of Liquid Crystals (Clarendon Press, Oxford, 1993), 2nd ed.  P. Chaikin and T. Lubensky, Principles of Condensed Matter Physics (Cambridge University Press, 1995), 1st ed.  P. G. de Gennes, Molecular Crystals and Liquid Crystals 12, 193 (1971).
 O. A. Liskovets, J. Di?. Eqs. 1, 1308 (1965).  R. J. LeVeque, Finite Di?erence Methods for Ordinary and Partial Di?erential Equations, Steady State and Time Dependent Problems (SIAM, 2007).  A. Cuetos and M. Dijkstra, Phys. Rev. Lett. 98, 095701 (2007).  G. Rien¨ acker, M. Kr¨ oger and S. Hess, Physica A: Statistical Mechanics and its Applications 315, 537 (2002).  L. N. Trefethen, Spectral Methods in Matlab (SIAM, 2000).  R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems (Cambridge University Press, 2002).  M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables (Dover, 1964).  B. Fornberg, SIAM Review 40, 685 (1998).  J. D. Bernal and I. Fankuchen, J. Gen. Physiol. 25, 111 (1941).  C. Herring, Phys. Rev. 82, 87 (1951).  S. Chandrasekhar, Molecular Crystals and Liquid Crystals 2, 71 (1966).  E. G. Virga, Variational Theories for Liquid Crystals (Chapman and Hall, London, 1994).  P. Prinsen and P. van der Schoot, Phys. Rev. E 68, 021701 (2003).  N. Schopohl and T. J. Sluckin, Phys. Rev. Lett. 59, 2582 (1987).  A. J. Bray, Adv. Phys. 43, 357 (1994).  S. Puri and A. J. Bray, Phys. Rev. Lett. 67, 2670 (1991).  C. Denniston, E. Orlandini, and J. M. Yeomans, Phys. Rev. E 64, 021701 (2001).  M. Zapotocky, P. M. Goldbart, and N. Goldenfeld, Phys. Rev. E 51, 1216 (1995).  S. Dutta and S. K. Roy, Phys. Rev. E 71, 026119 (2005).  A. K. Sen and D. E. Sullivan, Phys. Rev. A 35, 1391 (1987).  V. Popa-Nita and T. J. Sluckin, J. Phys. II (France) 6, 873 (1996).  V. Popa-Nita, T. J. Sluckin, and A. A. Wheeler, J. Phys. II (France) 7, 1225 (1997).  A. Rapini and M. J. Papoular, J. Phys. (Paris) Colloq. 30, C4 (1969).  T. W. B. Kibble, J. Phys. A 9, 1387 (1976).  I. Chuang, R. Durrer, N. Turok, and B. Yurke, Science 251, 1336 (1991).  Y. Oono and S. Puri, Phys. Rev. Lett. 58, 836 (1987).  C. Denniston, E. Orlandini, and J. M. Yeomans, Phys. Rev. E 63, 056702 (2001).
superheating spinodal line I?N transition line supercooling spinodal line
UN?BN transition line biaxial nematic phase isotropic phase
uniaxial nematic phase
?6 ?4 ?2 0 2 4 6
FIG. 1: (Color online) The mean ?eld phase diagram obtained from the Landau-de Gennes free energy (7). The phase boundaries are given by solutions of the following algebraic equations : 4B 4 E + B 2 C (3C 2 ? 144AE ′ ) = A(81C 4 ? 864AC 2 E ′ + 2304A2 E ′ ) for the isotropic to uniaxial nematic transition; 25A3 E ′ = ?18B 2 C 3 for the uniaxial to biaxial nematic transition; A=0 for the supercooling spinodal; 9B 4 E ′ + 8B 2 (C 3 ? 36ACE ′ ) = 192A(C 2 ? 4AE ′ )2 for the superheating spinodal.  S. Succi, Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Oxford U.P., 2001).  S. V. Lishchuk, C. M. Care, and I. Halliday, Phys. Rev. E 67, 036701 (2003).  S. V. Lishchuk, C. M. Care, and I. Halliday, J. Phys.:Condens. Matter 16, S1931 (2004).  A. J. Bray and S. Puri, Phys. Rev. Lett. 67, 2670 (1991).
0.045 0.04 0.035 10 8 6 4 2 0 0.02 0.015 0.01 0.005 0 110 0 0.05
Width Fit L1
S Fit T
Uniaxial nematic phase
125 130 135 140 145
FIG. 2: (Color online) Variation of the degree of alignment (S) and biaxiality (T) across the nematic - isotropic interface with planar anchoring. Symbols represent numerical data, solid curves are the de Gennes ansatz (26). The numerical parameters are A = 0.0035, B = -0.5, C = 2.67, E ′ = 0, L1 = 0.01, L2 = 0, Γ = 1/20, grid size 8 × 512 and nematic strip width = 256. The inset shows the variation of the interfacial width with the elastic constant L1 . The expected quadratic variation is accurately reproduced.
(a) t = 0
(b) t = 1000
(c) t = 2000
(d) t = 3000
(e) t = 0
(f) t = 300
(g) t = 600
(h) t = 900
(i) t = 0
(j) t = 500
(k) t = 1000
(l) t = 1500
FIG. 3: (Color online) Time evolution of the degree of alignment in the relaxation of a nematic droplet. The initially circular droplet develops an anisotropy and becomes elliptical, the twodimensional analogue of the ellipsoidal droplets (tactoids) seen in three dimensions. Figs. 3(a), 3(e) and 3(i) shows the circular droplet at time t = 0. The time evolution with L2 = 0 of the circular droplet are shown in Figs. 3(b), 3(c) and 3(d). The elliptical conformations with L2 > 0 are shown in Figs. 3(f), 3(g) and 3(h). These have aspect ratio 1.1191, 1.3046 and 1.5037 respectively. With L2 < 0 conformations are shown in Figs. 3(j), 3(k) and 3(l). These have aspect ratio 1.0122, 1.0718 and 1.1316 respectively. The numerical parameters are A = 0.001, B = -0.5, C = 2.67, E ′ = 0, L1 = 0.0236, Γ = 1.0, grid size 128 × 128 and the droplet radius = 20. The elastic constant L2 = 10L1 for the L2 > 0 and L2 = ?L1 for L2 < 0.
FIG. 4: (Color online) Schlieren textures in a coarsening nematic. Topological defects of both integer charge (top right in (c) and (d)) and half - integer charge (throughout) are observed in the simulations. (a) shows the random con?guration at time t = 0. (b), (c), (d), (e) and (f) show the dynamics at di?erent times t = 103 , 5 × 103 , 2 × 104 , 6 × 104 and 105 respectively. The parameters chosen are A = -0.1, B = -0.5, C = 2.67, E ′ = 0.0, L1 = 1.0, L2 = 0 and Γ = 1.0, grid size 256 × 256 for 105 time steps.
45 40 35 30 0.14 25 20 15 10 0.08 5 5 10 15 20 25 30 35 40 45 0.12 0.18
0.2 0.18 0.16 S T
0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0 5 10 15 20 25 30 35
FIG. 5: (Color online) The schlieren texture around a half integer charge defect is shown in (a), (b) shows the density plot of the uniaxial order parameter, while (c) exhibits the variation of the uniaxial and biaxial order parameter along 23a line passing through the defect core. Note the presence of strong biaxiality within the defect core as seen previously . The sharpness of the
0.8 0.7 0.6
0.6 0.4 0.2
0 0 10 20 30 t=600 t=700 t=800 t=900 BP function 40
0.5 0.4 0.3 0.2 0.1 0 0 0.5 1
t=100 t=200 t=300 t=400 t=500
FIG. 6: (Color online) Data collapse of the direct correlation function C (r ) with scaled distance r/L(t) for di?erent times. The symbol ? depicts the Bray - Puri function  for the O(2) model : fBP (x) = B 2 (0.5, 1.5)F [0.5, 0.5, 2; exp(?x2 )]exp(?x2 /2)/π . The inset shows the unscaled correlation function at di?erent times. Numerical parameters chosen were A = ?0.1, B = ?0.5; C = 2.67; E ′ = 0, L1 = 1.0, L2 = 0, Γ = 1/20 on a 256 × 256 grid
?1.5 ?2 ?2.5 t=300 t=400 t=500 t=600 t=700 t=800 10 ?3.5
ln[S(k, t) <k>2]
?4 ?4.5 ?5 ?5.5 ?2.5
10 1 10
?2 ?1.5 ?1 ?0.5 0 0.5 1
FIG. 7: (Color online) Data collapse of the structure function S (k, t) at di?erent times. The dashdot line has a slope of ?4 indicating the validity of generalised Porod’s law for O(n) systems. The departure from Porod’s law at high k is due to the ?nite core size of the defects as discussed in the text. The inset shows the time dependence of the correlation length L(t). The length grows as a power law with an exponent of 0.5. The maximum value of the correlation length is approximately 1/4-th the system size, ensuring the absence of ?nite-size artefacts. Numerical parameters chosen are A = ?0.025, B = ?0.5, C = 2.67, E ′ = 0, L1 = 0.1, L2 = 0, Γ = 1/20 on a 256 × 256 grid.