ABSTRACT

We investigate the flow induced by annihilation of a pair of defects in liquid crystals using the Doi theory with the Marrucci-Greco potential, in which the orientation state is described with the orientational distribution function. We have numerically studied both the transient behaviors of two defects with different structures and their velocity field, and estimated the magnitude of the induced velocity. A defect with positive strength moves faster than one with negative strength. The long-range order of the molecular orientation field has a large effect on the annihilation time, and the annihilation time is reduced by increasing the long-range order. We find that flows are induced during the annihilation of a pair of defects and that several vortices are generated in the vicinity of the defects. The maximum velocity is predicted to develop spatially between the two defects just after their annihilation in time. In our simulation, the maximum induced velocity reaches an order of 10 μm/s. The induced velocity increases with increasing long range-order and nematic potential strength.

Keywords:

Complex Fluid, Liquid Crystal, Defect, Annihilation, Generation of Flow

1. Introduction

At present, liquid crystals are widely used as displays (low-molar-mass liquid crystals) and high-strength engineering plastics (liquid crystalline polymers) by applying their anisotropy of material constants. On the other hand, because a liquid crystal has both solid and liquid properties, it has a possibly much wider range of applications, as do solid and liquid. Therefore, studies on new applications of liquid crystals are needed.

When a nematic liquid crystal is observed under a polarizing microscope, a region in which molecules orient parallel or perpendicular to the direction of polarization appears as a dark field, a region in which molecules orient ±45˚ to the direction of polarization appears as a bright field, and a region in which molecules orient in the other direction appears as a gray field. Figure 1(a) is a photograph of a nematic state right after the phase transition from an isotropic state of N-(p-methoxybenzylidene)-p’-butyl-aniline (MBBA), a typical low-molar-mass nematic liquid crystal. We observe the Schlieren texture peculiar to liquid crystals. There are some points in Figure 1(a) where black lines intersect with each other. These points are called defects [1] [2] , around which spatial distortions of molecular orientation fields are generated. Figure 1(b) shows the representative defect orientation states. At the defect cores, the direction of orientation is discontinuous, and the orientation state is singular.

The defects often generated in manufacturing liquid crystal products cause degradation of the productivity and performance of the products. However, it is experimentally known that a pair of defects such as those in Figure 1(b), which have different molecular orientations, attract and finally annihilate each other [2] . This defect annihilation changes the orientation direction, so that a flow is expected to be generated [3] - [8] . Applying thermal energy induces defect generation because the orientation order at defects is lower. Therefore, generating defects artificially at arbitrary locations by applying thermal energy would inevitably generate flows there, leading to the development of a new type of microactuator that converts thermal energy into kinetic energy.

A theory that is able to describe the molecular orientation state satisfactorily should include the effects of long-range and short-range molecular order as well as the effect of flow. Two theories, the Leslie-Ericksen theory [9] [10] and the Doi theory [11] [12] , are well-known and have been widely used for predicting the dynamics of liquid crystals. From the aforementioned viewpoint, however, the Leslie-Ericksen theory, which does not include the short-range order effect, and the Doi theory, which does not account for the long-range order effect, are out as candidates. Marrucci and Greco [13] have expanded the Maier-Saupe potential [14] used in the Doi theory to include both the short-range and long-range

order effects and applied it to the original Doi theory [15] . Feng et al. [16] have proposed an expression for the stress tensor and to complete the theory.

Our objectives are to study the flow induced by the annihilation of a pair of defects and to estimate its magnitude using the aforementioned theory as a constitutive equation. A number of studies have simulated the liquid crystal flow during the annihilation of paired defects [17] [18] [19] [20] . Experiments for the moving speed of defects under electric fields have also been performed [21] [22] . However, they do not aim to develop actuators but focus on the movement of defects rather than the liquid crystal flow.

To calculate the orientational distribution functions (ODFs) at many points in a physical space, we have to account for both computational accuracy and time. In this paper, we: 1) approximate the ODF with a series of spherical harmonic functions, 2) study the minimum number of expanded terms required to simulate the orientation state properly, and 3) finally estimate the magnitude of the induced flow during the annihilation of a pair of defects to obtain useful data that can contribute to developing new actuators.

2. Computations

2.1. Governing Equations

The evolution equation for the ODF f is written as

$\frac{\partial f}{\partial t}=\bar{D}{\nabla }_u\cdot \left({\nabla }_{u}f+\frac{f{\nabla }_{u}V}{kT}\right)-{\nabla }_u\cdot \left\{f\left(\kappa \cdot u-\kappa :uuu\right)\right\}$ . (1)

Here, t is the time, k the Boltzmann constant, T the absolute temperature, u the unit vector parallel to a liquid crystal molecule, ${\nabla }_u$ the differential operator on a unit sphere, and κ the velocity gradient tensor. $\bar{D}$ and V are the rotary diffusivity and Marrucci-Greco nematic potential [13] expressed as

$\bar{D}=D{\left(1-\frac{3}{2}S:S\right)}^{-2}$ , (2)

and

$V=-\frac{3}{2}UkT\left(S+\frac{l_i^2}{24}{\nabla }^{2}S\right):uu$ (3)

where D is the rotary diffusivity in an isotropic state, $\nabla$ the differential operator in physical space, and S the order parameter tensor defined as

