Phase Transitions in Nematics: Textures with Tactoids and Disclinations

TL;DR

This paper models the isotropic-to-nematic phase transition in liquid crystals using a generalized Landau-de Gennes framework, successfully reproducing textures with tactoids and disclinations observed experimentally.

Contribution

It introduces a specific elastic constant combination in the Landau-de Gennes theory to accurately simulate phase transition textures.

Findings

Qualitative reproduction of experimentally observed nematic textures.

Successful modeling of interfaces and topological defects.

Validation of the generalized Landau-de Gennes approach.

Abstract

We demonstrate that a first order isotropic-to-nematic phase transition in liquid crystals can be succesfully modeled within the generalized Landau-de Gennes theory by selecting an appropriate combination of elastic constants. The numerical simulations of the model established in this paper qualitatively reproduce the experimentally observed configurations that include interfaces and topological defects in the nematic phase.

Figures

Equations

$${\lambda }_1+{\lambda }_2+{\lambda }_3=0$$

$$Q=-\frac{{\lambda }_3}{2}\mathbf{l}\otimes \mathbf{l}-\frac{{\lambda }_3}{2}\mathbf{m}\otimes \mathbf{m}+{\lambda }_3\mathbf{n}\otimes \mathbf{n}=S\left(\mathbf{n}\otimes \mathbf{n}-\frac{1}{3}\mathbf{I}\right),$$

$${\lambda }_i\in [-1/3,2/3],\mathrm{for}i=1,2,3.$$

$$\begin{gathered} f_e(\nabla Q):=\frac{L_1}{2}{|\nabla Q|}^2+\frac{L_2}{2}{Q}_{ij,j}{Q}_{ik,k}+\frac{L_3}{2}{Q}_{ik,j}{Q}_{ij,k} \\ =\sum _{j=1}^3\left\{\frac{L_1}{2}{|\nabla Q_j|}^2+\frac{L_2}{2}{\left(\mathrm{div}{Q}_j\right)}^2+\frac{L_3}{2}\nabla Q_j\cdot \nabla Q_j^T\right\}. \end{gathered}$$

$$W(Q):=3a\mathrm{tr}\left(Q^2\right)-2b\mathrm{tr}\left(Q^3\right)+\frac{1}{4}{\left(\mathrm{tr}\left(Q^2\right)\right)}^2,$$

$$E[Q]:={\int }_{\Omega }\left\{f_e(Q,\nabla Q)+W(Q)\right\}dV.$$

$$A:=\left\{Q\in M^{3\times 3}:Q^T=Q,\mathrm{tr}Q=0\right\}$$

$$W({\lambda }_1,{\lambda }_2)=6a\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_1{\lambda }_2\right)+6b{\lambda }_1{\lambda }_2\left({\lambda }_1+{\lambda }_2\right)+{\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_1{\lambda }_2\right)}^2.$$

$$\frac{\partial W}{\partial {\lambda }_1}$$

$$\frac{\partial W}{\partial {\lambda }_2}$$

$$\frac{{\partial }^{2}W}{\partial {\lambda }_1^2}$$

$$\frac{{\partial }^{2}W}{\partial {\lambda }_2^2}$$

$$\frac{{\partial }^{2}W}{\partial {\lambda }_1\partial {\lambda }_2}$$

$$a+b\lambda +{\lambda }^2=0.$$

$$\lambda =\frac{-b\pm \sqrt{b^2-4a}}{2},$$

$${\lambda }_m=\frac{-b-\sqrt{b^2-4a}}{2}$$

$${\lambda }_s=\frac{-b+\sqrt{b^2-4a}}{2}$$

$$a=a_0:=2b^2/9.$$

$$W_m:=W({\lambda }_m,{\lambda }_m)=9a{b}^2-9a^2-\frac{3}{2}{b}^4-\frac{3}{2}b{\left(b^2-4a\right)}^{\frac{3}{2}}$$

$$W(Q)\to W(Q)-W_m.$$

$$E_{OF}[n]=\frac{1}{2}{\int }_{\Omega }\left\{K_1{\left(\mathrm{div}n\right)}^2+K_2{\left(n\cdot \mathrm{curl}n\right)}^2+K_3{\left|n\times \mathrm{curl}n\right|}^2\right\}dx$$

