1 Introduction↩︎

The structural characteristic of liquid crystals is the local anisotropy generated by the nonuniform orientational distribution of rigid, typically rod-like molecules. The local anisotropy has been expermimentally shown to bring about distinctive configurations both in stationary states and in dynamical evolutions, for which we refer to [1]–[5] and the references therein.

Theoretically, the dynamical models of liquid crystals are given by the evolutionary equations of order parameters that represent the local anisotropy, coupled with the Navier-Stokes equation for the velocity. These equations form an energy dissipative system. Different types of order parameters lead to models at different levels. For the most commonly observed uniaxial nematic phase, it is straightforward to use a unit vector as the order parameter. The corresponding dynamics of the unit vector field is the Ericksen-Leslie model [6], [7], which can be extended to incorporate flexibility in local anisotropy but is still limited to uniaxiality [8]–[10].

If one would like to fully characterize the local anisotropy and reach into molecular information, a natural choice is to consider the density function w.r.t. both the position and the orientation. The resulting dynamics is a kinetic equation, first proposed for spatially homogeneous cases in what is known as the Doi-Onsager theory [11] which has been studied extensively [12], [13]. Several attempts were later made to extend the Doi-Onsager theory to inhomogeneous cases [14]–[20]. For these models, the presence of orientational variables in combination with spatial variables brings significant computational challenges.

To balance the capability of describing local anisotropy with computational difficulties, tensor models are proposed. Tensor order parameters are able to distinguish different types of local anisotropy, so that tensor models can describe various phases and their dynamical behaviors. The interaction terms in the tensor models can also be derived from the microscopic level. Although the results obtained from tensor models of this type usually do not exactly reflect specific molecular architectures, they do reveal the changing trends of macroscopic phenomena as the molecular architecture varying, evidenced by several representative cases [21]–[23]. Based on the above features, the major points of focus in tensor models are the essential mathematical properties and overall structures, rather than the precision of specific terms. For rod-like molecules, the order parameter is chosen as a symmetric traceless tensor \(Q\). Dynamic tensor models can either be written down phenomenologically, such as Beris-Edwards [24] and Qian-Sheng models [25], or be derived from molecular models [21]. The latter type of models possesses clearer physical meanings in terms of molecular shape and interactions, and are suitable for molecules with more complex shapes [23]. However, a common problem for models of this type is that the evolutionary equations of the tensor involve higher-order tensors. It is necessary to supplement function relations between them and the tensor \(Q\) to close the system, which is called the closure approximation. Various closure approximations have been proposed, including some explicit functions [11], [26], [27], and the Bingham closure [28], [29] defined through the maximum entropy state. These closure approximations have been compared thoroughly [30], where the Bingham closure is found to perform better, which is attributed to the fact that the Bingham closure maintains the energy dissipation.

To comprehend how higher-order tensors affect the energy dissipation, we concentrate on the cases where the velocity is small. Under the reasonable approximation that the velocity is taken to be zero, the molecular models reduce to an evolutionary equation of the density function, while the tensor models reduce to an evolutionary equation of \(Q\). The equation of the density function can be written as a gradient flow, which ensures its energy dissipation. However, it is not always the case for the equation of the tensor \(Q\). It turns out that what underlies the fact that the Bingham closure maintains the energy dissipation is that it manages to preserve the gradient flow structure while other closure approximations do not. As we will briefly discuss in Sec. 3, when the Bingham closure is adopted, it introduces a singular entropy term in the free energy that constrains the eigenvalues of \(Q\) within \((-1/3,2/3)\). The dissipation operator is given by a fourth-order tensor that is negative definite, also guaranteed by the Bingham closure.

Although the Bingham closure maintains the gradient flow structure at the level of the tensor model, it defines the higher-order tensor as an implicit function through the maximum-entropy density function (also called the Bingham distribution) determined by \(Q\). It becomes a significant obstacle in computation, since implementation according to its definition results in high computational cost. To tackle this problem, fast algorithms for the computation of Bingham closure itself have been discussed in several works [31]–[37]. Nonetheless, they still cannot fill the gap towards solving the evolution equations efficiently, especially for spatially inhomogeneous cases.

Since it turns out that whether a tensor model works well relies heavily on the dissipation structure (but not necessarily on the specific form of the Bingham closure as we have explained above when introducing tensor models), it is worth seeking alternative simpler approaches for the closure approximation without breaking such a structure. In particular, it is better to avoid involving implicit functions through the density function and integrals on the unit sphere. Actually, this has successfully been done for the free energy of general liquid crystals (formed by rigid molecules with arbitrary shape) derived from molecular theory. The quasi-entropy, defined as the log-determinant of the covariance matrix, has been proposed for general rigid molecules to substitute the entropy term given by the maximum entropy density function (that is the Bingham distribution for rod-like molecules) [38]. It is shown that the quasi-entropy preserves the essential properties of the original entropy term: strict convexity; barrier function to constrain the covariance matrix positive definite; rotational invariance; consistency under symmetry reductions. The free energy with the quasi-entropy is able to capture the underlying physics correctly in representative cases. Specially, for rod-like molecules, the stationary points of the bulk energy are shown to be axisymmetric, which are further classified to describe the isotropic–uniaxial nematic phase transition. On the other hand, the fact that the quasi-entropy is an elementary function spontaneously reduces the computational complexity, so that preliminary discussions have been done on the design of efficient and accurate numerical methods [39], where the models do not contain higher-order tensors.

Because the free energy is a core ingredient in the dynamic model, we are now naturally in a position to discuss molecular-theory-based dynamic tensor models with the free energy including the quasi-entropy. In this paper, we propose the closure approximation using the quasi-entropy. The basic idea is similar to the treatment of the free energy, that is, to substitute the original entropy term with the quasi-entropy wherever it plays a role. Sec. 4 is dedicated to the tensor dynamics incorporating the quasi-entropy in the free energy and the closure approximation, which we outline below.

As we have mentioned, the Bingham closure is done by solving the maximum entropy density function first, then higher-order tensors are calculated. On the other hand, the original entropy term in the free energy is also calculated by the maximum entropy density function. But, such a formulation, with the density function as an intermediate, is difficult for us to carry out the quasi-entropy substitution. Thus, we first need to figure out an appropriate formulation for how the closure approximation is directly defined through the entropy term in the free energy. Actually, we are able to rewrite the Bingham closure as a constrained minimization problem. For the higher-order tensors appearing in the dynamic models, they can all be expressed by \(Q\) and a fourth-order symmetric traceless tensor \(Q_4\). We consider the maximum entropy density with both \(Q\) and \(Q_4\) fixed, to obtain a function of these two tensors. Then, the closure approximation is equivalently defined by minimizing this function of two tensors with \(Q\) fixed. Thus, the solution to this minimization problem expresses \(Q_4\) as a function of \(Q\). Using the formulation above, we substitute this function of \(Q\) and \(Q_4\) with the fourth-order quasi-entropy of \(Q\) and \(Q_4\), to establish the closure approximation by the quasi-entropy.

We shall show that the higher-order tensors obtained from the quasi-entropy closure approximation maintain several essential properties of the Bingham closure, due to the properties of the quasi-entropy. The strict convexity of the quasi-entropy guarantees the existence and uniqueness of the closure, with the higher-order tensor being a continuous function of \(Q\). The rotational invariance of the quasi-entropy leads to the fact that the higher-order tensor rotates together with the rotation of \(Q\). The fourth-order tensor satisfies the same symmetry as \(Q\) as a result of the consistency in symmetry reduction of the quasi-entropy. The domain of the quasi-entropy guarantees the negative definiteness of the dissipation operator. With these properties, the quasi-entropy closure approximation only requires to solve a convex minimization on an elementary function of three variables, together with tensor rotations.