$S={\displaystyle {\int }_{\left|u\right|=1}\left(uu-\frac{\delta }{3}\right)f\text{d}\Omega }\equiv \langle uu-\frac{\delta }{3}\rangle$ . (4)

δ is the unit tensor.

The conservation equations for the isothermal slow flow of liquid crystals are the continuity and linear momentum equations

$\nabla \cdot v=0$ , (5)

and

$\rho \frac{\partial v}{\partial t}=-\nabla p+\nabla \cdot \tau$ , (6)

where v is the velocity vector, ρ is the fluid density, and p is the pressure. τ is the extra stress tensor derived by Feng et al. [16] expressed as

$\begin{array}{l}\tau =3ckT[A-U\left(A\cdot A-A:Q\right)-\frac{U{l}_i^2}{24}\{A\cdot {\nabla }^{2}A-Q:{\nabla }^{2}A\underset{}{\overset{\stackrel{}{}}{{\displaystyle }}}\\ +\frac{\left(\nabla A\right):{\left(\nabla A\right)}^{{\rm T}}}{4}-\frac{\nabla \nabla A:A}{4}\}]+\frac{c}{2}\zeta \kappa :Q\end{array}$ (7)

Here, c is the number density of liquid crystal molecules, U the dimensionless nematic potential intensity, l_i a parameter indicating the long-range order effect of molecules, ζ the drag coefficient of rotary molecules. A and Q are the second and fourth moments of the ODF f, respectively, expressed as

$A={\displaystyle {\int }_{\left|u\right|=1}uu}f\text{d}\Omega \equiv \langle uu\rangle$ , (8)

and

$Q={\displaystyle {\int }_{\left|u\right|=1}uuuu}f\text{d}\Omega \equiv \langle uuuu\rangle$ . (9)

2.2. Computational Procedure and Boundary Conditions

Let us consider a two-dimensional square area with a side length of H, shown in Figure 2, and put a pair of defects with different orientational states on two points whose coordinates are P (0.3H, 0.5H) and Q (0.7H, 0.5H). For the coordinate system in Figure 2, the components of the velocity gradient tensor κ and order parameter tensor S are

$\kappa =\left|\begin{array}{ccc}\frac{\partial u}{\partial x}& 0& \frac{\partial u}{\partial z}\\ 0& 0& 0\\ \frac{\partial w}{\partial x}& 0& \frac{\partial w}{\partial z}\end{array}\right|$ , (10)

and

$S=\left|\begin{array}{ccc}{S}_{xx}& 0& S_{xz}\\ 0& S_{yy}& 0\\ S_{zx}& 0& S_{zz}\end{array}\right|$ , (11)

where u and w are the x and z components of the velocity vector v.

For computation of the orientation field, we substitute Equation (10) and Equation (11) into Equation (1) to obtain the ODF f. In this study, we approximate f with a finite series of spherical harmonic functions Y_lm(u) [23] [24] [25] .

$f\left(t,u,x,z\right)={\displaystyle \underset{l=0}{\overset{l_{\max }}{\sum }}{\displaystyle \underset{m=-l}{\overset{l}{\sum }}{C}_{lm}\left(t,x,z\right)Y_{lm}\left(u\right)}}$ . (12)

Here, C_lm(t, x, z) are coefficients, and l_max is the maximum of the azimuthal quantum number on which the number of terms of the series solution depends.

Since the head and tail of a liquid crystal molecule are indistinguishable, we have $f\left(t,u\right)=f\left(t,-u\right)$ . From the parity of the spherical harmonic functions, that is,

$Y_{lm}\left(-u\right)={\left(-1\right)}^l{Y}_{lm}\left(u\right)$ , (13)

the expression (?1)^l =1 is obtained, which restricts l to even values. We have non-dimensionalized Equation (1) with 1/D, multiplied the resulting equation by the complex conjugate of spherical harmonic functions, $Y_{lm}^{*}\left(={\left(–1\right)}^m{Y}_{lm}\right)$ , and integrated over the unit sphere to get the ordinary differential equations with respect to C_lm (see Appendix 1). The orthonormality of the spherical harmonic functions [26]

${\displaystyle {\int }_{\left|u\right|=1}{Y}_{lm}{Y}_{l^{\prime }{m}^{\prime }}^{*}\text{d}\Omega }={\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}$ (14)

has been used.

To compute the velocity field, we eliminate p from Equation (6) to obtain the vorticity transport equation

$\rho \frac{\partial \omega }{\partial t}=\frac{{\partial }^2{N}_1}{\partial x\partial z}-\frac{{\partial }^2{\tau }_{xz}}{\partial x^2}+\frac{{\partial }^2{\tau }_{zx}}{\partial z^2}$ , (15)

where $\omega \left(=\partial u/\partial z-\partial w/\partial x\right)$ is the vorticity, and N₁ (=τ_xx − τ_zz) is the first normal stress difference. Stresses N₁, τ_xz, and τ_zx are explained in the Appendix 2.

We have used the finite-difference scheme to discretize the equations and the implicit Euler method for time integration. Equation (15) is solved using an iterative procedure with the convergence criterion that the average relative error at each node is less than 10⁻⁵. We use periodic boundary conditions.