$$E_{OF}[n]=\frac{1}{2}{\int }_{\Omega }\left\{{\tilde{K}}_1{\left(\mathrm{div}n\right)}^2+{\tilde{K}}_2{\left|\nabla n\right|}^2\right\}dx,$$

$$Q={\lambda }_m\left(\mathrm{I}-3n\otimes n\right),$$

$$n\otimes n=\frac{1}{3}\left(\mathrm{I}-Q/{\lambda }_m\right)$$

$${|\nabla n|}^2=\frac{1}{18}{|\nabla (Q/{\lambda }_m)|}^2$$

$$-\frac{1}{3}\mathrm{div}(Q/{\lambda }_m)=\mathrm{div}\left(n\otimes n\right)=(\mathrm{div}n)n+\nabla nn.$$

$$\left(n\otimes n\right)\nabla nn=\left(n\cdot \nabla nn\right)n=\left(\nabla n^{T}n\cdot n\right)n=\frac{1}{2}\left(\nabla \left(|n{|}^2\right)\cdot n\right)n=0,$$

$$(\mathrm{div}n)n=\left(n\otimes n\right)\left((\mathrm{div}n)n+\nabla nn\right)=-\frac{1}{9}\left(\mathrm{I}-Q/{\lambda }_m\right)\mathrm{div}(Q/{\lambda }_m)$$

$${\tilde{K}}_1{\left(\mathrm{div}n\right)}^2+{\tilde{K}}_2{\left|\nabla n\right|}^2\sim \frac{L_1}{2}{|\nabla (Q/{\lambda }_m)|}^2+\frac{L_2}{2}{\left|\left(\mathrm{I}-Q/{\lambda }_m\right)\mathrm{div}(Q/{\lambda }_m)\right|}^2$$

$$F[Q]:={\int }_{\Omega }\left(\frac{L_1}{2}{|\nabla (Q/{\lambda }_m)|}^2+\frac{L_2}{2}{\left|\left(\mathrm{I}-Q/{\lambda }_m\right)\mathrm{div}(Q/{\lambda }_m)\right|}^2+w_{0}W(Q)\right)dx,$$

$${\left|\left(\mathrm{I}-Q/{\lambda }_m\right)\mathrm{div}(Q/{\lambda }_m)\right|}^2=\frac{1}{{\lambda }_m^2}{\left|\mathrm{div}Q\right|}^2-\frac{2}{{\lambda }_m^3}Q\mathrm{div}Q\cdot \mathrm{div}Q+\frac{1}{{\lambda }_m^4}{\left|Q\mathrm{div}Q\right|}^2.$$

$$\epsilon :=\frac{L_1}{L_2},\gamma :=\frac{4L_2}{w_0{l}^2{\lambda }_m^4}.$$

Full text

1 Introduction

Phase transitions between ordered and isotropic states in nematic liquid crystals are of interest both for technological reasons as well as because nematics offer one of the simplest experimental systems where interfaces and topological defects may coexist. In particular, a uniaxial nematic liquid crystal typically consists of asymmetrically-shaped molecules that under appropriate conditions tend to align roughly in the same direction. The preferred molecular orientation in a given region occupied by a uniaxial nematic is often described in terms of the unit vector field $n$, called the director. A more precise description takes into account the non-polar character of ordering, e.g., because the probabilities of finding the head or the tail of a molecule in a given direction are equal. This description replaces the director $n$ with a projection matrix $n\otimes n$ which can be identified with an element of the projective space and is invariant with respect to reversal of orientation $n\to-n$.

Two different types of topological defects can be present in a nematic in three dimensions: [math]-dimensional point defects also known as vortices or nematic hedgehogs and $1$-dimensional disclinations. Point defects may also be present on the surface of a nematic; such defects are known as boojums. There is a large body of literature devoted to the study of nematic defects and we refer the reader to [1] for a comprehensive exposition on the topic. Our principal interest in this work is to examine the interaction between the nematic defects and isotropic-to-nematic interfaces.

The interfaces form in the process of the first order isotropic-to-nematic phase transition that can be induced either by lowering temperature (thermotropic nematics) or increasing the concentration of asymmetric molecules in a solvent (lyotropic nematics). Typically, nuclei of the nematic phase form within the isotropic phase upon lowering the temperature. The nematic nuclei are separated from the isotropic regions by phase boundaries or interfaces. The interfaces subsequently propagate converting the isotropic phase to the nematic phase in the process. In addition to motion of interfaces, the resulting dynamics of the system involves formation, annihilation and propagation of various types of defects.