In this way, we arrive at a tensor gradient flow derived from molecular theory with the quasi-entropy determining several core ingredients. The quasi-entropy gives a singular term in the free energy that is able to constrain the eigenvalues of \(Q\). This enables us to include in the free energy a cubic elastic term bounded from below (involving spatial derivatives). The gradient flow is indeed energy dissipative, which is a direct consequence of the negative definiteness of the dissipation operator given by the quasi-entropy closure approximation.

After the structures of the tensor gradient flow are sufficiently comprehended, we go on to consider numerical aspects for the gradient flow in Sec. 5, with emphases on preserving these structures. Since the quasi-entropy closure approximation is extremely easy to implement, we acquire enough flexibility to deal with other terms. To be as simple as possible, we propose a first-order-in-time scheme constraining the eigenvalues of \(Q\) within \((-1/3,2/3)\) and preserving energy dissipation unconditionally, which is also suitable for the cubic elastic term. We carry out numerical simulations in 2D, which, to our knowledge, has not been done before for tensor models requiring closure approximation. Evolutions of the isotropic–nematic interface and nematic defects are examined, where differences from gradient flows without incorporating higher-order tensors are evident.

Although many works on numerical simulations, both for molecular models and tensor models, have been done before, in this paper we do not attempt to compare with previous results. The main reason is that the problem settings and points of focus are quite different between this work and previous works (and also between those works themselves). To be specific, the previous works concentrate more on the effect of fluid velocity that is not discussed in this work. The interaction between rod-like molecules, deduced from the variational derivative of the free energy, does not have elastic terms but is usually given by a spatially dependent kernel function. Moreover, the systems that those works investigate are not actually energy dissipative because the boundary conditions they impose give rise to energy inputs. Based on the quasi-entropy closure approximation proposed in this work, we expect in the near future to extend it to the coupled system with the Navier-Stokes equations (i.e. not taking the velocity to be zero as an approximation). Comparisons of numerical results will be done in these forthcoming works. The quasi-entropy closure approximation for other rigid molecules is also a problem of interests. We shall give further concluding remarks on the above aspects in Sec. 6.

2 Preliminaries↩︎

For a system of rod-like molecules, the orientation of a single molecule is represented by a unit vector \(\mathbf{m}\in \mathbb{S}^2\), whose coordinates in the reference right-handed frame \((\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3)\) are denoted by \(m_i,\,i=1,2,3\). Denote by \(\mathcal{R}= \mathbf{m}\times \nabla_{\mathbf{m}}\) the rotational gradient operator on the unit sphere \(\mathbb{S}^2\). When written in the coordinates, it reads \(\mathcal{R}_{i} = \epsilon_{ijk}m_{j}{\partial}/{\partial m_k}\), where \(\epsilon\) is the Levi-Civita symbol, and summations on repeated indices are adopted throughout this work unless stated otherwise. The uniform unit measure on \(\mathbb{S}^2\) is denoted by \(\mathrm{d}\mathbf{m}\). For functions \(g_1(\mathbf{m})\), \(g_2(\mathbf{m})\) on \(\mathbb{S}^{2}\), the integration by parts holds as \(\int_{\mathbb{S}^2}g_1\mathcal{R} _{i}g_2\, \mathrm{d}\mathbf{m}= -\int_{\mathbb{S}^2}g_2\mathcal{R}_{i}g_1\mathrm{d}\mathbf{m}\).

We introduce some notations and basic results for tensors, largely following [40]. An \(n\)th-order tensor \(U\) is expressed in terms of the basis generated by \(\mathbf{e}_i\), \[\begin{align} U = U_{i_1i_2\dots i_n} \mathbf{e}_{i_1} \otimes \cdots \otimes \mathbf{e}_{i_n}, \quad i_1, \dots, i_n \in \{1, 2, 3\}. \end{align}\] The dot product between two tensors of the same order is defined as the sum of the products of the corresponding components, \(U \cdot W = U_{i_1\dots i_n} W_{i_1\dots i_n} .\) In particular, the components of a tensor can be given by \(U_{i_1\dots i_n} = U \cdot (\mathbf{e}_{i_1} \otimes \dots \otimes\mathbf{e}_{i_n}).\) Denote the rotation \(\mathfrak{t}\in SO(3)\) on a vector \(\mathbf{q}\) as \(\mathfrak{t}\circ\mathbf{q}\). If expressed in coordinates, \(\mathfrak{t}\) is represented by an orthogonal matrix with determinant one and \(\mathfrak{t}\circ \mathbf{q}\) is exactly given by a conventional matrix-vector product. The rotation can also be imposed on the tensor \(U\) by rotating the vectors \(\mathbf{e}_i\) while keeping the components unchanged, i.e., \[\mathfrak t\circ U= U_{i_1\cdots i_n}\,(\mathfrak{t}\circ\mathbf{e}_{i_1})\otimes\cdots\otimes (\mathfrak{t}\circ\mathbf{e}_{i_n}).\] The dot product is rotationally invariant, \(U_1\cdot U_2 = (\mathfrak{t}\circ U_1)\cdot(\mathfrak{t}\circ U_2)\). The integral is invariant both under rotations and inversions, i.e., \(\int_{\mathbb{S}^2}g(\mathbf{m}) \, \mathrm{d}\mathbf{m}= \int_{\mathbb{S}^2} g(\mathfrak t\circ \mathbf{m}) \, \mathrm{d}\mathbf{m}=\int_{\mathbb{S}^2}g(-\mathbf{m})\, \mathrm{d}\mathbf{m}\). In other words, the integral is invariant under arbitrary orthogonal transformations.

A tensor \(U\) is symmetric if \(U_{i_1 \dots i_n} = U_{i_{\sigma(1)} \dots i_{\sigma(n)}}\) for any permutation \(\sigma\) of \(\{1, \dots, n\}\). For an \(n\)th-order tensor \(U\), we define its average over all index permutations as \[{\left(U_{\mathrm{sym}}\right)}_{i_1i_2\dots i_n} = \frac{1}{n!} \sum_{\sigma} U_{i_{\sigma(1)} i_{\sigma(2)} \dots i_{\sigma(n)}},\] which is a symmetric tensor. We introduce the monomial notation for symmetric tensors generated by basis vectors of an orthonormal frame \((\mathbf{n}_1,\mathbf{n}_2,\mathbf{n}_3)\), \[\begin{align} \mathbf{n}_1^{k_1} \mathbf{n}_2^{k_2} \mathbf{n}_3^{k_3} = {\left( \underbrace{\mathbf{n}_1 \otimes \cdots \otimes \mathbf{n}_1}_{k_1} \otimes \underbrace{\mathbf{n}_2 \otimes \cdots \otimes \mathbf{n}_2}_{k_2} \otimes \underbrace{\mathbf{n}_3 \otimes \cdots \otimes \mathbf{n}_3}_{k_3} \right)}_{\mathrm{sym}}. \end{align}\] In this way, a homogeneous polynomial of \(\mathbf{n}_1,\mathbf{n}_2,\mathbf{n}_3\) represents a symmetric tensor. Since the second-order identity tensor \(\mathbf{i}\) satisfies \(\mathbf{i}= \mathbf{n}_1^2 + \mathbf{n}_2^2 + \mathbf{n}_3^2,\) we can define \(\mathbf{n}_1^{k_1} \mathbf{n}_2^{k_2} \mathbf{n}_3^{k_3}\mathbf{i}^l\) in the same manner.