Values of the physical quantities are the fluid density ρ = 10³ kg/m³, the absolute temperature T = 320 K, the rotary diffusivity in an isotropic state D = 5.2 × 10³ s⁻¹, the number density of molecules c = 2.25 × 10²⁴ m⁻³, the drag coefficient of rotary molecules ζ = 8.89 × 10⁻²⁴ kgm²/s, and the side length of the computational domain H = 1 μm. The computational parameters we select are the nematic potential intensity U = 5.0, 5.5, and 6.0, and the long-range order effect $l_i^{\ast }$ (=l_i/H) = 0.02 − 0.1 with a step of 0.01. The choice of l_max is important, and we have determined it by accounting for both the computational accuracy and load. We report the details in the following chapter.

The mesh size is set to $\Delta x^{*}\left(=\Delta x/H\right)=\Delta z^{*}\left(=\Delta z/H\right)={10}^{-2}$ , and the time step Δt^* (=ΔtD) is varied depending on $l_i^{\ast }$ ; for example, Δt^* = 5 × 10⁻⁴ when $l_i^{\ast }=0.1$ .

2.3. Initial Values

The initial velocity vector v is 0. The initial values of C_lm are obtained as follows: We multiply Equation (12) by $Y_{lm}^{*}$ , integrate over the unit sphere, and use Equation (14) to get

$C_{lm}\left(0,x,z\right)={\displaystyle {\int }_{\left|u\right|=1}f\left(0,u,x,z\right)Y_{lm}^{*}\left(u\right)\text{d}\Omega }$ . (16)

Assuming that f is in an equilibrium state (no flow) at t^* = 0, f is expressed as [12]

$f\left(0,u,x,z\right)=\exp \left(-\frac{V}{kT}\right)/{\displaystyle {\int }_{\left|u\right|=1}\exp \left(-\frac{V}{kT}\right)\text{d}\Omega }$ . (17)

We use the denominator to normalize f in Equation (17). Since f has uniaxial symmetry in the absence of both flow and external field, the order parameter tensor, Equation (4), is rewritten as

$S=S\left(nn-\frac{\delta }{3}\right)$ , (18)

where n is a unit vector describing the average local molecular orientation, called the director, and S is the scalar order parameter ranging from 0 in a random orientation state to 1 in a perfect alignment, defined as

$S=\sqrt{\frac{3}{2}S:S}$ . (19)

The symbol “:” means the double dot product of two tensors. Substituting Equation (3) and Equation (18) into Equation (17), expanding exponential terms into a power series, and expressing the power by the spherical harmonic functions give

$C_{lm}\left(0,x,z\right)={\displaystyle \underset{p=0}{\overset{\infty }{\sum }}\frac{{\left(\frac{3}{2}US\right)}^p\left(2p\right)!}{p!\left(2p-l\right)!!\left(2p+l+1\right)!!}{Y}_{lm}^{*}\left(n\right)}/{\displaystyle \underset{q=0}{\overset{\infty }{\sum }}\frac{{\left(\frac{3}{2}US\right)}^q}{\left(2q+1\right)q!}}$ . (20)

where $\left(2n\right)!!=\left(2n\right)\left(2n-2\right)\cdots 4\cdot 2$ , $\left(2n+1\right)!!=\left(2n+1\right)\left(2n-1\right)\cdots 3\cdot 1$ , 0!! = 1 (n is 0 or even number). In Equation (20), n and S are functions of x and z. We have determined both values (see Appendix 3).

3. Results and Discussion

3.1. Determination of l_max

Before computing the velocity and orientation fields, we must determine the value of l_max. Let us consider a liquid crystal to which simple shear flow is applied. Setting l_i = 0, Equation (3) reduces to

$V=-\frac{3}{2}UkTS:uu$ , (21)

which is the Maier-Saupe nematic potential [14] . The coordinate system is shown in Figure 3. The flow is in the x direction, and the velocity gradient is in the z direction. The polar and azimuthal angles of the unit vector u are θ and ϕ, respectively. For simple shear flow, the velocity gradient tensor is expressed as

$\kappa =\left|\begin{array}{ccc}0& 0& \stackrel{\dot{}}{\gamma }\\ 0& 0& 0\\ 0& 0& 0\end{array}\right|$ , (22)

where $\stackrel{\dot{}}{\gamma }$ is the shear rate. Equation (11), Equation (21), Equation (22) are substituted into Equation (1) to obtain the ODF f. Equation (20) can be used as an initial value of C_lm, but it is a function only of time t. At t^* = 0 the director n is assumed to orient in the x direction, so that we set θ_m = π/2 and ϕ_m = 0, where θ_m and ϕ_m are the polar and azimuthal angles of the director, respectively.

The azimuthal angle of the director ϕ_m is always 0 because the director is in the x-z plane in Figure 3. The polar angle θ_m is obtained from S as follows:

$\tan 2{\theta }_m=\frac{2S_{xz}}{S_{zz}-S_{xx}}$ . (23)