In a recent work [2], the authors examine the interplay between the interfaces and defects that are present during phase transitions in lyotropic chromonic liquid crystals (LCLC). The principal feature of the isotropic-to-nematic phase transition in LCLC is that the interface provides an “easy direction” for nematic anchoring on the interface which in this case forces the director field to be tangent to this phase boundary. When combined with the anchoring (boundary) conditions on the walls of the container, the topology of the nematic configuration leads to formation of structural defects.

Our goal in the present paper is to demonstrate that the zoo of singularities observed in [2] can be described within the framework of the Landau-de Gennes model for $Q$-tensors related to the projection matrix descriptor of the nematic phase alluded to above. Critical to our modeling will be an assumption of large disparity between the values of the elastic constants appearing in the energy. In Section 2 we briefly describe the experimental observations that expand on some of the results presented in [2] by incorporating scenarios of the phase transition in which the isotropic phase regions, often called “negative tactoids” [3], [2], or simply “tactoids”, shrink into a uniform nematic state or a state with a topological defect, depending on the winding number of the tactoid. Note here that positive tactoids, i.e., regions of a nematic surrounded by an isotropic melt, have been previously theoretically treated in [4]-[9].

In Section 3 we review the basics of the Landau-de Gennes theory and then develop our model in Section 4. In Section 5 we describe our numerical results and compare them with experimental observations.

2 Experimental Results

The sample configuration consisted of the ITO glass that was spin-coated with a polyimide layer, SE7511, in which the directors of disodium chromoglycate (DSCG) are aligned parallel to substrates. Subsequently, two substrates were assembled into a cell with the thickness of 2 $\mu$m. $16$ wt% of DSCG solution was injected into the assembled cell.

The cell was cooled at the rate of $-0.4^{\circ}C/min$. As the temperature decreased from the isotropic phase, nematic domains appeared, grew, and coalesced. When many large nematic domains coalesced, they occasionally trapped isotropic islands, or tactoids, around which the director may have either zero or nonzero winding number. The snapshots of configurations with islands having different winding number of the director on their boundary are shown in Figs. 1-2. These are Polscope textures with color representing optical retardance, bars giving the orientation of the director, and circles indicating the point where orientation was measured.

3 Landau-de Gennes Model

3.1 The $Q$-tensor

Given a point $x\in\mathbb{R}^{3}$, the second moment of the orientational distribution of the rod-like nematic liquid crystal molecules near $x$ can be described by a $2$-tensor $Q(x)$ that takes the form of a $3\times 3$ symmetric, traceless matrix. By virtue of being symmetric and traceless, $Q$ has three real eigenvalues $\lambda_{1},\ \lambda_{2},$ and $\lambda_{3}$ that satisfy

[TABLE]

and a mutually orthonormal eigenframe $\left\{\mathbf{l},\mathbf{m},\mathbf{n}\right\}$. While a more detailed overview of the theory can be found in [10], below we outline the key elements that will be needed in the developments below.

Suppose that $\lambda_{1}=\lambda_{2}=-\lambda_{3}/2.$ Then the liquid crystal is in a uniaxial nematic state and

[TABLE]

where $S:=\frac{3\lambda_{3}}{2}$ is the uniaxial nematic order parameter and $\mathbf{n}\in\mathbb{S}^{2}$ is the nematic director. If there are no repeated eigenvalues, the liquid crystal is said to be in a biaxial nematic state, while it is in the isotropic state when all three eigenvalues and hence $Q$ itself vanish. The isotropic state of a liquid crystal is associated, for instance, with a high temperature regime.

From the modeling perspective it turns out that the eigenvalues of $Q$ must satisfy the constraints [11, 12]:

[TABLE]

3.2 Landau-de Gennes energy

Suppose that the bulk elastic energy density of a nematic liquid crystal is given by a frame-indifferent expression $f_{e}(Q,\nabla Q)$. A common expression that includes terms up to the second order in $Q$ reads

[TABLE]

The bulk Landau-de Gennes energy density is

[TABLE]