For an \(n\)th-order symmetric tensor \(W\), its trace is defined as an \((n-2)\)th-order tensor \({\operatorname{tr}(W)}_{i_1 \dots i_{n-2}}=W_{i_1 \dots i_{n-2} j j}\). If \(\operatorname{tr} W\) is the zero tensor, \(W\) is called symmetric traceless. For any \(n\)th-order symmetric tensor \(U\), there exists a unique symmetric traceless tensor \({(U)}_0=U-{(\mathbf{i}\otimes W)}_{\mathrm{sym}}\) (see, for example, Proposition 3.2 in [40]). We call \({(U)}_0\) the symmetric traceless tensor generated by \(U\). For example, \((\mathbf{m}^2)_0=\mathbf{m}^2-{\mathbf{i}}/{3}\) is the symmetric traceless tensor generated by \(\mathbf{m}^2\). The \(n\)th-order symmetric traceless tensors form a space of dimension \(2n+1\). For the second-order symmetric traceless tensors, an orthogonal basis can be chosen as \[\label{eq:the32basis32of32symmetric32traceless32tensors} \mathbf{s}_1 = \mathbf{n}_1^2 - \frac{\mathbf{i}}{3},\,\mathbf{s}_2 = \mathbf{n}_2^2 - \mathbf{n}_3^2,\,\mathbf{s}_3 = \mathbf{n}_1\mathbf{n}_2,\,\mathbf{s}_4 = \mathbf{n}_1\mathbf{n}_3,\,\mathbf{s}_5 = \mathbf{n}_2\mathbf{n}_3.\tag{1}\]

The orientation distribution at the position \(\mathbf{x}\) is described by a density function \(f\left(\mathbf{x},\mathbf{m}\right)>0\) satisfying the normalization condition \(\int_{\mathbb{S}^2} f\,\mathrm{d}\mathbf{m}= 1\). In tensor models, the order parameter is chosen as the symmetric traceless tensor generated by the second moment, defined as \(Q(\mathbf{x}) = \langle (\mathbf{m}^2)_0\rangle\), where the notation \(\langle \cdot \rangle\) represents averaging over \(\mathbb{S}^2\) w.r.t. the density \(f(\mathbf{x},\mathbf{m})\). When necessary, we use the notation \(\langle\cdot\rangle_{f}\) to specify over which density function \(f\) the average is taken. It is also natural to regard a second-order tensor as a \(3\times 3\) matrix, so that notions for a matrix, such as eigenvalues and eigenvectors, can be adopted. The definition of \(Q\) implies that the eigenvalues \(\lambda(Q)\) lie within the open interval \(\left(-{1}/{3}, {2}/{3}\right)\). We denote by \(\mathscr{Q}_{\mathrm{phys}}\) the second-order symmetric traceless tensors satisfying the above eigenvalue constraints, \[\begin{align} \mathscr{Q}_{\mathrm{phys}}= \left \{ Q : Q_{ij} = Q_{ji}, \, Q_{ii} = 0, \, \lambda(Q) \in \left(-{1}/{3}, {2}/{3}\right) \right \} . \end{align}\]

3 Tensor dynamics from kinetic equation with the Bingham closure↩︎

Our starting point is the kinetic equation of the density function \(f(\mathbf{x},\mathbf{m},t)\), which includes both spatial and orientational convection and diffusion terms. The spatial diffusion term is often omitted because it can be easily recognized as a higher-order infinitesimal under rescaling. As a result, the equation is written in the following form, \[\label{eq:the32molecular32model32of32Smoluchowski32equation} \frac{\partial f}{\partial t}+\nabla\cdot(\mathbf{v}f)=\mathcal{R}\cdot\bigl(\mathcal{R}f+f\mathcal{R}\mu_{\mathrm{r}}\bigr)+\mathcal{R}\cdot\bigl((\mathbf{m}\cdot\nabla)\mathbf{v}\times\mathbf{m}f\bigr),\tag{2}\] where \(\mathbf{v}(\mathbf{x},t)\) is the fluid velocity, \(\mu_{\mathrm{r}}=\delta F_{\mathrm{r}}/\delta f\) is the interaction potential given by the variational derivative of the interaction energy \(F_{\mathrm{r}}\). For most cases, the interaction energy \(F_{\mathrm{r}}\) is a functional of \(Q\). Consequently, the interaction potential \(\mu_{\mathrm{r}}\) can be expressed as \[\label{eq:mean-field32potential} \mu_{\mathrm{r}} = \frac{\delta F_{\mathrm{r}}}{\delta f} = \frac{\delta F_{r}}{\delta Q}\cdot\frac{\delta Q}{\delta f}=V_Q\cdot{\left(\mathbf{m}^2\right)}_0,\tag{3}\] where \(V_Q=\delta F_{\mathrm{r}}/\delta Q\) depends solely on \(Q\). The velocity \(\mathbf{v}\) is either prescribed (such as steady shear flows studied extensively [19], [20], [41]–[43]), or obeys a Navier-Stokes equation, which is not written down here since it is beyond the focus of this work.

The dynamic tensor model can be derived from [eq:the molecular model of Smoluchowski equation] by multiplying it with \({\left(\mathbf{m}^2\right)}_0\) and integrating on \(\mathbb{S}^2\), leading to an equation of \(Q\), \[\label{eq:general32equation32of32tensor32model} \begin{align} \frac{\partial Q_{ij}}{\partial t} + \mathbf{v}_k\partial_{k} Q_{ij}= & -6Q_{ij}-\bigl(2Q_{ik}{(V_Q)}_{kj}+2Q_{jk}{(V_Q)}_{ki}+\frac{4}{3}{(V_Q)}_{ij}-4\langle\mathbf{m}^4\rangle_{ijkl}{(V_Q)}_{kl}\bigr) \\ & +Q_{ik}\partial_{k}\mathbf{v}_j +Q_{jk}\partial_{k}\mathbf{v}_i +\frac{1}{3}\bigl(\partial_i\mathbf{v}_{j}+\partial_j\mathbf{v}_{i}\bigr)-2\langle\mathbf{m}^4\rangle_{ijkl}\partial_{k}\mathbf{v}_l. \end{align}\tag{4}\] Notably, the fourth-order moment \(\langle\mathbf{m}^4\rangle\) appears on the right-hand side. It is then necessary to supplement a closure approximation, i.e. to express the fourth moment as a function of \(Q\). It is worth pointing out that the fourth moment also appears in the Navier-Stokes equation, if it is included as part of the whole system.

It is thus clear that the closure approximation is an essential constituent of the tensor model deduced from the kinetic equation. As we have mentioned, various closure approximations have been proposed and the Bingham closure performs better, which is deemed owing to the maintenance of dissipation structure. Let us figure out below the dissipation structure and other significant properties resulting from the Bingham closure.

Since the dissipation structure mainly comes from the diffusion term, we consider in what follows the case of small velocity and set \(\mathbf{v}=0\) that eliminates the convection terms.

Under this assumption, only the diffusion term remains in [eq:general equation of tensor model], \[\label{eq:the32molecular32model32of32Smoluchowski32equation32with32zero32velocity} \frac{\partial f}{\partial t}=\mathcal{R}\cdot\bigl(\mathcal{R}f+f\mathcal{R}\mu_{\mathrm{r}}\bigr).\tag{5}\] Noticing that \(\mathcal{R}f=f\mathcal{R}(1+\log \! f)\), we rewrite \(\mathcal{R}f+f\mathcal{R}\mu_{\mathrm{r}}=f\mathcal{R}\mu\), where \(\mu=\delta F/\delta f\) is the variational derivative of the free energy \(F\), given by the sum of the interaction energy \(F_{\mathrm{r}}\) and the entropy term, \[\label{eq:free32energy32of32molecular} F=\int \, f\log f\,\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{m}+F_{\mathrm{r}}[f].\tag{6}\] When \(f>0\), [eq:the molecular model of Smoluchowski equation with zero velocity] is a gradient flow of the energy \(F\) with the dissipation operator \(\mathcal{R}\cdot\bigl(f\mathcal{R}(\cdot)\bigr)\). Therefore, assuming that the boundary terms vanish in integration by parts, we deduce the energy dissipation law, \[\frac{\mathrm{d}F}{\mathrm{d}t} = -\int f \left|\mathcal{R}\mu\right|^2\,\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{m}\leq 0.\]