Figure 4 shows the transient behaviors of θ_m and S at l_max = 4 − 10 (no convergent solution was obtained at l_max = 2) when the potential intensity U = 5 and the dimensionless shear rate ${\stackrel{\dot{}}{\gamma }}^{*}\left(=\stackrel{\dot{}}{\gamma }/D\right)=1$ and 40. In Figure 4(a), the director rotates and the order parameter oscillates about the equilibrium value S = 0.615. The behaviors at l_max = 8 are almost the same as those at l_max = 10, and the behaviors at l_max = 6 have smaller periods compared with those at l_max = 8 and 10. At l_max = 4, steady-state instead of periodic behaviors are obtained. In Figure 4(b), the director and the order parameter reach their steady state values in a short time at l_max = 6, 8, and 10, whereas the predictions at l_max = 4 are even qualitatively different from the results at other l_max. Thus, l_max = 4 is not available, l_max = 6 is qualitatively acceptable, and l_max ≥ 8 is quantitatively

satisfactory. Expansion terms required to approximate the ODF depend on the potential intensity and shear rate. When U and ${\stackrel{\dot{}}{\gamma }}^{*}$ are large, more terms are necessary because the ODF becomes steep and non-uniaxial.

3.2. Orientation Fields

Since the flows induced by annihilation of a pair of defects are not large, and the selected potential intensity U in the present study is not high (5, 5.5, 6), we have used l_max = 6 for the following computations.

Figure 5 shows a sequence of predictions of transmitted light intensity observed under crossed nicols for $l_i^{\ast }=0.\text{1}$ . The transmitted light intensity I is evaluated using

$I=I_0{\sin }^2\left(2{\theta }_m\right)$ . (24)

Here, I₀ is the incident light intensity, and θ_m is the polar angle of the director obtained by Equation (23). A pair of defect cores attracts each other with time, and brushes connecting the defects become short. Finally the defects are annihilated at t^* = 3.85. Since the initial orientation profile is symmetric with respect to z^* = 0.5, the defect cores move along this line. After finishing the annihilation, the orientation state does not become homogeneous immediately, but slightly gray areas are discernible at t^* = 5. A completely homogeneous orientation state is achieved at t^* = 6.05.

Figure 6 shows the profiles of S for the same parameters and times as

Figure 5. Since the order parameters near the defect cores are low, the profile of S looks like an inverted cone. The two inverted cones gradually approach and coalesce laterally with each other as time advances. The two sharp ends merge into a single end at t^* = 3.85, which coincides with the time the defects

annihilate in Figure 5. Just after the annihilation, we can see slight inhomogeneity in S.

Figure 7(a) plots the locations of a pair of defect cores before annihilation for several $l_i^{\ast }$ . The defect cores approach at a constant speed, but they slow down just before the annihilation. The defect core with the strength s = +1/2 (at point P in Figure 2) moves slightly faster than that with the strength s = −1/2 (at point Q in Figure 2). Thus, the annihilation takes place at a position slightly over x^* = 0.5.

Figure 7(b) shows the relationship between the time $t_a^{\ast }$ required for the

annihilation process and the long-range order parameter $l_i^{\ast }$ . A large $l_i^{\ast }$ gives a short annihilation time because the attractive force between defects becomes strong for large $l_i^{\ast }$ . For example, $t_a^{\ast }=190$ at $l_i^{\ast }=0.02$ , and $t_a^{\ast }=3.85$ at $l_i^{\ast }=0.\text{1}$ (refer to Figure 5 and Figure 6); that is, increasing $l_i^{\ast }$ by a factor of 5 decreases $t_a^{\ast }$ by a factor of 50. Therefore, $l_i^{\ast }$ has a large effect on annihilation.

3.3. Velocity Fields

Figure 8 shows the velocity distributions at t^* = 1, 3.85, and 5 for $l_i^{\ast }=0.\text{1}$ as well as the locations of defect cores. An arrow at the bottom of each figure is the reference velocity vector. It is confirmed that flows are induced during annihilation of a pair of defects. Since the molecular orientation state is symmetric with respect to z^* = 0.5, the velocity distribution is also symmetric. Flows to the right are observed in the vicinity of z^* = 0.5. As a result, a counterclockwise vortex in the upper half region and a clockwise one in the lower half region are generated. Four vortices (two in each half region) are generated at t^* = 1 and 3.85. By comparing the velocity vectors near the left and right defect cores at t^* = 1, we find that the velocity vector and vortex on the left side are larger, and that the velocity distributions depend on the defect structure. Both right and left vortices rotate in the same direction at t^* = 1, so that they coalesce into a larger vortex at t^* = 3.85. Since the induced velocity vector is maximum spatially on the line connecting the defect cores (z^* = 0.5), it is efficient to use the flow between defects when using the flow induced during their annihilation. We have checked the velocity profiles at times different from those of Figure 8 to verify that the velocity vector has its maximum value when annihilation finishes (t^* = 3.85).

We define v_max as the maximum velocity value in the computational region at the time when annihilation finishes. Figure 9(a) plots the relationship between v_max and $l_i^{\ast }$ . When $l_i^{\ast }$ increases, the induced velocity also increases. However, the effect of $l_i^{\ast }$ on the velocity is not large compared with that on the annihilation

time $t_a^{\ast }$ ; the relationship between v_max and $l_i^{\ast }$ is almost linear. Figure 9(b) plots v_max against the nematic potential intensity U for $l_i^{\ast }=0.\text{1}$ . The induced velocity increases with increasing U. We explain this result as follows: Since annihilation stems mainly from the spatial gradient of S, an increase in U corresponds to an increase in S, resulting in a steep spatial gradient for S. Thus, a liquid crystal with higher liquid crystallinity (a low-temperature liquid crystal) may generate a faster flow.