cf. [10]. Here $Q_{j},\,j=1,2,3$ is the $j$-th column of the matrix $Q$ and $A\cdot B=\mathrm{tr}\,{\left(B^{T}A\right)}$ is the dot product of two matrices $A,B\in M^{3\times 3}.$ Further, the coefficient $a$ is temperature-dependent and in particular is negative for sufficiently low temperatures. One readily checks that the form (4) of this potential implies that in fact $W$ depends only on the eigenvalues of $Q$, and due to the trace-free condition, therefore depends only on two eigenvalues. The form of $W$ guarantees that the isotropic state $Q\equiv 0$ yields a global minimum at high temperatures while a uniaxial state of the form (1) gives the minimum when temperature (i.e. the parameter $a$) is reduced below a certain critical value, cf. [13, 10]. We remark for future use that $W$ is bounded from below and can be made non-negative by adding an appropriate constant.

We now turn to the behavior of the nematic on the boundary of the sample. Here two alternatives are possible. First, the Dirichlet boundary conditions on $Q$ are referred to as strong anchoring conditions in the physics literature: they impose specific preferred orientations on nematic molecules on surfaces bounding the liquid crystal. An alternative is to specify the anchoring energy on the boundary of the sample; then orientations of the molecules on the boundary are determined as a part of the minimization procedure. This approach is known as weak anchoring.

Putting the energy densities together, cf. (3), (4), we arrive at a Landau-de Gennes type model to be analyzed in this study, given by

[TABLE]

4 Model Development

In this section we derive a version of the Landau-de Gennes model that is appropriate for the modeling of nematic systems with disparate elastic constants. In particular, we are interested in the case when the elastic constant corresponding to splay deformations is larger than those for bend and twist so that the splay of the director is relatively expensive. Note that this situation can be found in experimental systems, such as chromonic lyotropic liquid crystals [14] shown in Figures 1-2 and in thermotropic nematics of certain molecular shape, such as dimers [15]. To this end, let

[TABLE]

and consider the Landau-de Gennes potential $W(Q)$ defined in Eq. 4. As long as $Q\in A$, one finds that this potential depends only on the eigenvalues of $Q$, say $\lambda_{1}$ and $\lambda_{2}$ with $\lambda_{3}=-(\lambda_{1}+\lambda_{2})$, and with a slight abuse of notation we arrive at

[TABLE]

We are interested in describing the extrema of $W(\lambda_{1},\lambda_{2})$. Taking the derivatives of $W$ gives

[TABLE]

It is easy to check that the critical points of $W$ are $(\lambda_{1},\lambda_{2})=(0,0)$ as well as $(\lambda_{1},\lambda_{2})=(\lambda,-2\lambda)$, $(\lambda_{1},\lambda_{2})=(-2\lambda,\lambda)$, and $(\lambda_{1},\lambda_{2})=(\lambda,\lambda)$, where $\lambda$ solves

[TABLE]

The point $(\lambda_{1},\lambda_{2})=(0,0)$ corresponds to a local minimum of $W$ if $a>0$ and to a local maximum of $W$ if $a<0$. The solutions of Eq. 6 are given by

[TABLE]

hence there are no nematic critical points whenever $a>b^{2}/4$ and $W$ has a single extremum corresponding to the global minimum at $(\lambda_{1},\lambda_{2})=(0,0)$. By considering second derivatives of $W$, we also find that

[TABLE]

is a point of a local minimum and

[TABLE]

is a saddle point of $W$ if $a<b^{2}/4$. We conclude that the nematic minima at $(\lambda_{m},-2\lambda_{m})$, $(-2\lambda_{m},\lambda_{m})$, and $(\lambda_{m},\lambda_{m})$ coexist with the isotropic minimum at $(0,0)$ as long as $0<a<b^{2}/4$. An easy computation shows that all of these minima have the same depth of [math] when

[TABLE]

Because here we will be interested in the regime $0\leq a\leq a_{0}$ when the energy value in the isotropic state is greater than or equal to the minimum energy in a nematic state, we subtract

[TABLE]

from $W$ to ensure that the global minimum value of the Landau-de Gennes energy is [math]. Hence, from now on

[TABLE]

Note that decreasing $a$ from $a_{0}$ to [math] corresponds to increased undercooling and a larger thermodynamic force driving the isotropic-to-nematic phase transition.

Let $\Omega\subset\mathbb{R}^{3}$ denote the region occupied by the liquid crystal sample. To derive an expression for the elastic energy, recall first that the Oseen-Frank energy defined over director fields $n\in H^{1}\left(\Omega,\mathbb{S}^{2}\right)$ is given by