Now let us turn to the equation 4 of the tensor \(Q\) and assume \(\mathbf{v}=0\). Define a fourth-order tensor \({M}\), \[\label{eq:the32structure32of32dissipative32operator32M} \begin{align} {M}_{iji'j'} & =\langle\mathcal{R}_k{(\mathbf{m}^2)}_0\otimes\mathcal{R}_k{(\mathbf{m}^2)}_0\rangle_{iji'j'} = 4\mathcal{H}{\left(\langle\mathbf{m}^2\otimes(\mathbf{i}-\mathbf{m}^2)\rangle\right)}_{iji'j'} \\ & = -4\langle\mathbf{m}^4\rangle_{iji'j'} + 4\mathcal{H}{\left(Q\otimes \mathbf{i}\right)}_{iji'j'} + \frac{4}{3}\mathcal{H}{\left(\mathbf{i}\otimes\mathbf{i}\right)}_{iji'j'}, \end{align}\tag{7}\] where \({\mathcal{H}(U)}_{ijkl} = (U_{ikjl}+U_{iljk}+U_{jkil}+U_{jlik})/4\). For a second-order tensor \(W\), denote in short by \(MW\) the second-order tensor given by \((MW)_{ij}=M_{iji'j'}W_{i'j'}\). The equation of the tensor then becomes \[\label{eq:tensor32model32with32Bingham32closure} \frac{\partial Q}{\partial t} = -6Q - {M}V_Q.\tag{8}\] The above form itself does not imply that the equation is a gradient flow.

When the Bingham closure is adopted, [eq:tensor model with Bingham closure] can be rewritten in an equivalent form with the gradient flow structure, which we explain below. The Bingham closure is defined by solving the maximum entropy state from \(Q\), i.e., \[\label{eq:the32Bingham32closure} \min_{f} \; \int_{\mathbb{S}^2} f(\mathbf{m})\log f(\mathbf{m})\,\mathrm{d}\mathbf{m} \qquad \text{s. t.} \qquad \int_{\mathbb{S}^2} {(\mathbf{m}^2)}_0 f(\mathbf{m})\,\mathrm{d}\mathbf{m}= Q, \quad \int_{\mathbb{S}^2} f(\mathbf{m})\,\mathrm{d}\mathbf{m}= 1.\tag{9}\] The maximum entropy state is sometimes referred to as the Bingham distribution, denoted by \(f_Q\). The existence and uniqueness of \(f_Q\) have been established (see, e.g., [22], [40], [44], [45]). Together with symmetry properties of \(f_Q\), the results are summarized in the following.

The fourth-order tensor \({M}\) is then calculated from \(f_Q\). Since \(f_Q\) is uniquely determined by \(Q\), the tensor \({M}\) is expressed as a function of \(Q\). We denote by \({M}^{\mathrm{Bin}}\) the tensor \({M}\) determined from the Bingham closure.

Armed with the density \(f_Q\), we define its entropy \({\zeta}\), which is a function of \(Q\), \[\label{eq:the32free32energy32of32the32tensor32model32with32Bingham32closure} {\zeta}(Q) = \int_{\mathbb{S}^2} f_Q(\mathbf{m})\log f_Q(\mathbf{m})\mathrm{d}\mathbf{m}= B_Q \cdot Q - \log Z_Q.\tag{13}\]

For the closure approximation to be successfully carried out, one requires \(\lambda(Q)\in (-1/3,2/3)\) in the equation of \(Q\). This can actually be constrained by the term \({\zeta}(Q)\) in the free energy, which is a barrier function. The properties of \({\zeta}(Q)\) are also discussed in the literature (cf. [38], [44], [45]), which are summarized below.

The fourth-order tensor \({M}^{\mathrm{Bin}}\) also satisfies a few symmetry properties resulting from those of \(f_Q\).

To summarize, the Bingham closure is carried out by solving the maximum entropy state \(f_Q\) from \(Q\), followed by calculating higher-order tensors using \(f_Q\). The gradient flow structure is recognized by rewriting \(6Q={M} ^{\mathrm{Bin}}(\partial \zeta/\partial Q)\). The singular term \(\zeta(Q)\) constrains \(\lambda(Q)\in (-1/3,2/3)\), so that \(f_Q\) can be uniquely solved. Furthermore, \({M}^{\mathrm{Bin}}\) is positive definite as a result of \(f_Q>0\), and enjoys several symmetry properties inherited from \(f_Q\).

On the other hand, it is important that these structures and properties are maintained in computations, which brings about difficulties. If we carry out the Bingham closure through \(f_Q\), we need to handle implicit functions involving integrals on \(\mathbb{S}^2\). Naive implementations of these functions lead to high computational cost [16], [34]. To this end, several fast algorithms have been developed for the Bingham closure using various approximation formulae [32]–[34], [36], [37] or neural networks if high accuracy is not needed [31], [35], [46], [47], including the mappings between \(Q\) and \(B_Q\) [31], [36], as well as the computation of \(\zeta(Q)\), \(Z_Q\) [34] and the fourth-order tensor \({M}^{\mathrm{Bin}}\) from \(Q\) or \(B_Q\) [32], [33], [37], [48]. However, there is still a large gap towards efficient numerical methods that preserve the desired structures and properties. One may choose to change the independent variables to \(B_Q\) from which \(Q\) and \({M}^{\mathrm{Bin}}\) can be calculated [32]. This gives rise to much more complicated terms for the functional derivatives of the interaction energy \(F_{\mathrm{r}}\). An alternative way is to utilize the form 8 that only requires to compute \({M}^{\mathrm{Bin}}\) from \(Q\). But it literally disregards the gradient flow structure and the eigenvalue constraints, so that they are difficult to be guaranteed in the schemes. If schemes are built based on ?? with fast algorithms from \(Q\) to \(B_Q\), the main obstacle lies in the singularities for \(\lambda(Q)\) close to the boundary of \(\mathscr{Q}_{\mathrm{phys}}\), which restricts sufficient accuracy, and in turn impair the convexity of \(\zeta(Q)\) up to certain precision required for the solvability of the schemes.

The implementations of the Bingham closure bring limitations in the numerical simulations for tensor dynamics derived from the molecular models. Indeed, numerical simulations have been carried out only for 1D, or 2D with the orientation restrained on the unit circle \(\mathbb{S}^1\) [16], [20]. The difficulties become worse if rigid molecules of complex shapes are considered, where more order parameters and higher-order tensors are present [49].

As we have mentioned at the beginning of this article, for tensor models the emphasis is the essential properties and structures, instead of the precisions of specific terms, especially the Bingham closure on which many works have put efforts. This can also be recognized from the interaction energy \(F_{\mathrm{r}}\), usually given by a few terms in the gradient expansion that are only able to roughly depict the molecular information. Thus, given that the essential properties and structures are maintained, it is desirable if we are able to seek alternative ways to treat the terms involving \(f_Q\) (including \(B_Q\), \(\zeta(Q)\), and \({M}\)).

Actually, in the free energy, the quasi-entropy has been proposed as a substitute for \(\zeta(Q)\). In the following section, we introduce the quasi-entropy and discuss how to use it in the closure approximation, so that we construct a tensor dynamics based on the quasi-entropy.

4 Tensor dynamics with the quasi-entropy↩︎

The quasi-entropy [38] is a class of elementary functions of tensors averaged on \(SO(3)\). It is proposed for general choice of tensor order parameters, to act as the entropy term in the free energy. Generally, the quasi-entropy is defined as the log-determinant of the covariance matrix for the tensors up to certain even order. It is shown that the quasi-entropy possesses the essential properties of the entropy defined from the maximum entropy state, including strict convexity, enforcing positive definiteness of the covariance matrix, invariance under rotations, and consistency under symmetry reduction. Moreover, for various shapes of molecules, the free energy incorporating the quasi-entropy is able to capture the underlying physics correctly. We shall see shortly that the above properties are also crucial for the closure approximation.

Below, we start from briefing the results for the case of one tensor \(Q\). Then, we discuss the closure approximation using the quasi-entropy, which is built on a reformulation of the closure approximation as a minimization problem w.r.t. a fourth-order tensor.

4.1 Quasi-entropy in the free energy↩︎