4. Conclusion

In this study we have predicted flows induced by the annihilation of a pair of liquid crystal defects using the Doi theory with the Marrucci-Greco potential and the constitutive equation of Feng et al. The long-range order effect on the time required for the annihilation process of a pair of defects is remarkable; when the long-range order is large, the annihilation time becomes short. We have shown that a flow is induced by the annihilation and that several vortices are generated in the vicinity of the defects. The maximum flow is obtained on the line connecting the two defect cores in space, and at the time, the annihilation is just finished. The maximum value of the induced velocity is on the order of 10 μm/s in our study. The induced velocity becomes large when the long-range order and nematic potential strength are high.

Acknowledgements

This work was partially supported by the Japan Society for the Promotion of Science KAKENHI (Grant No. 25289035).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Chono, S. and Tsuji, T. (2018) Numerical Simulation of Liquid Crystal Flow Induced by Annihilation of a Pair of Defects. Open Journal of Fluid Dynamics, 8, 343-360. https://doi.org/10.4236/ojfd.2018.84022

References

Appendix 1

The time evolution equation for coefficients C_lm is expressed as

$\begin{array}{c}\frac{\partial C_{lm}}{\partial t}=-u\frac{\partial C_{lm}}{\partial x}-w\frac{\partial C_{lm}}{\partial z}-l\left(l+1\right){\left(1-S^2\right)}^{-2}{C}_{lm}-\frac{3}{4}U{\left(1-S^2\right)}^{-2}\\ \times [\left\{\sqrt{\frac{\text{2π}}{15}}\left(D_{22}+D_{2-2}\right)-\sqrt{\frac{\text{4π}}{45}}{D}_{20}\right\}{\left(-1\right)}^m\sqrt{\frac{\text{8π}}{15}}{\displaystyle \underset{l^{\prime }{m}^{\prime }}{\sum }{C}_{l^{\prime }{m}^{\prime }}\{A\left(\Vert l,-m+1;2,1;l^{\prime },m^{\prime }\Vert -\Vert l,-m+1;2,-1;l^{\prime },m^{\prime }\Vert \right)}\\ -B\left(\Vert l,-m-1;2,1;l^{\prime },m^{\prime }\Vert -\Vert l,-m-1;2,-1;l^{\prime },m^{\prime }\Vert \right)-2m\left(\Vert l,-m;2,2;l^{\prime },m^{\prime }\Vert -\Vert l,-m;2,-2;l^{\prime },m^{\prime }\Vert \right)\}\\ -\left\{\sqrt{\frac{\text{2π}}{15}}\left(D_{22}+D_{2-2}\right)+\sqrt{\frac{\text{4π}}{45}}{D}_{20}\right\}{\left(-1\right)}^m\sqrt{\frac{\text{8π}}{15}}{\displaystyle \underset{l^{\prime }{m}^{\prime }}{\sum }{C}_{l^{\prime }{m}^{\prime }}\{-A\left(\Vert l,-m+1;2,1;l^{\prime },m^{\prime }\Vert +\Vert l,-m+1;2,-1;l^{\prime },m^{\prime }\Vert \right)}\\ -B\left(\Vert l,-m-1;2,1;l^{\prime },m^{\prime }\Vert +\Vert l,-m-1;2,-1;l^{\prime },m^{\prime }\Vert \right)+2m\left(\Vert l,-m;2,2;l^{\prime },m^{\prime }\Vert -\Vert l,-m;2,-2;l^{\prime },m^{\prime }\Vert \right)\}\end{array}$

$\begin{array}{l}+\sqrt{\frac{\text{16π}}{45}}{D}_{20}{\left(-1\right)}^m\sqrt{\frac{\text{32π}}{15}}{\displaystyle \underset{l^{\prime }{m}^{\prime }}{\sum }{C}_{l^{\prime }{m}^{\prime }}\left(A\Vert l,-m+1;2,-1;l^{\prime },m^{\prime }\Vert +B\Vert l,-m-1;2,1;l^{\prime },m^{\prime }\Vert \right)}\\ -2\sqrt{\frac{\text{2π}}{15}}\left(D_{21}-D_{2-1}\right){\left(-1\right)}^m\sqrt{\frac{\text{8π}}{15}}{\displaystyle \underset{l^{\prime }{m}^{\prime }}{\sum }{C}_{l^{\prime }{m}^{\prime }}\{A\Vert l,-m+1;2,-2;l^{\prime },m^{\prime }\Vert -B\Vert l,-m-1;2,2;l^{\prime },m^{\prime }\Vert }\\ +m\left(\Vert l,-m;2,1;l^{\prime },m^{\prime }\Vert +\Vert l,-m;2,-1;l^{\prime },m^{\prime }\Vert \right)-\sqrt{\frac{3}{2}}\left(A\Vert l,-m+1;2,0;l^{\prime },m^{\prime }\Vert -B\Vert l,-m-1;2,0;l^{\prime },m^{\prime }\Vert \right)\}]\\ -\frac{\partial u}{\partial x}{\displaystyle \underset{l^{\prime }{m}^{\prime }}{\sum }{C}_{l^{\prime }{m}^{\prime }}(\sqrt{\frac{\text{π}}{30}}{A}^{\prime }{\Vert l,m;2,1;l^{\prime },m^{\prime }-1\Vert }^{*}-\sqrt{\frac{\text{π}}{30}}{B}^{\prime }{\Vert l,m;2,1;l^{\prime },m^{\prime }+1\Vert }^{*}-\sqrt{\frac{\text{π}}{30}}{A}^{\prime }{\Vert l,m;2,-1;l^{\prime },m^{\prime }-1\Vert }^{*}}\\ +\sqrt{\frac{\text{π}}{30}}{B}^{\prime }{\Vert l,m;2,-1;l^{\prime },m^{\prime }+1\Vert }^{*}-\sqrt{\frac{\text{2π}}{15}}{m}^{\prime }{\Vert l,m;2,2;l^{\prime },m^{\prime }\Vert }^{*}+\sqrt{\frac{\text{2π}}{15}}{m}^{\prime }{\Vert l,m;2,-2;l^{\prime },m^{\prime }\Vert }^{*}\end{array}$

$\begin{array}{l}+\sqrt{\frac{\text{4π}}{5}}{\Vert l,m;2,0;l^{\prime },m^{\prime }\Vert }^{*}-\sqrt{\frac{\text{6π}}{5}}{\Vert l,m;2,2;l^{\prime },m^{\prime }\Vert }^{*}-\sqrt{\frac{\text{6π}}{5}}{\Vert l,m;2,-2;l^{\prime },m^{\prime }\Vert }^{*})\\ -\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right){\displaystyle \underset{l^{\prime }{m}^{\prime }}{\sum }{C}_{l^{\prime }{m}^{\prime }}\{\sqrt{\frac{\text{4π}}{45}}{B}^{\prime }{\Vert l,m;2,0;l^{\prime },m^{\prime }+1\Vert }^{*}-\sqrt{\frac{\text{4π}}{45}}{A}^{\prime }{\Vert l,m;2,0;l^{\prime },m^{\prime }-1\Vert }^{*}}\\ +\sqrt{\frac{\text{2π}}{15}}{m}^{\prime }\left({\Vert l,m;2,1;l^{\prime },m^{\prime }\Vert }^{*}+{\Vert l,m;2,-1;l^{\prime },m^{\prime }\Vert }^{*}\right)+\sqrt{\frac{\text{6π}}{5}}\left({\Vert l,m;2,1;l^{\prime },m^{\prime }\Vert }^{*}-{\Vert l,m;2,-1;l^{\prime },m^{\prime }\Vert }^{*}\right)\}\\ -\left(\frac{1}{6}\frac{\partial u}{\partial z}-\frac{1}{3}\frac{\partial w}{\partial x}\right)\left(A{C}_{lm-1}-B{C}_{lm+1}\right)+\frac{\partial w}{\partial z}{\displaystyle \underset{l^{\prime }{m}^{\prime }}{\sum }{C}_{l^{\prime }{m}^{\prime }}(\sqrt{\frac{\text{2π}}{15}}{A}^{\prime }{\Vert l,m;2,1;l^{\prime },m^{\prime }-1\Vert }^{*}}\\ +\sqrt{\frac{\text{2π}}{15}}{B}^{\prime }{\Vert l,m;2,-1;l^{\prime },m^{\prime }+1\Vert }^{*}+\sqrt{\frac{\text{16π}}{5}}{\Vert l,m;2,0;l^{\prime },m^{\prime }\Vert }^{*})\end{array}$