[TABLE]

where for simplicity we assume that the admissible competitors are subject to the same Dirichlet boundary data on $\partial\Omega$. Assuming in this paper that $K_{2}=K_{3}<K_{1}$, the energy can be written (up to null Lagrangian terms) as

[TABLE]

where $\tilde{K}_{1}=K_{1}-K_{2}$ and $\tilde{K}_{2}=K_{2}$ are nonnegative elastic constants. We want to derive a well-posed variational problem for $Q$-valued fields such that the elastic energy for $Q$ formally reduces to Eq. 10 whenever $Q$ is in a nematic state that minimizes the potential energy $W$. To this end, the minimum of $W$ is achieved whenever $Q$ has eigenvalues $\lambda_{m}$, $\lambda_{m}$, and $-2\lambda_{m}$ so that $Q$ can be written as

[TABLE]

where $n$ is a unit eigenvector corresponding to the eigenvalue $-2\lambda_{m}$. From this equation,

[TABLE]

and it can be easily checked that

[TABLE]

while

[TABLE]

Making use of the fact that

[TABLE]

it follows that

[TABLE]

and

[TABLE]

where $L_{1}:=\frac{\tilde{K}_{2}}{9}$ and $L_{2}:=\frac{2\tilde{K}_{1}}{81}$. Note that here we switched the indices of the elastic constants in order to conform with the standard notation (cf. Eq. 3).

The total energy of a nematic configuration as described within the Landau-de Gennes $Q$-tensor theory in this work will be given by

[TABLE]

where $w_{0}$ is a constant that has units of energy per unit volume.

Remark 4.1**.**

Expanding the second elastic term in Eq. 11 we obtain

[TABLE]

*Comparing the terms appearing in this expression with the elastic invariants of the generalized Landau-de Gennes theory in [16], we see that these terms correspond to $L_{2}^{(2)}$-, $L_{3}^{(3)}$-, and $L_{6}^{(4)}$-invariants. *

Now let $l>0$ denote a characteristic length of the problem and set

[TABLE]

Before proceeding further, we introduce the following rescalings

[TABLE]

and drop the bar notation to obtain

[TABLE]

where

[TABLE]

In this scaling, the potential $W$ given by (14) is now minimized by any symmetric traceless matrix with eigenvalues $-1,\ -1,\ 2$, and the global minimum value of $W$ is equal to zero.

In our simulations we consider a simplified form of a $Q$-tensor that can be obtained via a dimension reduction procedure for thin nematic films [17]. In the corresponding ansatz, one eigenvector of admissible $Q$-tensors must be perpendicular to the plane of the film. Then in the system of coordinates in which the normal to the film is parallel to the $z$-axis, the $Q$-tensor is independent of $z$ and can be written [18] as

[TABLE]

for a scalar-valued function $\beta$ and vector-valued function $u=(u_{1},u_{2})$. Let

[TABLE]

so that

[TABLE]

where $\mathrm{I}_{2}$ is the $2\times 2$ identity matrix. In terms of $\beta$ and $u$ the contributions to the energy density are as follows

[TABLE]

[TABLE]

and

[TABLE]

Suppose that the film occupies the region $\omega\in\mathbb{R}^{2}$. The system of Euler-Lagrange PDEs corresponding to the energy functional Eq. 13 then takes the form

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

are the Pauli matrices.

When the liquid crystal is in the energy minimizing nematic state while the director $n=(\cos{\theta},\sin{\theta},0)$ lies in the plane of the film, we have

[TABLE]

It follows from this computation that

[TABLE]

hence the vector $u$ winds twice as fast along a given curve in $\bar{\omega}$, compared to the director $n$ when $n$ is planar.

More generally, the ansatz Eq. 15 splits the nematic component of the minimal set of the bulk energy $W$ into two disconnected components

[TABLE]

and

[TABLE]

Here the first component corresponds to a constant nematic state with the director perpendicular to the surface of the film. The second component includes all configurations with the director lying in the plane of the film—these configurations can be nontrivial and, in particular, they include the director fields that carry a nonzero winding number.

Since we are interested in regimes when both nematic and isotropic phases coexist, we recall that whenever

[TABLE]

the set

[TABLE]