Recall that for the free energy given in ?? , the entropy term \(\zeta(Q)\) is defined in 13 . For the tensor \(Q\), we use the second-order quasi-entropy \(\Xi_2(Q)\) to substitute \(\zeta(Q)\), which reads \[\begin{align} \label{eq:second-order32quasi-entropy} \Xi_2(Q)=-\log\det\left(Q+\frac{\mathbf{i}}{3}\right)-2\log\det\left(\frac{\mathbf{i}}{3}-\frac{Q}{2}\right). \end{align}\tag{15}\] The domain of \(\Xi_2\) consists of all \(Q\) such that \(Q+\mathbf{i}/3\) and \(\mathbf{i}/3-Q/2\) are positive definite, which is exactly \(\mathscr{Q}_{\mathrm{phys}}\). The quasi-entropy \(\Xi_2\) possesses the essential properties of \(\zeta\).

For the proof, see Proposition 2.1 in [39] (also Theorem 4.8 in [38] for general cases).

To further illustrate the feature of \(\Xi_2\), we consider the bulk energy containing a quadratic term, \[\label{eq:bulk32energy} \Xi_2(Q)-\frac{1}{2}c_{02}\lvert Q\rvert^2.\tag{16}\] When \(\zeta(Q)\) takes the place of \(\Xi_2(Q)\), we recover exactly the Maier-Saupe energy. It is known that the stationary points of the Maier-Saupe energy are uniaxial. The same results hold when using the quasi-entropy.

Based on the uniaxial expression, the stationary points have been completely classified [38], for which two critical values of \(c_{02}\) describe the first-order isotropic–nematic phase transition. When the elastic energy (terms with spatial derivatives) is considered, the defect morphology has also been validated in an \(L^2\) gradient flow without including the higher-order tensors [39]. Since the free energy is a core ingredient of the dynamic model, it is natural to discuss the treatment of higher-order tensors using quasi-entropy in a consistent way with the free energy.

4.2 Closure approximation by the quasi-entropy↩︎

It is clear from the discussions in the previous section that the closure approximation is also determined through entropy. However, it is done by solving the maximum entropy state, making a substitution like the free energy not quite evident. In what follows, we reformulate the Bingham closure as a minimization problem w.r.t. an entropy function of \(Q\) and \({M}\).

By far, we have not gone into a basic problem: the independent variables of the tensor \({M}\), which is unimportant with a density function in hand. Actually, the independent variables of a tensor can be given by symmetric traceless tensors [40]. For \({M}\), defined in 7 , we express it by \(Q\) and a fourth-order symmetric traceless tensor \(Q_4\), \[\label{eq:the32dissipative32operator32M} {M}= -4Q_4 - \frac{4}{7}\mathcal{A}(Q) - \frac{2}{15}{E},\tag{17}\] where \[\label{eq:symtrls-M} \begin{align} & Q_4 = \langle {(\mathbf{m}^4)}_0 \rangle = \langle \mathbf{m}^4 - \frac{6}{7} \mathbf{m}^2 \mathbf{i}+ \frac{3}{35} \mathbf{i}^2 \rangle, \\ & {E}_{ijkl} = 2\delta_{ij} \delta_{kl} - 3\delta_{ik} \delta_{jl} - 3\delta_{il} \delta_{jk}, \\ & \mathcal{A}{(Q)}_{ijkl} = Q_{ij} \delta_{kl} + Q_{kl} \delta_{ij} - \frac{3}{4} (Q_{ik} \delta_{jl} + Q_{il} \delta_{jk} + Q_{jk} \delta_{il} + Q_{jl} \delta_{ik}). \end{align}\tag{18}\] We are now ready to reformulate the Bingham closure. Define \[\label{eq:the32equvialent32form32of32Bingham32closure32Step321} \begin{align} & \widetilde{\zeta}(Q,Q_4)=\min \int_{{\mathbb{S}}^2} f \log f\,\mathrm{d}\mathbf{m}\quad \text{s. t. } \quad \int_{{\mathbb{S}}^2} f \,\mathrm{d}\mathbf{m}=1, \quad \langle {\left(\mathbf{m}^2\right)}_0 \rangle = Q, \quad \langle{\bigl(\mathbf{m}^4\bigr)}_0 \rangle = Q_4. \end{align}\tag{19}\]

The convenience of the formulation 19 is that there is an entropy function \(\widetilde{\zeta}(Q,Q_4)\). In a similar way of dealing with the free energy, we now substitute \(\widetilde{\zeta}(Q,Q_4)\) with the quasi-entropy. Since the entropy function depends on the fourth-order tensor \(Q_4\), we use the fourth-order quasi-entropy, denoted by \(\Xi_4\), in the substitution. The function \(\Xi_4\) can be found in previous works for cases with more tensors for other symmetries [38], [50]. A more straightforward version is given by [38], for which we set \(M^4_1=0\) by symmetry reduction and arrive at the expressions below.

Let \(U^{[j]}\) be a \(j\)th-order tensor and \((\mathbf{n}_1,\mathbf{n}_2,\mathbf{n}_3)\in SO(3)\), and define vectors and matrices as follows, \[\begin{align} &\Phi_{2}{\left(U^{[2]}\right)}_{i}=U^{[2]}\cdot\mathbf{s}_i, \quad \Psi_{2}{\left(U^{[2]}\right)}_{ij}=U^{[2]}\cdot\mathbf{n}_i\otimes\mathbf{n}_j, \\ &\Psi_{3}{\left(U^{[3]}\right)}_{ij}=U^{[3]}\cdot\mathbf{n}_i\otimes\mathbf{s}_j, \quad \Psi_{4}{\left(U^{[4]}\right)}_{ij}=U^{[4]}\cdot\mathbf{s}_i\otimes\mathbf{s}_j, \end{align}\] where \(\mathbf{s}_{j},j=1,\dots,5\), are defined in 1 . Define \[\label{eq:the32three32matrices32of32Psi954} \begin{align} & C_1=Q_4-\frac{4}{21}\mathcal{A}(Q)-\frac{1}{45}{E},\; C_2=\frac{1}{8}Q_4+\frac{1}{7}\mathcal{A}(Q)-\frac{1}{60}{E}, \; & C_3=-\frac{1}{2}Q_4-\frac{1}{14}\mathcal{A}(Q)-\frac{1}{60}{E}, \end{align}\tag{20}\] and the fourth-order quasi-entropy \(\Xi_4\) is then written as \[\label{Simplified32Xi954} \begin{align} \Xi_4\left(Q,\,Q_4\right)=\, & -\log\det\Psi_2\left(Q+\mathbf{i}/3\right)-2\log\det \Psi_4(C_2) -\log\det(\Psi_4(C_1)-\Phi_2(Q)^t\Phi_2(Q)) \\ & -2\log\det\begin{pmatrix} \Psi_2\left(\mathbf{i}/3-Q/2\right) & \frac{1}{4}\Psi_3\left(\mathcal{B}\left(Q\right)\right) \\ \frac{1}{4}\Psi_3\left(\mathcal{B}\left(Q\right)\right)^t & \Psi_4(C_3) \end{pmatrix}, \end{align}\tag{21}\] where \(\mathcal{B}{\left(Q\right)}_{ijk}=\epsilon_{ijs}Q_{ks}+\epsilon_{iks}Q_{js}\) and recall that \(\epsilon\) is the Levi-Civita symbol. The domain of \(\Xi_4\) consists of \((Q,Q_4)\) making each the matrix after a log-determinant in 21 are positive definite, enforced by the barrier function property given by the log-determinant. The domain contains all the \((Q,Q_4)\) calculated from the average under some density function, since the density function ensures that the covariance matrix is positive definite. Let us restate the properties of \(\Xi_4\), which were discussed for general cases in [38], Proposition 4.1 and Theorem 4.6.