where

$D_{lm}=C_{lm}+\frac{l_i^2}{24}\left(\frac{{\partial }^2{C}_{lm}}{\partial x^2}+\frac{{\partial }^2{C}_{lm}}{\partial z^2}\right)=C_{lm}+\frac{l_i^2}{24}{\nabla }^2{C}_{lm}$

$\Vert l_1,m_1;l_2,m_2;l_3,m_3\Vert ={\displaystyle {\int }_{\left|u\right|=1}{Y}_{l_1{m}_1}{Y}_{l_2{m}_2}{Y}_{l_3{m}_3}\text{d}\Omega }$

${\Vert l_1,m_1;l_2,m_2;l_3,m_3\Vert }^{*}={\displaystyle {\int }_{\left|u\right|=1}{Y}_{l_1{m}_1}^{*}{Y}_{l_2{m}_2}{Y}_{l_3{m}_3}\text{d}\Omega }$

$S^2=\frac{3}{2}S:S$

$A=\sqrt{\left(l+m\right)\left(l-m+1\right)}$

$B=\sqrt{\left(l-m\right)\left(l+m+1\right)}$

$A^{\prime }=\sqrt{\left(l^{\prime }+m^{\prime }\right)\left(l^{\prime }-m^{\prime }+1\right)}$

and

$B^{\prime }=\sqrt{\left(l^{\prime }-m^{\prime }\right)\left(l^{\prime }+m^{\prime }+1\right)}$

The non-zero components of the order parameter tensor S and fourth-order tensor Q are expressed in terms of C_lm as

$S_{xx}=\sqrt{\frac{\text{2π}}{15}}\left(C_{22}+C_{2-2}\right)-\sqrt{\frac{\text{4π}}{45}}{C}_{20}$