yields a local minimum of $W$ corresponding to the isotropic phase. The corresponding local minimum energy value is greater than or equal to the global minimum value of $W$. When $a=\frac{2b^{2}}{9}$, the set

[TABLE]

gives the third connected component of the minimal set of the bulk energy.

Given that $W$ has a multi-component minimal set, the energy Eq. 13 is of Allen-Cahn-type as long as the factor in front of $W$ is large. To this end, in what follows we assume that $\gamma=\varepsilon$ and $\varepsilon>0$ is small. We will consider gradient flow dynamics associated with this model, but before we proceed further it is worth comparing our present situation to the more familiar one of gradient flow dynamics for a multi-well Allen-Cahn type potential with diffusion given simply by the Laplacian. In this scaling, when the diffusion is given by the Laplacian, formal asymptotics suggest that for Allen-Cahn dynamics one should expect an interface propagating by curvature flow that separates the different phases defined by components of the zero set of the potential. In the scalar setting of Allen-Cahn where a double well potential vanishes at two points, there are by now numerous rigorous proofs of this fact based on maximum principles, barriers and/or comparison principles, see e.g. [19, 20, 21, 22] as well as the energetic argument in the radial setting in [23]. For the vector setting of time-dependent Allen-Cahn, formal asymptotics based on multiple time-scale expansions again suggest that mean curvature flow emerges as the governing equation for the interface, [24, 25].

What distinguishes the dynamics in the present study, however, is the anisotropy of the diffusive terms indicated by the presence of the divergence terms in (16)–(18). While we again anticipate that in the regime $\varepsilon\ll 1$ an interface separating the different phases of the $Q$ tensor will evolve by a law involving curvature, the process will be significantly affected by the interaction between the director associated with $Q$ and the normal to the interface.

Our present model is closely related to the investigation in [26] of the $\varepsilon\to 0$ asymptotics for a director-like model based on an $\mathbb{R}^{2}$-valued order parameter $u$. There the elastic energy is similarly anisotropic and is coupled to a Chern-Simons-Higgs-type potential $|u|^{2}(|u|^{2}-1)^{2}$. Thus, the structure of the energy functional in [26] is similar to the setup proposed here in that it involves a potential with minima at the isotropic and a nematic state as well as the elastic terms that are quadratic in the gradient of the order parameter field. The term penalizing splay deformations in [26] dominates other elastic terms so that the divergence of the director is very expensive. The principal difference between the model considered in [26], which from now on we will refer to as the CSH-director model, and the present work is that here we consider non-orientable tensor fields, while the admissible fields in [26] are orientable.

Our numerical results indicate that similarity between the models leads to similar properties of critical points. In particular, one of the principal observations in [26] is that the director field is parallel to the interface in its immediate vicinity when the parameter $\varepsilon$ is small. The same behavior is exhibited by the numerically obtained critical points of the Landau-de Gennes energy Eq. 11, as will be demonstrated in the next section.

We introduce dynamics into the problem by assuming that evolution of the isotropic-to-nematic transition proceeds via gradient flow

[TABLE]

where $\mu_{u}>0$ and $\mu_{\beta}>0$ are inverses of the scalar mobilities of $u$ and $\beta$, respectively. This gives the following systems of PDEs

[TABLE]

[TABLE]

[TABLE]

In order to simulate experimentally available configurations, we will set the initial configuration of $(u,\beta)$ to consist of an isotropic region $\omega_{I}$ embedded in a nematic phase corresponding to the component $C_{2}^{N}$ of the minimal set of $W$. That is, we will assume that the initial condition is an appropriate mollification of $(u,\beta)$ such that

[TABLE]

Here $\chi_{\omega\backslash\omega_{I}}$ is the characteristic function of the region $\omega\backslash\omega_{I}$ occupied by the nematic phase. We also impose boundary conditions on $\partial\omega$ with values in $C_{2}^{N}$; we will assume that the boundary data for $u$ may have a nonzero winding number.

In the next section we present the result of simulations for the system Eq. 20-Eq. 22 and compare the numerical outcomes to experimental observations. We remind the reader that our goal in the present paper is to demonstrate numerically that behavior of interfaces and singularities observed in experiments with chromonic liquid crystals can be qualitatively recovered within the framework of the Landau-de Gennes model. We leave for future work both formal and rigorous analysis of our results, as well as any quantitative considerations and associated modifications of the model.