A convenient choice of \((\mathbf{n}_1,\mathbf{n}_2,\mathbf{n}_3)\) would be the eigenframe of \(Q\). The closure approximation by the quasi-entropy is given by substituting \(\widetilde{\zeta}(Q,Q_4)\) with \(\Xi_4(Q,Q_4)\) in the minimization problem in ?? , which reads \[\label{eq:the32quasi-entropy32closure} Q_4={\arg\min}\,\Xi_4\left(Q,Q_4\right),\quad Q \in \mathscr{Q}_{\mathrm{phys}}\;\text{fixed}.\tag{22}\] The properties of the \(Q_4\) obtained by the above quasi-entropy closure approximation result from the properties of the \(\Xi_4\).

In addition, the rotational invariance of \(\Xi_4\) leads to the fact that \(Q_4\) rotates together with \(Q\).

The fourth-order tensor \({M}\) is defined by 17 with \(Q_4\) obtained from 22 , denoted as \({M}^{\mathrm{qent}}\). Its existence, uniqueness and continuity follow from those of \(Q_4\). Moreover, \({M}^{\mathrm{qent}}\) has the same form as \({M}^{\mathrm{Bin}}\) when expressed by the eigenframe of \(Q\).

In the end, the domain of the quasi-entropy \(\Xi_4\) guarantees the positive definiteness of \({M}^{\mathrm{qent}}\).

When \((\mathbf{n}_1,\mathbf{n}_2,\mathbf{n}_3)\in SO(3)\) takes an eigenframe of \(Q\), let us write down the quasi-entropy \(\Xi_4\). We calculate \(\Psi_4(E)\), \(\Psi_4\big(\mathcal{A}(\mathbf{n}_1^2)_0\big)\), \(\Psi_4\big(\mathcal{A}(\mathbf{n}_2^2-\mathbf{n}_3^2)\big)\), \(\Psi_4\big((\mathbf{n}_1^4)_0\big)\), \(\Psi_4\big((\mathbf{n}_2^4)_0\big)\), \(\Psi_4\big((\mathbf{n}_1^2\mathbf{n}_2^2)_0\big)\), yielding \[\label{eq:the32coordinate32matrices32of326tensors} \begin{align} & X_1=\operatorname{diag}\left(\begin{pmatrix} -9 & 0 \\ 0 & -3 \end{pmatrix},-12,-12,-12\right),\, X_2=\operatorname{diag}\left(\begin{pmatrix} -\frac{3}{2} & 0 \\ 0 & \frac{1}{2} \end{pmatrix},-1,-1,2\right) \\ & X_3=\operatorname{diag}\left(\begin{pmatrix} 0 & \frac{3}{2} \\ \frac{3}{2} & 0 \end{pmatrix},-3,3,0\right),\, X_4=\operatorname{diag}\left(\begin{pmatrix} \frac{18}{35} & 0 \\ 0 & \frac{1}{35} \end{pmatrix},-\frac{16}{35},-\frac{16}{35},\frac{4}{35}\right) \\ & X_5=\operatorname{diag}\left(\begin{pmatrix} \frac{27}{140} & -\frac{3}{28} \\ -\frac{3}{28} & \frac{19}{140} \end{pmatrix},-\frac{16}{35},\frac{4}{35},-\frac{16}{35}\right),\, X_6=\operatorname{diag}\left(\begin{pmatrix} -\frac{9}{35} & \frac{3}{28} \\ \frac{3}{28} & -\frac{1}{70} \end{pmatrix},\frac{18}{35},-\frac{2}{35},-\frac{2}{35}\right). \end{align}\tag{25}\] Define \[\begin{align} & \omega={\left(s,b,0,0,0\right)}^T,\quad W_1\overset{\text{def}}{=}\Psi_2\left(Q+\frac{1}{3}\mathbf{i}\right)=\operatorname{diag}\left(\frac{1}{3}\left(2s+1\right),\frac{1}{3}\left(1-s\right)+b,\frac{1}{3}\left(1-s\right)-b\right), \\ & W_2=W_3\overset{\text{def}}{=}\Psi_2\left(\frac{1}{3}\mathbf{i}-\frac{1}{2}Q\right)=\operatorname{diag}\left(\frac{1}{3}\left(1-s\right),\frac{1}{6}\left(2+s\right)-\frac{1}{2}b,\frac{1}{6}\left(2+s\right)+\frac{1}{2}b\right), \\ & T_2=\begin{pmatrix} 0 & 0 & 0 & 0 & -b \\ 0 & 0 & 0 & \frac{1}{2}\left(s+b\right) & 0 \\ 0 & 0 & \frac{1}{2}\left(-s+b\right) & 0 & 0 \end{pmatrix}. \end{align}\] Then, the quasi-entropy \(\Xi_4\) is written as \[\label{Xi432ss} \begin{align} \Xi_4\left(s,b,a_1,a_2,a_3\right) = & -\log\det\left(-\frac{1}{45}X_1-\frac{4}{21}\left(sX_2+bX_3\right)+a_1X_4+a_2X_5+a_3X_6 -\omega\omega^T\right) \\ & -2\log\det \left(-\frac{1}{15}X_1+\frac{4}{7}\left(sX_2+bX_3\right)+\frac{1}{2}\left(a_1X_4+a_2X_5+a_3X_6\right)\right) \\ & -2\log\det\begin{pmatrix} W_2 & T_2 \\ T_2^T & -\frac{1}{60}X_1-\frac{1}{14}\left(sX_2+bX_3\right)-\frac{1}{2}\left(a_1X_4+a_2X_5+a_3X_6\right) \end{pmatrix} \\ & -\log\det W_1 + c_0, \end{align}\tag{26}\] for some constant \(c_0\). In the first two lines, each matrix is \(5\times 5\) consisting of diagonal blocks of \(2\times 2\) and three \(1\times 1\). In the third line, the matrix is \(8\times 8\) and can be rearranged so that it consists of four \(2\times 2\) diagonal blocks.

By Theorem 10 the closure approximation can be done according to the following three steps.

1. Given \(Q\), find the \(s\), \(b\), and the eigenframe \((\mathbf{n}_1,\mathbf{n}_2,\mathbf{n}_3)\in SO(3)\).

2. Solve \((a_1,a_2,a_3)\) by minimizing \(\Xi_4(s,b,a_1,a_2,a_3)\) given in 26 with \((s,b)\) fixed.

3. Assemble \(Q_4\) and \({M}^{\mathrm{qent}}\) using 23 and 17 .

In the second step, we only need to minimize an elementary function, which is easy to implement. In particular, the block diagonal structure of 26 makes it convenient to evaluate the function value and its derivatives and to determine whether certain \((s,b,a_1,a_2,a_3)\) lies within the domain of \(\Xi_4\).

4.3 Gradient flow with the quasi-entropy↩︎

In ?? , we use \(\Xi_2(Q)\) to substitute \(\zeta(Q)\) in the free energy and let the tensor \({M}\) be given by \({M}^{\mathrm{qent}}\). In this way, we obtain the following gradient flow, \[\label{eq:DTM32with32quasi-entropy} \begin{align} \frac{\partial Q}{\partial t} & = -{M}^{\mathrm{qent}} \mu_Q,\quad \mu_Q=\frac{\delta F}{\delta Q}. \end{align}\tag{27}\] The free energy is given by \[\label{eq:free32energy32of32the32DTM32with32quasi-entropy} F[Q]=\int \Xi_2\left(Q\right)\mathrm{d}\mathbf{x}+F_{\mathrm{r}}[Q],\tag{28}\] where the interaction energy \(F_r\) contains the bulk and the elastic energy (involving spatial derivatives), \[\label{eq:interaction32energy} F_{\mathrm{r}}=\frac{1}{2}\int -c_{02}\left|Q\right|^2+c_{21}\left|\nabla Q\right|^2+c_{22}\partial_{i}Q_{ik}\partial_{j}Q_{jk}+c_{24}Q_{ij}\partial_{i}Q_{kl}\partial_{j}Q_{kl}\,\mathrm{d}\mathbf{x},\tag{29}\] with the short notations \(|Q|^2=Q\cdot Q\), \(|\nabla Q|^2=\partial_iQ_{jk}\partial_iQ_{jk}\). We assume that, after integration by parts, the boundary terms either vanish or contribute only a constant. This assumption can be satisfied under the periodic or Dirichlet boundary conditions, for which we omit the detailed discussions. It is worth noting that the second-order quasi-entropy \(\Xi_2\) gives a singular term in the free energy that is able to constrain the eigenvalues of \(Q\). This enables us to include in the free energy a cubic elastic term while the free energy is bounded from below. Without this constraint, it has been shown that this cubic term is unbounded from below [44].

