Lattice-Boltzmann simulation of free nematic-isotropic interfaces

We use a hybrid method of lattice Boltzmann and finite differences to simulate flat and curved interfaces between the nematic and isotropic phases of a liquid crystal described by the Landau-de Gennes theory. For the flat interface, we measure the interfacial velocity at different temperatures around the coexistence. We show that the interface is completely static at the coexistence temperature and that the profile width is in line with the theoretical predictions. The interface is stable in a range of temperatures around coexistence and disappears when one of the two phases becomes mechanically unstable. We stabilize circular nematic domains by a shift in temperature, related to the Laplace pressure, and estimate the spurious velocities of these lattice Boltzmann simulations.


Introduction
Liquid crystals are systems that flow as liquids but where the particles are partially organized in crystal-like ways [1,2]. They are formed by elongated particles the orientations of which may point in a preferential direction. At low temperatures, these liquid crystals are in the nematic phase (orientationally ordered) while at high temperatures, they are in the isotropic phase, where the orientation of the particles is random (orientationally disordered). At the nematic-isotropic temperature, the two phases coexist and a stable interface may be observed. The standard theory that describes this static interface is based on a tensor order parameter expansion, with gradient terms that account for the order parameter variations, proposed by de Gennes [1].
To simulate the dynamics of liquid crystals, it is common to solve numerically the continuum equations of Beris-Edwards and Navier-Stokes [3][4][5]. A hybrid method of lattice Boltzmann [6] and finite differences [7] is used. Special care has to be taken in order to simulate interfaces due to spurious currents that arise in these numerical methods, which can lead to interfacial motion or changes in the director field [6,8].
We analyze the applicability of the method of Ref. [4] to simulate free interfaces, which are very sensitive to these spurious numerical effects. We show that it is possible to simulate static flat interfaces at coexistence, with the profiles predicted by the Landau-de Gennes theory. We estimate the magnitude of the spurious velocities at nearly static circular interfaces stabilized by temperature shifts, related to the Laplace pressure arising from the curvature.

Equations of motion
To simulate the liquid crystal dynamics, we solve the Beris-Edwards equation coupled to the Navier-Stokes equation. The first describes the evolution of the order parameter Q αβ = S (n α n β − δ αβ /3) Here, ξ stands for the alignment parameter, Γ is the rotational diffusive constant and the corotational term is given by: where A 0 is a constant, γ is a temperature related parameter and K is the elastic constant. Thus, the molecular field, H αβ = −δF /δQ αβ + (δ αβ /3) Tr(δF /δQ γǫ ), is The continuity and Navier-Stokes equations describe the evolution of the density and velocity of the fluid where, the stress-tensor is P 0 = ρc 2 s is the hydrostatic pressure. Using Eq. (3), we obtain δF /δ(∂ β Q γν ) = K∂ β Q γν , which can be replaced in Eq 6.

Hybrid method
Our numerical scheme solves Eq. (1) using finite differences and Eq. (5) using lattice Boltzmann. Both methods are solved using the same grid (spatial discretization) and the same time step.

Lattice Boltzmann
The lattice Boltzmann method is a numerical technique which solves the Boltzmann equation and recovers the Navier-Stokes equation in the macroscopic limit. The space is discretized in a regular grid and the velocity space is discretized according to the lattice (Gaussian quadrature). Here, we use the D3Q19 lattice, which has 19 velocity vectors c i isotropically distributed in three dimensions. The discrete version of the Boltzmann equation is implemented in the method where f i is the distribution function corresponding to the c i vector. The simplest collision operator is the Bhatnagar-Gross-Krook (BGK) one, which assumes that the fullequilibrium distribution f i relaxes to the equilibrium f eq i with a characteristic relaxation time τ: (∂ f i /∂t) coll = −(∆t/τ)( f − f eq ). The multi-relaxation time operator, used in this work, is a generalization of the BGK: ( where the matrix M transforms from the distribution function space to the hydrodynamic moments space and R is the relaxation matrix. It assumes that the hydrodynamic moments may relax with different relaxation times. This approach is known to improve the accuracy and stability in many problems [6]. The equilibrium distribution is the second order expansion in Hermite polynomials of the Maxwell-Boltzmann distribution [9]: where c s is the speed of sound in the given lattice. For the D3Q19, it is c s = 1/ √ 3. The density and velocity fields are calculated as follows: The source term, for the BGK operator, reads: For the MRT operator, this source term must be relaxed in the moments space as the distribution function (see Refs. [4,6] for more details of the MRT method). The force is calculated using the stress tensor of Eq. (6): F α = ∂ β (Π αβ + P 0 δ αβ ).

Finite differences
The Eq. (1) is solved explicitly using a predictor-corrector finite difference method. All the differences are second order accurate. For instance, the first derivative of a vector V in the x-direction is: ∂ x V(x) = [V(x+∆x)−V(x−∆x)]/(2∆x)+O(∆x 2 ). To calculate the time evolution of Q αβ , we follow the following steps: 1) Calculate the time derivative, (∂Q αβ /∂t) ′ , using Eq. (1) with the current Q αβ ; 2) Calculate the predictor: Q P αβ = Q αβ + ∆t(∂Q αβ /∂t) ′ ; 3) Calculate the time derivative using the predictor, ( ; 5) Calculate the corrector: Q αβ (t + ∆t) = Q αβ (t) + ∆t ∂Q αβ ∂t ; 6) Return to step 1. Notice that the velocity used in Eq. (1) is the actual one (and not the lattice Boltzmann one: ρu LB = i c i f i ), given by Eq. (9), which is second order accurate.