$S_{yy}=-\sqrt{\frac{\text{2π}}{15}}\left(C_{22}+C_{2-2}\right)-\sqrt{\frac{\text{4π}}{45}}{C}_{20}$

$S_{zz}=\sqrt{\frac{\text{16π}}{45}}{C}_{20}$

$S_{xz}=-\sqrt{\frac{\text{2π}}{15}}\left(C_{21}-C_{2-1}\right)$

$\begin{array}{c}{Q}_{xxxx}=\sqrt{\frac{\text{2π}}{315}}\left(C_{44}+C_{4-4}\right)-\sqrt{\frac{\text{8π}}{2205}}\left(C_{42}+C_{4-2}\right)+\sqrt{\frac{\text{4π}}{1225}}{C}_{40}\\ +\sqrt{\frac{\text{72π}}{735}}\left(C_{22}+C_{2-2}\right)-\sqrt{\frac{\text{16π}}{245}}{C}_{20}+\sqrt{\frac{\text{4π}}{25}}{C}_{00}\end{array}$

$Q_{xxxz}=-\sqrt{\frac{\text{π}}{315}}\left(C_{43}-C_{4-3}\right)+\sqrt{\frac{\text{π}}{245}}\left(C_{41}-C_{4-1}\right)-\sqrt{\frac{\text{6π}}{245}}\left(C_{21}-C_{2-1}\right)$

$Q_{xxyy}=-\sqrt{\frac{\text{2π}}{315}}\left(C_{44}+C_{4-4}\right)+\sqrt{\frac{\text{4π}}{11025}}{C}_{40}-\sqrt{\frac{\text{16π}}{2205}}{C}_{20}+\sqrt{\frac{\text{4π}}{225}}{C}_{00}$

$\begin{array}{c}{Q}_{xxzz}=\sqrt{\frac{\text{8π}}{2205}}\left(C_{42}+C_{4-2}\right)-\sqrt{\frac{\text{64π}}{11025}}{C}_{40}+\sqrt{\frac{\text{2π}}{735}}\left(C_{22}+C_{2-2}\right)\\ +\sqrt{\frac{\text{4π}}{2205}}{C}_{20}+\sqrt{\frac{\text{4π}}{225}}{C}_{00}\end{array}$

$Q_{xyyz}=\sqrt{\frac{\text{π}}{315}}\left(C_{43}-C_{4-3}\right)+\sqrt{\frac{\text{π}}{2205}}\left(C_{41}-C_{4-1}\right)-\sqrt{\frac{\text{2π}}{735}}\left(C_{21}-C_{2-1}\right)$

$Q_{xzzz}=-\sqrt{\frac{\text{16π}}{2205}}\left(C_{41}-C_{4-1}\right)-\sqrt{\frac{\text{6π}}{245}}\left(C_{21}-C_{2-1}\right)$

$\begin{array}{c}{Q}_{yyzz}=-\sqrt{\frac{\text{8π}}{2205}}\left(C_{42}+C_{4-2}\right)-\sqrt{\frac{\text{64π}}{11025}}{C}_{40}-\sqrt{\frac{\text{2π}}{735}}\left(C_{22}+C_{2-2}\right)\\ +\sqrt{\frac{\text{4π}}{2205}}{C}_{20}+\sqrt{\frac{\text{4π}}{225}}{C}_{00}\end{array}$

and

$Q_{zzzz}=\sqrt{\frac{\text{256π}}{11025}}{C}_{40}+\sqrt{\frac{\text{64π}}{245}}{C}_{20}+\sqrt{\frac{\text{4π}}{25}}{C}_{00}$

Appendix 2

The first normal stress difference N₁ and the shear stresses τ_xz and τ_zx are expressed in terms of S and Q as

$\begin{array}{c}\frac{N_1}{3ckT}=S_{xx}-S_{zz}-U(S_{xx}^2-S_{zz}^2+\frac{1}{3}{S}_{xx}-\frac{1}{3}{S}_{zz}-S_{xx}{Q}_{xxxx}+S_{xx}{Q}_{xxzz}-S_{yy}{Q}_{xxyy}\\ \underset{}{\overset{}{}}+S_{yy}{Q}_{yyzz}-S_{zz}{Q}_{xxzz}+S_{zz}{Q}_{zzzz}-2S_{xz}{Q}_{xxxz}+2S_{xz}{Q}_{xzzz})\\ -\frac{U{l}_i}{24}[S_{xx}{\nabla }^2{S}_{xx}-S_{zz}{\nabla }^2{S}_{zz}+\frac{1}{3}{\nabla }^2\left(S_{xx}-S_{zz}\right)-\left(Q_{xxxx}-Q_{xxzz}\right){\nabla }^2{S}_{xx}\\ -\left(Q_{xxyy}-Q_{yyzz}\right){\nabla }^2{S}_{yy}-\left(Q_{xxzz}-Q_{zzzz}\right){\nabla }^2{S}_{zz}-2\left(Q_{xxxz}-Q_{xzzz}\right){\nabla }^2{S}_{xz}\end{array}$

$\begin{array}{l}+\frac{1}{4}\{{\left(\frac{\partial S_{xx}}{\partial x}\right)}^2-{\left(\frac{\partial S_{xx}}{\partial z}\right)}^2+{\left(\frac{\partial S_{yy}}{\partial x}\right)}^2-{\left(\frac{\partial S_{yy}}{\partial z}\right)}^2+{\left(\frac{\partial S_{zz}}{\partial x}\right)}^2-{\left(\frac{\partial S_{zz}}{\partial z}\right)}^2\\ +2{\left(\frac{\partial S_{xz}}{\partial x}\right)}^2-2{\left(\frac{\partial S_{xz}}{\partial z}\right)}^2-\left(\frac{{\partial }^2{S}_{xx}}{\partial x^2}-\frac{{\partial }^2{S}_{xx}}{\partial z^2}\right)S_{xx}-\left(\frac{{\partial }^2{S}_{yy}}{\partial x^2}-\frac{{\partial }^2{S}_{yy}}{\partial z^2}\right)S_{yy}\\ -\left(\frac{{\partial }^2{S}_{zz}}{\partial x^2}-\frac{{\partial }^2{S}_{zz}}{\partial z^2}\right)S_{zz}-2\left(\frac{{\partial }^2{S}_{xz}}{\partial x^2}-\frac{{\partial }^2{S}_{xz}}{\partial z^2}\right)S_{xz}\}]\\ +\frac{\zeta }{6kT}\left\{\frac{\partial u}{\partial x}\left(Q_{xxxx}-Q_{xxzz}\right)+\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\left(Q_{xxxz}-Q_{xzzz}\right)+\frac{\partial w}{\partial z}\left(Q_{xxzz}-Q_{zzzz}\right)\right\}\end{array}$