The energy dissipation law follows from the fact that \({M}^{\mathrm{qent}}\) is positive definite.

5 Numerical simulations↩︎

As we have discussed above, in tensor dynamics (no matter with the Bingham closure or the quasi-entropy) the eigenvalue constraints and the dissipation structure play key roles. Maintaining these properties thus becomes the major target computationally. When adopting the quasi-entropy, we are able to construct numerical schemes in a natural way to take care of these properties. To illustrate this convenience, we construct a temporal first-order scheme for 27 satisfying these properties. In computations, we shall pay attention to these properties and the effects of the fourth-order tensor \({M}\) and the \(c_{24}\) term.

5.1 Numerical methods↩︎

The most significant point in a numerical scheme for 27 is to keep \(Q\in \mathscr{Q}_{\mathrm{phys}}\), because it is necessary for the free energy and the closure approximation to be defined. In the gradient flow, it is achieved by the barrier given by the quasi-entropy \(\Xi_2\). Therefore, we choose to discretize \(\partial \Xi_2/\partial Q\) implicitly. Because the fourth-order tensor \({M}^{\mathrm{qent}}\) requires solving a minimization problem, it is preferable to deal with it explicitly in the time discretization. In this way, the closure approximation is decoupled from the solving procedure of the scheme, so that its easy implementation described in the previous section is inherited in the scheme.

Moreover, since the fourth-order tensor \({M}^{\mathrm{qent}}\) from the explicit discretization is positive definite, it suffices to find an appropriate way to discretize \(\mu_Q=\delta F/\delta Q\) to maintain the dissipation structure. Here, we choose the convex splitting method since it ensures the existence and uniqueness of the scheme. Consider the splitting of the free energy \(F[Q] = F_{+}[Q] - F_{-}[Q]\) in the following form, \[\label{eq:the32convex32splitting} \begin{align} F_{+}[Q] = & \int \Xi_2(Q) +\frac{1}{2}\Bigl( \gamma_1\left|\nabla Q\right|^2 +c_{22}\partial_i Q_{ik}\partial_j Q_{jk} \Bigr) \\ & \quad +\frac{c_{24}}{2} Q_{ij}\partial_i Q_{kl}\partial_j Q_{kl} +\frac{\gamma_2}{2}(\left|Q\right|^2 +\frac{1}{2}\left|\nabla Q\right|^4) +\frac{\gamma_3}{2}\left|\nabla Q\right|^2 \,\mathrm{d}\mathbf{x}\\ F_{-}[Q] = & \frac{1}{2}\int (c_{02}+\gamma_2)\left|Q\right|^2 +\frac{\gamma_2}{2}\left|\nabla Q\right|^4 +(\gamma_1+\gamma_3-c_{21})\left|\nabla Q\right|^2 \,\mathrm{d}\mathbf{x}, \end{align}\tag{32}\] where \(\gamma_i,\,(i=1,2,3)\) are to be determined to make \(F_{\pm}[Q]\) convex w.r.t. \(Q\). Denote \(\mu_{\pm}=\delta F_{\pm}/\delta Q\). The scheme is given by \[\label{eq:first32order32scheme32derivation} \frac{Q^{n+1} - Q^{n}}{\delta t} = - {M}^{\mathrm{qent},n}\mu^{n+1},\quad \mu^{n+1}=\mu_{+}(Q^{n+1}) - \mu_{-}(Q^{n}).\tag{33}\]

It is clear that each term in \(F_{-}\) is convex given that the coefficient of each term is positive. For \(F_{+}\), the convexity is established under certain conditions on the coefficients.

Under appropriate assumptions on \({M}^{\mathrm{qent},n}\), the resulting scheme admits a unique solution at each time step.

To compute the quasi-entropy closure approximation, it is necessary to provide an initial guess of \((a_1,a_2,a_3)\) for given \((s,b)\) such that they lie within the domain of \(\Xi_4\), which is provided below.

5.2 Numerical results↩︎

We examine the case where the system is homogeneous in the \(z\)-direction, so that \(Q(\mathbf{x})\) depends only on \((x,y)\). The computational domain \(\Omega = [0,2\pi]\times[0,2\pi]\) is equipped with periodic boundary conditions. For \(\gamma_1\), \(\gamma_2\), \(\gamma_3\) in the scheme 32 –33 , they are chosen so that inequalities hold in Theorem 14. The discretization in space is carried out using the Fourier expansion with \(N=32\) Fourier modes. The nonlinear equations resulting from the scheme are solved by Newton’s iteration with the tolerance \(10^{-6}\).

For the influences of the fourth-order tensor \({M}\), we will compare the results with the \(L^2\) gradient flow without incorporating the fourth-order tensor \({M}\), written as \[\label{eq:L2gf} \frac{\partial Q}{\partial t}=-\mathcal{P}\mu_Q.\tag{35}\] Here, \(\mu_Q\) is still the variational derivative of \(F\) given in 28 –29 , and \(\mathcal{P}\) is the projection onto the symmetric traceless tensors, i.e. \((\mathcal{P}U)_{ij}=(U_{ij}+U_{ji})/2-U_{kk}\delta_{ij}/3\). The first-order scheme proposed in [39] is adopted to solve 35 . Meanwhile, we are also interested in the influence of the cubic elastic \(c_{24}\) term. Accordingly, we define labels for different cases: \({M}^{\mathrm{qent}}\) or \(L^2\) indicates whether the equation 27 or 35 we are solving; Cub or Quad indicates \(c_{24}\ne 0\) or \(c_{24}=0\) (if \(c_{24}\) is specified as a nonzero value, using the label Quad means that \(c_{24}\) is set to zero). For example, \({M}^{\mathrm{qent}}\)–Cub means the results for 27 with some nonzero \(c_{24}\).

Accuracy. We first carry out an accuracy test as a validation. Choose the coefficients \(c_{02}=15\), \(c_{21}=0.16\), \(c_{22} = 0.02\), \(c_{24} = 0.015\), and the initial condition, \[Q(x,y)=\big(s_0+0.05\sin(x)\sin(y)\big)\left(\mathbf{e}_1^2-\mathbf{i}/3\right),\] where \(s_0\) is the minimizer of the bulk energy 16 . The solution at \(t=1\) is computed using various time steps. The reference solution is obtained by \(\delta t=2^{-10}\). The error plotted in Figure 1 indicates first-order accuracy.