Free interface
In this section, we show the results from numerical simulations of unconfined nematicisotropic interfaces. First, we consider a flat interface in two dimensions and then we simulate a circular interface. In the simulations, we apply periodic boundary conditions in both directions. The initial velocity is set to zero and the density to ρ = 1 everywhere. Other parameters are: τ = 1.5, K = 0.04, ξ = 0.7, A 0 = 0.1 and Γ = 0.34. Our results are given in lattice units: the distance between nodes is ∆x = 1 and the time step is ∆t = 1.

Flat interface
In order to apply periodic conditions, we initialize our system as follows: the nematic phase, with directors in the vertical direction, is in the center of the domain with isotropic phase elsewhere. Thus, we simulate two interfaces, but we focus our analysis on the right interface. At the coexistence temperature (γ NI = 2.7) the interface is completely static and there are no spurious velocities. Eq. (3) allows an interface solution of the form [10] where S n is the nematic order parameter and λ = 27K/A 0 γ is the correlation length. At γ = 2.7, λ = 2 and S n = 1/3. Fig. 1 D shows that the simulation reproduces accurately the interfacial profile. Next, we change the temperature slightly away from coexistence. In this case, the interface moves with constant velocity, which can be to the left (negative) or to the right (positive). Figs. 1 A and C, show the interfaces at two different temperatures after an interval of time. For γ < γ NI the domain of the nematic phase shrinks and for γ > γ NI the nematic domain expands. In the inset of Fig. 1D, we plot the interface velocity at different temperatures. It is approximately linear close to γ NI but is deviates from the linear behavior as the shift from coexistence increases. Above a certain shift from γ NI one of the two phases becomes unstable and the entire domain becomes nematic or isotropic in the absence of interfacial propagation. Below γ < = 2.65 the system becomes isotropic and above γ > = 3.1 it becomes nematic.

Circular interface and spurious velocities
Here we estimate the spurious velocities generated at a circular interface. The spurious velocities arise at curved interfaces as a result of the discretization of the velocity space. However, it is difficult to distinguish the spurious velocities from the physical ones in this multiphase model, because curved interfaces are usually unstable. Fig. 2E shows the interface position at three different temperatures. In order to identify the spurious velocities, we iteratively search for the temperature where the interface velociy is of the order of the spurious ones. We found that for γ = 2.707681 the interface velocity at t = 35000 is approximately u i ≈ 1.89 × 10 −7 . Fig. 2A and B show the oder parameter and velocity field of this almost static interface. From the vorticity field, we notice that the there are vortices between the directions of the velocity vectors of the D3Q19 lattice, which is an indication that these are spurious currents. In Fig. 2 D we plot the velocity field for the same system rotated by 20 • . One notes that the intensity of the vortices changed due to the physical component of these velocities (the small interface velocity), but the positions of the vortices are the same confirming that they are mostly due to the velocity discretization. From this analysis, we estimate that the spurious velocities are around 10 −6 in lattice units for this setup. If necessary, one could reduce the spurious velocities by increasing the isotropy of the lattice [11] or using more advanced collision operators [12].

Conclusions
We examined the applicability of the hybrid method of lattice Boltzmann and finite differences for the simulation of nematic-isotropic interfaces. The shape and width of the interface is correctly reproduced and the free interface is static at the coexistence temperature. At temperatures slightly different from the coexistence, the interface moves with constant velocity that depends on the temperature. The interface exists only for a certain range of temperatures: for temperatures much lower or higher than the coexistence, one of the two phases becomes unstable and the interface disappears without propagation. The circular interfaces are not static at the coexistence temperature. Thus, we stabilize them by a temperature shift in a way that the interface velocity is comparable to the spurious velocities. This allows us to visualize the profile of theses spurious currents and to estimate their magnitude.