5 Simulations vs Experiment

The simulations in this section are performed using COMSOL Multiphysics® [27] with the domain $\omega$ taken to be a disk of nondimensional radius $1/4$, centered at the origin. In what follows, we set dimensional $b=1$ and let $\gamma=\varepsilon=0.06$. The initial data in all simulations is specified by choosing $\omega_{I}=\left\{(x,y)\in\omega:r(x,y)\leq\frac{1}{6}\right\}$ in Eq. 23. Both here and below $\left(r(x,y),\theta(x,y)\right)$ are polar coordinates of the point $(x,y)\in\mathbb{R}^{2}$.

Recall that when $a$ and $b$ are related via Eq. 8 (either in dimensional or nondimensional setting), the energies of the nematic and isotropic states are the same. Then, if in dimensional variables we set

[TABLE]

the parameter $\alpha\in\left(0,2b^{2}/9\right)$ measures the degree of undercooling in the system. That is, when $\alpha$ increases, the bulk energy density of the isotropic state grows with respect to that of the nematic state and thermodynamic forces driving the isotropic-to-nematic transition become stronger. Setting $\bar{\alpha}=\alpha/\lambda^{2}$ and dropping the bar, the nondimensional $a$ and $b$ are given by

[TABLE]

First, we verify the conjecture that the behavior of interfaces in our Landau-de Gennes system of PDEs with highly anisotropic elastic constants resembles that of interfaces in the related CSH-director system of [26]. To this end, we assume that $a$ is given by Eq. 8 so that $\alpha=0$ and both the nematic and the isotropic states are the global minimum of the Landau-de Gennes potential $W$. We simulate the system Eq. 20-Eq. 22 subject to the boundary conditions

[TABLE]

Taking into account Eq. 19, the director field on the boundary can be chosen as

[TABLE]

where the topological degree of $n$ is equal to $-1$ on $\partial\omega$.

The results of the gradient flow simulation for the degree $-1$ boundary data are shown in Fig. 3. An initially circular interface is seen to evolve into the shape given in the figure by a thick red curve that is also a contour line of $\beta=1/4$.

The system in Fig. 3 has reached a steady state that is very similar to what is observed for the CSH-director model in [26]. Indeed, for small $\varepsilon$, the energy of the isotropic/nematic configuration in both cases should include a penalty for divergence of the director in the nematic phase, as well as the cost for the perimeter of the interface. The anisotropy of the elastic constants also forces the director to align with the interface. The resulting director configuration and the shape of the interface are coupled; in particular, unlike the standard curvature flow, shrinking the perimeter and, thus, size of the isotropic region leads to unbounded growth of the elastic energy. The competition between the perimeter and elastic energy contributions enforces the equilibrium between the two phases observed in Fig. 3.

One is then tempted to ask whether there is any difference between the Landau-de Gennes and the CSH-director model [26] in terms of how the director field extracted from solutions of these two models should behave as $\varepsilon\to 0$. We expect that the difference would manifest itself when the size of an isotropic region is small enough so that this region can be thought of as a core of a topological vortex.

Suppose, for example, that the interface is a single smooth closed curve. Due to the director tangency condition on the interface and because the director field is orientable, the director has the winding number $1$ around the interface. When the island shrinks to a small size and because the director is $\mathbb{S}^{1}$-valued in the nematic phase, topological constraints would keep the isotropic island from disappearing completely in order to prevent the nematic configuration from having an infinite energy. The resulting degree $1$ vortex is stable in the CSH-director setting and would persist for a finite time, perhaps until it annihilates with another vortex of the opposite sign or collides with the boundary of $\omega$.

Now recall that within the Landau-de Gennes theory, a nematic state is described by a $Q$-tensor that is a translation and dilation of a projection matrix $n\otimes n$. The field $n\otimes n$ is not orientable and can be associated with an element of a projective plane $\mathbb{RP}^{2}$. In other words, by working with $n\otimes n$ instead of $n$, we identify the opposite directions $-n$ and $n$ as being the same. As the result, the smallest “quantum” of the winding number for a nematic tensor is $1/2$. Hence the degree $1$ vortex in the director description contains two “quanta” of degree if we switch to the Landau-de Gennes framework. It is a well-established fact that a higher degree vortex is unstable with respect to splitting into several vortices of smaller degrees since the energetic cost of a degree $d$ vortex is proportional to $d^{2}$, [28]. We expect that the degree $1$ CSH vortex would split into two degree $1/2$ vortices when considered in the sense of Landau-de Gennes. Our experimental observations conform to the latter picture; hence the description that disposes with orientability is necessary when considering an isotropic-to-nematic phase transition problem.