Isotropic–uniaxial nematic interfaces. We next examine the evolution of the interface between the isotropic and the uniaxial nematic phases. To this end, we need to carefully choose the bulk energy coefficient \(c_{02}\) such that two phases coexist. The value is \(c_{02}\approx 13.1117\), and the corresponding uniaxial nematic solution gives \(s_0\approx 0.1836\). The elastic coefficients are chosen as \(c_{21}=1.6\times 10^{-3}\), \(c_{22}=2.0\times 10^{-3}\), \(c_{24}=1.5\times 10^{-3}\). The initial condition represents a sharp square-shaped isotropic–nematic interface (with a rotation of \(\pi/4\)), given by \[\label{eq:interface-Q} Q(\mathbf{x},0)=s(\mathbf{x}) \left(\mathbf{e}_1^2-\frac{\mathbf{i}}{3}\right), s(\mathbf{x})=\left \{ \begin{array}{ll} s_0, & |x-\pi|+|y-\pi|\le 1.5, \\ 0, & \text{otherwise}. \end{array}\right.\tag{36}\] The preferred direction of the uniaxial nematic phase is the \(x\)-direction \(\mathbf{e}_1\).

We compute the four cases (with and without \({M}\), zero and nonzero \(c_{24}\)) with \(\delta t=5\times 10^{-2}\). We observe that in all four cases, the interface gradually shrinks and eventually arrive at a globally isotropic state. The stages of evolution are similar for the four cases, which we illustrate in Figure 2 using a few snapshots of \({M}^{\mathrm{qent}}\)–Cub. The colorbar in Figure 2 represents the largest eigenvalue of \(Q\). After the sharp corners of the region become rounded, the square-shaped region gradually transforms into an ellipse. The long axis of the ellipse coincides with the direction of the uniaxial nematic phase, also drawn in Figure 2. The ellipse region then shrinks and finally vanishes.

On the other hand, the rate of interface evolution is different for the four cases. More specifically, it turns out that the presence of \({M}^{\mathrm{qent}}\) greatly accelerates the evolution, while the difference generated by the cubic \(c_{24}\) term is not that evident. This can be acquired by plotting the area occupied by the nematic phase, defined as the \(0.1\)-contour of the maximum eigenvalue of \(Q\) (Figure 3 (a)). The evolution of the free energy (Figure 3 (b)) further reflects evolution rate after the maximum eigenvalue of \(Q\) grows less than \(0.1\). We can see that the energy dissipation is met for the four cases. After the energy of the two \({M}^{\mathrm{qent}}\) cases has reached that of the uniform isotropic state, for the two \(L^2\) cases it still takes some time. Meanwhile, the rate of evolution is slightly faster for the two Cub cases. To further comprehend the effects caused by \({M}^{\mathrm{qent}}\), we plot the eigenvalues of \(\Psi_4({M}^{\mathrm{qent}})\) in Figure 3 (c). The acceleration shall be closely related to the fact that the maximum eigenvalue is roughly five times of the minimum eigenvalue.

Defect evolutions. Set \[c_{02} = 20,\,c_{21} = 1.6,\,c_{22} = 2.0,\,c_{24} = 1.5,\,s_0=0.63.\] Let \[\mathbf{u}(x,y)=\Big(\cos\big(2(x-\frac{3}{16}\pi)\big)-1,\cos(y-\pi)-1,0\Big)^t,\quad \mathbf{n}(x,y)=\frac{\mathbf{u}}{|\mathbf{u}|}.\] It generates discontinuities of \(\mathbf{n}(\mathbf{x})\) at two points \(({3\pi}/{16},\pi)\) and \(({19\pi}/{16},\pi)\). We choose an initial value \(Q_0\) such that it possesses the principal eigenvector \(\mathbf{n}(x,y)\) and takes zero at two points where \(\mathbf{n}\) is discontinuous, \[Q_0(x,y)=s_0 \big(1-\exp(-10|\mathbf{u}|)\big)(\mathbf{n}^2-\frac{\mathbf{i}}{3}).\] The computation is carried out with \(\delta t=10^{-5}\).

We first examine the principal eigenvector of \(Q\). In the four cases (with and without \({M}\), zero and nonzero \(c_{24}\)), we find that the evolution of the principal eigenvector largely goes through the same stages. As an example, the \({M}^{\mathrm{qent}}\)–Cub case is drawn in Figure 4, illustrating the disengagement of two discontinuities.

Despite the similarities for the four cases given by the evolution stages of the principal eigenvector, there are significant differences between the configurations. To reveal the differences, we investigate the biaxiality quantified by \(1-6{{\left(\operatorname{tr}(Q^3)\right)}^2}/{{\left(\operatorname{tr}(Q^2)\right)}^3}\in [0,1]\) for nonzero \(Q\). According to the chosen initial value, the biaxiality is zero at \(t=0\) and emerges with \(t\) increasing. At \(t=0.02\), the biaxiality for the four cases is presented in Figure 5, where it is easy to notice that the two \(L^2\) cases have greater biaxiality than the two \({M}^{\mathrm{qent}}\) cases. This difference actually does not come from the different rates of defect evolution. Indeed, we plot the area enclosed by the \(0.1\)-contour as a function of time \(t\) (Figure 6 (a)), and find that the area is significantly smaller for the two \({M}^{\mathrm{qent}}\) cases. In other words, one major effect of \({M}^{\mathrm{qent}}\) is that it leads to smaller biaxial regions during the defect evolution.

The free energy evolution is plotted in Figure 6 (b), which still decreases with time. Combined with Figure 6 (a), it is clear that the disengagements of defects are much faster for the two \({M}^{\mathrm{qent}}\) cases. In contrast, whether to include a cubic \(c_{24}\) term in the free energy appears not to make a big difference. In addition, since the numerical scheme is nonlinear, we also look into the number of nonlinear iterations for Newton’s method (Figure 6 (c)). For most time steps, only two iterations are needed, thus nonlinearity does not seem to affect the efficiency.

6 Concluding Remarks↩︎

For the tensor dynamics derived from the molecular models, we propose the quasi-entropy closure approximation, which is an elementary function proposed for the entropy term in the free energy. It is done by reformulating the Bingham closure as a minimization problem of the original entropy function, followed by substituting the original entropy with the quasi-entropy. The quasi-entropy closure approximation has the same symmetry properties as the Bingham closure. Together with the quasi-entropy in the free energy, in the low velocity approximation we construct a tensor gradient flow with the quasi-entropy. Such a tensor model maintains the gradient flow structure of the molecular model and the tensor model with the Bingham closure. In particular, the \(Q\) is constrained within \(\mathscr{Q}_{\mathrm{phys}}\), the dissipation operator given by the higher-order tensor is guaranteed to be positive definite, and the elastic energy may possess a cubic term that is bounded from below as a result of \(Q\in\mathscr{Q}_{\mathrm{phys}}\).

The quasi-entropy closure approximation can be done by minimizing an elementary function w.r.t. three variables \(a_1\), \(a_2\), \(a_3\). As a result, we write down a first-order-in-time scheme preserving the eigenvalue constraints and energy dissipation, which can be implemented easily from the fact that the discretization of the closure approximation is explicit and thus decoupled from the scheme. We examine the evolution of interface and defects, finding that the fourth-order tensor would significantly affect the dynamical behaviors.

When the velocity is not discarded, the tensor dynamics forms a coupled system with the Navier–Stokes equations, where higher-order tensors play different roles. The quasi-entropy closure approximation can also be incorporated in this system, for which it requires to study whether the essential structures are maintained, and how to construct efficient and accurate numerical methods. After that, it would be available to systematically carry out numerical simulations and make meaningful comparisons with previous models. The quasi-entropy closure approximation can also be extended to other rigid molecules without axisymmetry. In this case, more order parameter tensors are included [22], [23], [40], [49], [51], leading to many more higher-order tensors to be handled by the closure approximation. We expect to investigate these problems in the near future.

Acknowledgements↩︎

This work is partially supported by Beijing Natural Science Foundation (No. JQ25002), National Natural Science Foundation of China (Nos. 12288201, 12371414, 12571412, 12171041), National Key R&D Program of China (No. 2023YFA1008802), the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDB0510201).

References↩︎

---

1. Laboratory of Mathematics and Complex Systems and School of Mathematical Sciences, Beijing Normal University, Beijing 100875, China. Email:yongyong.cai@bnu.edu.cn↩︎

2. SKLMS & NCMIS, Institute of Computational Mathematics and Scientific/Engineering Computing (ICMSEC), Academy of Mathematics and Systems Science (AMSS), Chinese Academy of Sciences, Beijing, China. Email:xujie@lsec.cc.ac.cn↩︎

3. Laboratory of Mathematics and Complex Systems and School of Mathematical Sciences, Beijing Normal University, Beijing 100875, China. Email: haixinzhang@mail.bnu.edu.cn↩︎