$\begin{array}{c}\frac{{\tau }_{xz}}{3ckT}=S_{xz}-U\left\{\left(S_{xx}+S_{zz}\right)S_{xz}+\frac{1}{3}{S}_{xz}-S_{xx}{Q}_{xxxz}-S_{yy}{Q}_{xyyz}-S_{zz}{Q}_{xzzz}-2S_{xz}{Q}_{xxzz}\right\}\\ -\frac{U{l}_i^2}{24}\{\underset{\_}{\underset{\_}{S_{xx}{\nabla }^2{S}_{xz}+S_{xz}{\nabla }^2{S}_{zz}}}+\frac{1}{3}{\nabla }^2{S}_{xz}-Q_{xxxz}{\nabla }^2{S}_{xx}-Q_{xyyz}{\nabla }^2{S}_{yy}\\ -Q_{xzzz}{\nabla }^2{S}_{zz}-2Q_{xxzz}{\nabla }^2{S}_{xz}+\frac{1}{4}(\frac{\partial S_{xx}}{\partial x}\frac{\partial S_{xx}}{\partial z}+\frac{\partial S_{yy}}{\partial x}\frac{\partial S_{yy}}{\partial z}+\frac{\partial S_{zz}}{\partial x}\frac{\partial S_{zz}}{\partial z}\\ +2\frac{\partial S_{xz}}{\partial x}\frac{\partial S_{xz}}{\partial z}-\frac{{\partial }^2{S}_{xx}}{\partial x\partial z}{S}_{xx}-\frac{{\partial }^2{S}_{yy}}{\partial x\partial z}{S}_{yy}-\frac{{\partial }^2{S}_{zz}}{\partial x\partial z}{S}_{zz}-2\frac{{\partial }^2{S}_{xz}}{\partial x\partial z}{S}_{xz})\}\\ +\frac{\zeta }{6kT}\left\{\frac{\partial u}{\partial x}{Q}_{xxxz}+\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)Q_{xxzz}+\frac{\partial w}{\partial z}{Q}_{xzzz}\right\}\end{array}$

and

$\frac{{\tau }_{zx}}{3ckT}=\cdots -\frac{U{l}_i^2}{24}\left\{\underset{\_}{\underset{\_}{S_{xz}{\nabla }^2{S}_{xx}+S_{zz}{\nabla }^2{S}_{xz}}}+\cdots \right\}+\cdots$

The difference between τ_xz and τ_zx is only the underlined term.

Appendix 3

With respect to the director n, let us define the angle between n and the x axis as θ_c, so that the angular momentum equation of the Leslie-Ericksen theory with the one-constant approximation of the elastic constants in the molecular field reduces simply to

${\nabla }^2{\theta }_c=0$ (A1)

in an equilibrium state [2] . When defects are included, a fundamental solution of Equation (A1) is

${\theta }_c=s{\tan }^{-1}\left(\frac{z}{x}\right)$ , (A2)

where s is the defect strength. Since Equation (A1) is linear, the superposition of solutions is effective. Thus, when a defect with s = +1/2 exists at P (x₁, z₁) and a defect with s = −1/2 exists at Q (x₂, z₂) in Figure 2, the director distribution around the defect cores is

${\theta }_c=\frac{1}{2}{\tan }^{-1}\frac{z-z_1}{x-x_1}-\frac{1}{2}{\tan }^{-1}\frac{z-z_2}{x-x_2}$ . (A3)

Finally, we have modified the values of Equation (A3) so that they fit the periodic boundary condition. We denote the initial distribution of the director by line segments in Figure 2.

For the scalar order parameter, we set S = 0 at defect cores and S = S^eq at the other region. S^eq depends only on the nematic potential intensity U, and is obtained from Equation (1) without flow terms. For example, S^eq = 0.615 at U = 5.