We now conduct numerical experiments to test whether evolution of tactoids observed in physical experiments can be captured within the Landau-de Gennes model.

5.1 Degree $1$ tactoid

Here we suppose that the director is parallel to the boundary of the disk $\omega$, e.g.,

[TABLE]

so that the topological degree of the orientable director field $n$ is equal $1$ on $\partial\omega$. In the non-orientable, $Q$-tensor setting the corresponding condition can be expressed as

[TABLE]

If the simulations are conducted at zero undercooling, when $a$ is given by Eq. 8, we observed both here and for the Ginzburg-Landau-type model [26] that, similar to Fig. 3, the isotropic domain evolving via gradient flow shrinks down to a certain size and then stabilizes. In order to drive this size down, we enforce larger undercooling by choosing $\alpha=0.2.$

The evolution of the nematic/isotropic interface is shown in Fig. 4.

While the isotropic tactoid is not too small, it maintains its circular shape and the director in the nematic region remains parallel to the boundary. Indeed, in this configuration the divergence contribution vanishes, the director is parallel to the interface, and the perimeter of the isotropic tactoid is minimal given its shape.

When the tactoid eventually shrinks to a vortex size seen in the right inset in the second row in Fig. 4, it loses stability and splits into two degree $1/2$ vortices that drift away from each other with time. The overall behavior is qualitatively similar to that of the right tactoid in Fig. 2.

5.2 Degree $-1$ tactoid

We now return to the boundary conditions

[TABLE]

considered in the beginning of this section for zero undercooling. Here we assume instead that $\alpha=-0.1$ and look for the effects on evolution of the bias between the values of local minima of the potential energy corresponding to the nematic and isotropic state. This bias is associated with lowering the temperature below that of the nematic-to-isotropic transition. The corresponding numerical results are shown in Fig. 5.

The first three figures essentially demonstrate the development of a configuration shown in Fig. 4. The larger thermodynamic forces driving the phase transition in the present case, however, push the size of the tactoid further down essentially to that of a vortex core. At this point, the isotropic region loses stability and splits into two vortices of degree $-1/2$, similar to what can be seen in Fig. 1.

5.3 Degree [math] tactoid

Next, we impose the constant boundary conditions

[TABLE]

and suppose that $\alpha=0$. The numerically computed evolution of a degree zero tactoid that results is shown in Fig. 6 and qualitatively resembles the behavior of a similar tactoid in the experiment as depicted in Fig. 2.

A similar shape is also seen in evolution of degree zero interfaces in the CSH-director model in [26] and is explained by the fact that the director has to be parallel to the interface. In particular, the interface cannot be smooth, for if it were, then it would carry a nonzero topological degree different from the degree on the boundary of the domain.

One interesting issue that we observed in the course of our simulations is illustrated in Fig. 7. When the undercooling is large ($\alpha=0.2$), the director appears to be orthogonal to the moving interface rather than being parallel to it as would be expected. A possible explanation for this effect is that the velocity of the interface is relatively large for larger undercoolings and the mobility of the director might not be sufficient for it to relax in a proper direction. We plan to investigate this behavior further in a future work.

5.4 Coalescence of Nematic Tactoids

Finally, the Landau - de Gennes model can also be used to simulate the reverse situation when positive nematic tactoids nucleate in the isotropic phase, then grow and coalesce to form the nematic phase with embedded topological defects (cf. [2]). In Fig. 8 the simulations were conducted subject to Neumann boundary data on $\partial\omega$ and assuming that $\omega$ has radius $1/2$, while $\alpha=0.01$. Three circular tactoids of different orientations were assumed to be present at the time $t=0$; in the course of the simulation, tactoids merged generating a single degree $-1/2$ defect. This situation closely resembles the original Kibble’s model of strings formed in early universe through coalescence of domains with different “phase” [29].

6 Acknowledgements

DG acknowledges the support from DMREF NSF

DMS-1729538. OL acknowledges the support from DMREF NSF DMS-1729509. PS acknowledges the support from the Simons Collaboration grant 585520.