Landau-de Gennes Modelling of Confinement Effects and Cybotactic Clusters in Bent-Core Nematic Liquid Crystals

TL;DR

This paper develops a generalized Landau-de Gennes model for bent-core nematic liquid crystals, capturing phase transitions, confinement effects, and cybotactic clusters, with numerical insights into director profiles and cluster properties.

Contribution

It introduces a new two-order-parameter model for BCNs, extending previous models to include cluster effects and external field influences.

Findings

Identifies two phase transitions in 3D BCN systems.

Demonstrates tunable director profiles via the model.

Provides numerical analysis of confined BCN systems.

Abstract

We study bent-core nematic (BCN) systems in two-dimensional (2D) and three-dimensional (3D) settings, focusing on the role of cybotactic clusters, phase transitions, confinement effects and applied external fields. We propose a generalised version of Madhusudana's two-state model for BCNs in [Madhusudhana NV, Physical Review E, 96(2), 022710] with two order parameters: $\mathbf{Q}_g$ to describe the ambient ground-state (GS) molecules and $\mathbf{Q}_c$ to describe the additional ordering induced by the cybotactic clusters. The equilibria are modelled by minimisers of an appropriately defined free energy, with an empirical coupling term between $\mathbf{Q}_g$ and $\mathbf{Q}_c$. We demonstrate two phase transitions in spatially homogeneous 3D BCN systems at fixed temperatures: a first-order nematic-paranematic transition followed by a paranematic-isotropic phase transition driven by the GS-cluster coupling. We also numerically compute and give heuristic insights into solution landscapes of confined BCN systems on 2D square domains, tailored by the GS-cluster coupling, temperature and external fields. This benchmark example illustrates the potential of this generalised model to capture tunable director profiles, cluster properties and potential biaxiality induced by antagonistic $\mathbf{Q}_g$ and $\mathbf{Q}_c$-profiles.

Figures

Equations

$${\mathbf{Q}}_g=\sqrt{\frac{3}{2}}{S}_g\left({\mathbf{n}}_g\otimes {\mathbf{n}}_g-\frac{I_3}{3}\right)$$

$${\mathbf{Q}}_g=\sqrt{2}{S}_g\left({\mathbf{n}}_g\otimes {\mathbf{n}}_g-\frac{I_2}{2}\right),$$

$${\mathbf{Q}}_c=\sqrt{2}{S}_c\left({\mathbf{n}}_c\otimes {\mathbf{n}}_c-\frac{I_2}{2}\right),$$

$$F({\mathbf{Q}}_g,{\mathbf{Q}}_c)={\int }_{\Omega }{K}_g|\nabla {\mathbf{Q}}_g{|}^2+K_c|\nabla {\mathbf{Q}}_c{|}^2+(1-a_x)$$

$$\times \{\frac{a_g}{2}(T-T^{\ast })\text{tr}{\mathbf{Q}}_g^2-\frac{B_g}{3}\text{tr}{\mathbf{Q}}_g^3+\frac{C_g}{4}(\text{tr}{\mathbf{Q}}_g^2{)}^2$$

$$-\sqrt{2}{E}_{el}{\mathbf{e}}^T{\mathbf{Q}}_g\mathbf{e}\}$$

$$+\frac{a_x}{N_c}\left\{-(1-a_x)\gamma \text{tr}({\mathbf{Q}}_g{\mathbf{Q}}_c)+\frac{{\alpha }_c}{2}\text{tr}{\mathbf{Q}}_c^2+\frac{{\beta }_c}{4}(\text{tr}{\mathbf{Q}}_c^2{)}^2\right\}$$

$$-\sqrt{2}{a}_{x}J{E}_{el}{\mathbf{e}}^T{\mathbf{Q}}_c\mathbf{e}d\mathbf{x},$$

$$a_x=\exp (-E_{ex}/k_{B}T)/\left[1+\exp \left(\frac{-E_{ex}}{k_{B}T}\right)\right],$$

$$\tilde{F}({\mathbf{Q}}_g,{\mathbf{Q}}_c)=\frac{{\lambda }^{2-d}F}{K_g}={\int }_{\tilde{\Omega }}|\tilde{\nabla }{\mathbf{Q}}_g{|}^2+|\tilde{\nabla }{\mathbf{Q}}_c{|}^2+2{\tilde{\lambda }}^2{f}_{b}d\tilde{\mathbf{x}}.$$

$$f_b$$

$$+\frac{M}{2}\text{tr}{\mathbf{Q}}_c^2+\frac{N}{4}(\text{tr}{\mathbf{Q}}_c^2{)}^2-D\text{tr}({\mathbf{Q}}_g{\mathbf{Q}}_c)-\sqrt{2}PE{\mathbf{e}}^T{\mathbf{Q}}_c\mathbf{e}$$

$$A=(1-a_x)a_g(T-T^{\ast })/C_g,B=(1-a_x)B_g/C_g,$$

$$C=(1-a_x),D=\frac{a_x}{N_c}(1-a_x)\gamma /C_g,E=(1-a_x)E_{el}/C_g,$$

$$M=\frac{a_x}{N_c}{\alpha }_c/C_g,N=\frac{a_x}{N_c}{\beta }_c/C_g,P=a_{x}J/(1-a_x),$$

$${\tilde{\lambda }}^2=C_g{\lambda }^2/(2K_g),$$

$$f_b$$

$$+\frac{M}{2}{S}_c^2+\frac{N}{4}{S}_c^4-D{S}_g{S}_c-PE{S}_c,$$

$$A{S}_g-B{S}_g^2+C{S}_g^3-D{S}_c-E=0,$$

$$M{S}_c+N{S}_c^3-D{S}_g-PE=0.$$

$$\left(\begin{array}{cc}A-2B{S}_g+3C{S}_g^2& -D\\ -D& M+3N{S}_c^2\end{array}\right)$$

$$\begin{gathered} f_b=\frac{A}{2}\text{tr}{\mathbf{Q}}_g^2-\frac{B}{3}\text{tr}{\mathbf{Q}}_g^3+\frac{C}{4}(\text{tr}{\mathbf{Q}}_g^2{)}^2 \\ +\frac{M}{2}\text{tr}{\mathbf{Q}}_c^2+\frac{N}{4}(\text{tr}{\mathbf{Q}}_c^2{)}^2-D\text{tr}({\mathbf{Q}}_g{\mathbf{Q}}_c) \end{gathered}$$

$$A{\mathbf{Q}}_g-B\left({\mathbf{Q}}_g^2-\frac{1}{3}\text{tr}{\mathbf{Q}}_g^2\mathbf{I}\right)+C\left(\text{tr}{\mathbf{Q}}_g^2\right){\mathbf{Q}}_g=D{\mathbf{Q}}_c,$$

$$\left(M+N\text{tr}{\mathbf{Q}}_c^2\right){\mathbf{Q}}_c=D{\mathbf{Q}}_g.$$

$$A\sigma {\mathbf{Q}}_c-B{\sigma }^2\left({\mathbf{Q}}_c^2-\frac{1}{3}\text{tr}{\mathbf{Q}}_c^2\mathbf{I}\right)+C{\sigma }^3\left(\text{tr}{\mathbf{Q}}_c^2\right){\mathbf{Q}}_c=D{\mathbf{Q}}_c.$$

$$\Omega \times [0,h]$$

$$\mathbf{Q}$$

$$+s_3(2\mathbf{z}\otimes \mathbf{z}-\mathbf{x}\otimes \mathbf{x}-\mathbf{y}\otimes \mathbf{y}).$$

$$\mathbf{Q}=\left[\begin{array}{cc}{s}_1& s_2\\ s_2& -s_1\end{array}\right],$$

$${\mathbf{Q}}_g=\left[\begin{array}{cc}{q}_1& q_2\\ q_2& -q_1\end{array}\right],{\mathbf{Q}}_c=\left[\begin{array}{cc}{p}_1& p_2\\ p_2& -p_1\end{array}\right],\mathbf{e}=(e_1,e_2{)}^T,$$

$$q_1=S_g\cos 2\theta /\sqrt{2},q_2=S_g\sin 2\theta /\sqrt{2},$$

$$p_1=S_c\cos 2\varphi /\sqrt{2},p_2=S_c\sin 2\varphi /\sqrt{2}.$$

$$\tilde{F}({\mathbf{Q}}_g,{\mathbf{Q}}_c)={\int }_{\tilde{\Omega }}2(|\tilde{\nabla }{q}_1{|}^2+|\tilde{\nabla }{q}_2{|}^2+|\tilde{\nabla }{p}_1{|}^2+|\tilde{\nabla }{p}_2{|}^2)$$

$$+2{\tilde{\lambda }}^2(A(q_1^2+q_2^2)+C(q_1^2+q_2^2{)}^2$$

$$-\sqrt{2}E(q_1{e}_1^2+2q_2{e}_1{e}_2-q_1{e}_2^2)$$

$$+M(p_1^2+p_2^2)+N(p_1^2+p_2^2{)}^2-2D(p_1{q}_1+p_2{q}_2)$$

$$-\sqrt{2}PE(p_1{e}_1^2+2p_2{e}_1{e}_2-p_1{e}_2^2))d\tilde{x},$$

$$\tilde{\Delta }{q}_1$$

$$\tilde{\Delta }{q}_2$$

Full text

††preprint: APS/123-QED

I Introduction

Liquid crystals (LCs) are fascinating mesophases with physical properties intermediate between the conventional solid and liquid phases [2]. LCs are partially ordered materials, i.e., they are more ordered than liquids and less ordered than crystalline solids. The partial ordering can manifest in special material directions, layered structures, chirality, columnar structures etc. Nematic liquid crystals (NLCs) are the simplest and perhaps, the most common type of LCs composed of rod-like molecules that tend to preferentially align along some distinguished directions, referred to as directors, rendering long-range orientational order [2]. Smectic LCs are effectively layered LCs, i.e., the LC molecules organise into layers and there is a preferred director or nematic ordering within each layer. Then, we have cholesterics which are twisted nematics such that the director naturally twists in space and imparts macroscopic chirality. Molecular shape can have strong consequences on macroscopic LC properties and recently, there has been tremendous interest in bent-core nematics (BCNs) [11]. In non-technical terms, a BCN molecule is composed of two connected rods with an opening angle. This kinked shape opens a plethora of possibilities - new phases, biaxiality, natural chirality, non-Newtonian rheology and exotic morphologies [11, 3, 4].

BCNs were discovered in the 1990’s; they are remarkable materials that can exhibit different types of order — polar, chiral, octupolar, biaxial order etc. They can exhibit structural transitions to more ordered phases —such as twist‑bend and smectic—depending on temperature and molecular architecture  [11]. These BCN phases offer tremendous potential for electro-optic devices, rheological applications and novel functional materials. Cybotactic clusters are an intrinsic feature of BCNs; they are smectic-like nanoclusters such that BCN molecules get locked into clusters because of their kinked shape. The existence of cybotactic clusters has been confirmed by multiple experimental methods, e.g., SAXS, cryo-TEM methods, electro-optic measurements [26, 24, 5] etc. Cybotactic clusters are small, usually tens of nanometers in lateral dimensions and persist in the nematic and sometimes, even in the isotropic BCN phase. Cybotactic clusters render additional ordering to the conventional nematic phase and this ordering clearly manifests in optical images (four-lobed structures, characteristic of local smectic‑C‑like ordering, in SAXS measurements), rheological properties (non-Newtonian properties), elastic properties (suppressed twist elastic constants leading to spontaneous chirality), macroscopic biaxiality [11] etc. As such, it is essential to quantify the additional ordering induced by the cybotactic clusters and how this ordering modulates the conventional ambient nematic ordering.

In [17], Madhusudana proposes a two-state model for BCNs, based on the hypothesis that there are two types of molecules: the ground-state (GS) molecules that constitute the out-of-cluster ambient environment and the excited molecules (ES) that constitute the cybotactic clusters. There are two scalar order parameters, $S_{g}$ and $S_{c}$, to describe the orientational ordering of the GS molecules and the ES molecules respectively. Madhusudana assumes spatially uniform $S_{g}$ and $S_{c}$ profiles, i.e., he works with spatially homogeneous BCN systems along with various other assumptions — fixed time-independent and temperature-independent cluster sizes, an empirical coupling term between $S_{g}$ and $S_{c}$ and no explicit information about the BCN directors. Madhusudana proposes a phenomenological free energy comprising two Ginzburg-Landau type potentials for $S_{g}$ and $S_{c}$ with various temperature-dependent and material-dependent parameters and an empirical GS-ES coupling term that is parameterised by a coupling coefficient. The physically observable spatially homogeneous BCN configurations are modelled by minimisers of this two-state free energy. The most important finding concerns the existence of a new paranematic phase with small values of $S_{g}$ and $S_{c}$, just before the transition to the disordered isotropic phase. This model also effectively captures other experimental trends reported in [28] — enhanced order parameters induced by the cybotactic clusters, a nematic-paranematic phase transition followed by a paranematic-isotropic phase transition with increasing temperature and non-linear dependence of order parameters on external magnetic fields.

However, Madhusudana’s model cannot account for confinement effects or spatial inhomogeneities. In [21], the authors propose a one-dimensional model for spatially inhomogeneous BCN samples, including elastic energy terms for $S_{g}$ and $S_{c}$. They study BCN samples in channel geometries with fixed boundary conditions for $S_{g}$ and $S_{c}$; their numerical results demonstrate that confinement or geometric frustration promotes interior clustering or enhances interior ordering and that the GS-ES coupling can promote ordered director profiles even at high temperatures contrary to the expectation of isotropy or disorder at high temperatures. In [19], the authors study the effects of doping on cybotactic clusters along with the emergence of polar order in such BCN systems. They perform extensive experiments on phase transitions, birefringence, order parameter measurements and dielectric spectroscopy data and propose a simple Madhusudana-type free energy to partially describe their experimental results. In this paper, we build on these works to propose a generalised free energy for confined BCN systems in two-dimensional (2D) and three-dimensional (3D) settings. By analogy with Landau-de Gennes theory, we define two tensor order parameters: $\mathbf{Q}_{g}$ to describe the orientational ordering of the GS molecules and $\mathbf{Q}_{c}$ to describe the additional ordering induced by the ES molecules within cybotactic clusters [18]. The key advantage is that $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ not only contain information about $S_{g}$ and $S_{c}$ but also contain information about the GS and ES directors, that are the preferred directions of average alignment of the GS molecules and new distinguished material directions induced by the cybotactic clusters. The free energy has an elastic energy to penalise spatial inhomogeneities and a bulk energy that determines the GS and ES ordering as a function of temperature, GS-ES coupling and other parameters, along with terms that account for external (electric or magnetic) fields. The bulk energy is a direct generalisation of Madhusudana’s free energy and is composed of Ginzburg-Landau type potentials for $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ respectively and a $\mathbf{Q}_{g}$-$\mathbf{Q}_{c}$ coupling term. The Ginzburg-Landau potentials dictate that $S_{g}$ and $S_{c}$ increase with decreasing temperatures and approach zero with increasing temperature (the isotropic or disordered phase). The GS-ES coupling and external field terms can strongly boost ordering, even for high temperatures. The physically observable textures are modelled by minimisers of this generalised energy, subject to the imposed boundary conditions. The minimisers and the critical points of the free energy are solutions of boundary-value problems for the associated Euler-Lagrange equations — typically a system of coupled and nonlinear partial differential equations and the rest of this manuscript is devoted to the study of critical points of this generalised free energy in physically motivated settings.

The rest of the paper is organised as follows. In Section II, we introduce a generalised Landau-de Gennes type free energy for BCN systems. In Section III, we study phase transitions in spatially homogeneous BCN samples, with and without external fields, in 2D and 3D. Notably, we find two phase transitions in spatially homogeneous 3D BCN samples driven by increasing GS-ES coupling strength, at high temperatures. The GS-ES coupling strength is measured in terms of a parameter $\gamma$. If $\gamma=0$, then $S_{g}=S_{c}=0$ for spatially homogeneous samples at high temperatures (defined by $A>0$ in (5)). As $\gamma$ increases, the energy minimisers have positive values of $S_{g}$ and $S_{c}$ even for high temperatures. At a fixed high temperature, increasing $\gamma$ has a strong ordering effect that counteracts the disordering effects of high temperature. This is analogous to the nematic-paranematic phase transition and the paranematic-isotropic phase transition with increasing temperature, at fixed $\gamma$, as reported in [17]. The external fields typically boost the values of $S_{g}$ and $S_{c}$, shift the bifurcation points and we lose the paranematic-isotropic phase transition. Given that the values of $S_{g}$ and $S_{c}$ are very small in the paranematic phase, it is unclear if the paranematic phase would be distinguishable from the conventional isotropic phase in experiments. In Section IV, we study solution landscapes of confined BCN systems on a square domain with Dirichlet boundary conditions for $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$. This is a benchmark well-studied example for conventional nematics; see the papers [30, 16, 14]. A solution landscape is a network of critical points of the free‑energy surface, revealing possible pathways between the different stable solutions and unstable saddle points of the free energy. The unstable saddle points play a crucial role in the system dynamics. We study the interplay between domain size, temperature, the GS-ES coupling and external field strength on the energy minimisers and the associated solution landscapes. We recover some familiar $\mathbf{Q}_{g}$ profiles from the studies on conventional nematics, some new $\mathbf{Q}_{g}$ profiles tailored by the additional physics in this model, co-aligned $(\mathbf{Q}_{g},\mathbf{Q}_{c})$ profiles with co-aligned GS and ES directors and antagonistic $(\mathbf{Q}_{g},\mathbf{Q}_{c})$-profiles with mutually perpendicular GS and ES directors. Our simple examples illustrate the tremendous potential of tuning model parameters to get exotic morphologies, some of which could lead to macroscopic biaxiality e.g., the antagonistic $(\mathbf{Q}_{g},\mathbf{Q}_{c})$ profiles. We conclude with a discussion of the limitations of the model and future improvements in Section VI.

II Landau-de Gennes Modelling for Bent-Core Nematic Systems

In this paper, we propose a two-state generalised Landau-de Gennes (LdG) type free energy for confined BCN systems building on Madhusudana’s model in [17]. We define a GS LdG order parameter, $\mathbf{Q}_{g}$, which is a symmetric, traceless $d\times d$ matrix, where $d$ is the spatial dimension. The eigenvectors of $\mathbf{Q}_{g}$ model the GS nematic directors, interpreted as the locally preferred directions of averaged GS molecular alignment in space and the associated eigenvalues measure the degree of orientational ordering about the GS directors. We define the GS director to the eigenvector of $\mathbf{Q}_{g}$ with the largest positive eigenvalue. The GS is said to be in the isotropic phase if $\mathbf{Q}_{g}=0$. It is in the uniaxial phase if $\mathbf{Q}_{g}$ has two equal non-zero eigenvalues for $d=3$ and $\mathbf{Q}_{g}$ can be written as

[TABLE]

where $\mathbf{n}_{g}$ is the eigenvector with the non-degenerate eigenvalue and $S_{g}$ is proportional to the eigenvalue associated with $\mathbf{n}_{g}$. If $S_{g}>0$, then $\mathbf{n}_{g}$ is the GS director. $\mathbf{Q}_{g}$ can be biaxial for $d=3$ (and not for $d=2$) when it has three distinct eigenvalues. In 2D,

[TABLE]

where $\mathbf{n}_{g}$ is the GS director and $S_{g}$ is positive. When $d=2$, the defect set is identified with the nodal set of $S_{g}$. For completeness, $I_{d}$ is the identity matrix in $d$ dimensions.

We treat the cybotactic clusters as uniformly distributed sub-nanometric/very small clusters, such that the size of the cluster is much smaller than the spacing between the clusters and the volume fraction of the clusters is much smaller than unity. This is the so-called “dilute limit” in homogenization theory [1], consistent with the assumptions in Madhusudana’s paper. In this case, the highly localised cybotactic clusters interact with the surrounding GS molecules by means of an averaged field, the $\mathbf{Q}_{c}$ - tensor field, referred to as the cluster order parameter in this paper. In principle, $\mathbf{Q}_{c}$ contains information about the cluster shapes and the GS-ES interactions, assuming that the clusters do not interact with each other. For example, in $d=2$, we have

[TABLE]

where $\mathbf{n}_{c}$ is referred to as the ES/cluster director and $S_{c}$ is positive. As an immediate consequence, we have two tensor order parameters, $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$, as opposed to two scalar order parameters, $S_{g}$ and $S_{c}$ in [17], that can describe complex information about the GS and cluster directors, structural transitions and confinement effects in inhomogeneous systems which are outside the scope of the two-state model in [17].

The corresponding generalised LdG-type free energy is

[TABLE]

where $\Omega$ is a $d$-dimensional bounded domain with Lipschitz boundary and $d=2,3$. The first two terms are the elastic energy density terms associated with $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ respectively, i.e., we assume a one-constant elastic energy density to penalise spatial inhomogeneities in $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ such as strong GS-ES interactions, interior defects, geometric frustration etc. $K_{g}$ and $K_{c}$ denote the positive elastic constants associated with $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ respectively. We note that other choices of anisotropic elastic energy densities are also possible, see [7]. For BCN materials, the bend elastic constant $K_{33}$ is often considerably smaller than the splay ($K_{11}$) or twist ($K_{22}$) constants, reflecting an intrinsic bend‑softening associated with the molecular geometry. The mole fraction of the ES molecules $a_{x}$ is given by

[TABLE]

as $E_{ex}$ is the excitation energy of the ES molecules, $k_{B}$ is Boltzmann’s constant and $T$ is the temperature. $a_{x}=0.1$ throughout the paper consistent with [21]. Further, there maybe an intrinsic $a_{x}$-dependence of the elastic constants. In the interest of brevity and accounting for uncertainties in the properties of the cybotactic clusters, including their elastic constants, we use the one-constant elastic energy density as a first approximation. The parameters $a_{g}$, $B_{g}$, $C_{g}$, and $T^{*}$ are the usual LdG parameters that describe the first-order nematic-isotropic transition for the GS molecules and $\gamma$ is the coupling parameter between the GS molecules and the clusters [21, 17]. We note that the term $\textrm{tr}\left(\mathbf{Q}_{g}\mathbf{Q}_{c}\right)$ is a generic coupling term for systems with two order parameters that is consistent with frame-indifference and material symmetry requirements [1]. $\alpha_{c}$ and $\beta_{c}$ are coefficients for the saturation terms to ensure that $|\mathbf{Q}_{c}|^{2}$ is reasonably bounded, for at least a certain range of $\gamma$. There is no evidence of first-order isotropic-nematic phase transitions within cybotactic clusters to date, so the term $tr\mathbf{Q}_{c}^{3}$ is absent. We also note that Madhusudana’s original work does not include a cubic term in $S_{c}$, analogous to a $tr\mathbf{Q}^{3}_{c}$ term. In the subsequent numerical simulations, we adjust the values of $\alpha_{c}$ and $\beta_{c}$ so that the resulting values of $S_{c}$ are reasonable and not too large. $N_{c}$ is the number of ES molecules within each cluster. We also consider the effects of an external field, $\mathbf{E}=E_{a}\mathbf{e}$, whose direction is modelled by a unit-vector $\mathbf{e}$, and the magnitude is $E_{a}$. $E_{el}$ is the external field energy density $(\frac{1}{2}\epsilon_{0}\Delta\epsilon E_{a}^{2})$ where $\epsilon_{0}$ is the free-space permittivity, $\Delta\epsilon$ is a material anisotropy parameter (dielectric or magnetic) for GS molecules and $J$ accounts for the larger shape anisotropy of ES molecules, i.e., $J>1$. In particular, the GS molecules tend to align with $\mathbf{E}$ if $\Delta\epsilon>0$ and orthogonal to $\mathbf{E}$ if $\Delta\epsilon<0$ [2].

We nondimensionalise the free energy above with $\tilde{x}=x/\lambda$ and assume $K_{g}=K_{c}$, where $\lambda$ is a characteristic length scale of the $d$-dimensional domain $\Omega$,

[TABLE]

where

[TABLE]

where

[TABLE]

and $\tilde{\Omega}$ is the rescaled domain. For comparison with previous work on one-dimensional BCN systems in [21], we set $A=\pm 0.04$, $B=0.34$, $C=0.9$, $D=2\times 10^{-3}$, $E=1.6\times 10^{-4}$, $M=9.7\times 10^{-5}$, $N=1.78\times 10^{-3}$, $P=0.133$ which are derived from $K_{g}=K_{c}=K=15\mathrm{pN}=15\times 10^{-12}\mathrm{N}=15\times 10^{-7}\mathrm{dyn}\text{(under one constant approximation)}$; $a_{g}=0.04\times 10^{7}/4\ \mathrm{ergs/(cm^{3}K)}$, $B_{g}=1.7\times 10^{7}/4\ \mathrm{ergs/cm^{3}}$, $C_{g}=4.5\times 10^{7}/4\ \mathrm{ergs/cm^{3}}$, $\alpha_{c}=0.22\times 10^{7}/4\ \mathrm{ergs/cm^{3}}$, $\beta_{c}=4.0\times 10^{7}/4\ \mathrm{ergs/cm^{3}}$, and $\gamma=5.0\times 10^{7}/4\ \mathrm{ergs/cm^{3}}$; $T^{*}=355K$, $N_{c}=50$, $J=1.2$, $E_{el}=2000\mathrm{ergs/cm^{3}}$, $E_{ex}=1.1\times 10^{-13}\mathrm{ergs}$, and $T=360\mathrm{K}$ (for the high temperature case with $A>0$ ) and $350\mathrm{K}$ (for the low temperature case with $A<0$). In subsequent sections, we focus on phase transitions in spatially homogeneous systems and confinement effects in square domains as a function of temperature $A$, the GS-ES coupling strength measured by $D$ and external field strength, $E$.

For clarity, we temporarily refer to the $F$ from Eq. (4), where $\Omega$ is a $d$-dimensional region, as $F_{d}$. In three dimensions, the free energy $F_{3}$ has the unit of energy (erg), whereas in two dimensions, $F_{2}$ has the unit of force (N). This effectively assumes a uniform extension along the third dimension, so the true three‑dimensional energy can be written as $F_{3}=F_{2}\times(\text{thickness})$. After nondimensionalization in Eq. (5), however, the factor $\lambda^{2-d}$, which depends on the spatial dimension $d$, renders the dimensionless free energy $\tilde{F}$ unit‑free in both two and three dimensions.

III Phase Transitions in the Generalised LdG model in 3D and 2D settings

In [17], the author studies phase transitions in spatially homogeneous BCN systems as a function of temperature, based on the two-state free energy in terms of $S_{g}$ and $S_{c}$. The generalised LdG energy in (4) reduces to Madhusudana’s two-state model if we assume that $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ are spatially constant uniaxial tensors and the directors $\mathbf{n}_{g}=\mathbf{n}_{c}=\mathbf{e}\equiv constant$, i.e., the bulk energy density in (5) reduces to

[TABLE]

where $B=0$ in 2D settings. Equation (7) is the two-state free energy studied in [17]. The critical points of this bulk energy density are the solutions, $(S_{g},S_{c})$, of the following system of bivariate polynomial equations:

[TABLE]

We use an open-source computational knowledge engine, WolframAlpha [10], to obtain the numerical solutions for the nonlinear system of equations (8)-(9). Since (9) yields $S_{g}$ as a cubic function of $S_{c}$, we substitute it into (8) to obtain a ninth-degree polynomial in $S_{c}$. This results in nine solutions for $S_{c}$. Here we only consider the real ones. The stability of the critical points is determined by the eigenvalues of the Hessian matrix

[TABLE]

If both eigenvalues are positive, the critical point is stable; if there is at least one negative eigenvalue, the critical point is unstable.

In the 3D setting, with a zero external field ($E=0$), the isotropic phase $S_{g}=S_{c}=0$ is always a solution of the coupled algebraic equations (8)-(9), for any reduced temperature $A$ (Fig. 1(a-b)). At the characteristic high temperature $A=0.04$ (Fig. 1(a)), the isotropic solution, $(S_{g},S_{c})=(0,0)$, is the stable solution with the lowest energy for $D$ sufficiently small. When $D\approx 1.88\times 10^{-3}$, the isotropic solution loses stability and undergoes a pitchfork bifurcation into two stable solutions: a weakly ordered solution [17] with small $S_{g}>0$ and $S_{c}>0$ (solid blue line) and a negative solution branch with $S_{g}<0$ and $S_{c}<0$ (solid purple line). We interpret this weakly ordered stable solution as the paranematic-type phase. A bifurcation point is a point at which a small, smooth change in a parameter value causes a sudden qualitative change in the behaviour of the system, such as the appearance or disappearance of critical points and changes in their stability.

As $D$ further increases, we observe a saddle-node bifurcation accompanied by the creation of a stable nematic solution with large $S_{g}>0$ and $S_{c}>0$ (solid yellow line) and one unstable solution (dashed orange line). When $D\approx 3\times 10^{-3}$, the orange unstable branch and blue paranematic stable branch merge and disappear. Then the number of solutions changes from five to three for $D$ large enough. According to the energy plot in Fig. 1(a), there are two phase transitions as $D$ increases. The phase transition occurs at a value of $D$ for which there are two solution branches with equal free energies in (7). When $D\approx 1.88\times 10^{-3}$, there is a second-order phase transition between the isotropic phase and the weakly ordered paranematic-type phase. When $D\approx 2.4\times 10^{-3}$, there is a first-order phase transition between the weakly ordered paranematic-type phase and ordered nematic phase. For a first-order phase transition, the order parameters $S_{g}$ or $S_{c}$ are discontinuous. For a second-order phase transition, $S_{g}$ and $S_{c}$ are continuous, while the derivatives of $S_{g}$ or $S_{c}$ with respect to $D$ are discontinuous. If we informally interpret $D$ as promoting order or as being an analogue of low temperatures, then this is consistent with the picture in [17] — a first-order nematic-paranematic transition with increasing temperature (decreasing $D$) followed by the conventional isotropic-(para)nematic phase transition with increasing temperature (decreasing $D$). We note that $S_{c}$ is significantly larger than $S_{g}$ in magnitude for all solution branches.

At the characteristic low temperature $A=-0.04$ (Fig. 1(b)), we always have three solution branches — the unstable isotropic solution with $S_{g}=S_{c}=0$ (dashed line in orange), positive stable solution branch with positive $S_{g}$ and $S_{c}$ (solid yellow line) and negative stable solution branch with negative $S_{g}$ and $S_{c}$ (solid blue line). The positive solution branch always has the lowest energy in the solution set, corresponding to well-oriented configurations consistently observed in experiments. We note that $S_{g}<0$ (or $S_{c}<0$) describes a state where the GS (ES) molecules are relatively randomly oriented in the plane perpendicular to $\mathbf{n}_{g}$ ($\mathbf{n}_{c}$), where $\mathbf{n}_{g}$ is the eigenvector of $\mathbf{Q}_{g}$ with the negative non-degenerate eigenvalue.

When an external field is applied ($E=1.6\times 10^{-4}$), the isotropic solution branch no longer exists (Fig. 1(c-d)). For $A=0.04$ (Fig. 1(c)), the bifurcation structure and the number of solutions change compared to the zero external field (with $E=0$). When $D$ is small enough, there is a unique weakly-ordered paranematic-type solution with positive $S_{g}$ and $S_{c}$ close to zero. When $D\approx 2.1\times 10^{-3}$, a stable solution branch with relatively large positive values of $S_{g}$ and $S_{c}$, and an unstable solution branch appear through a saddle-node bifurcation. When $D\approx 2.8\times 10^{-3}$, the stable weakly-ordered paranematic solution and the unstable solution merge and disappear. When $D\approx 4\times 10^{-3}$, we note the simultaneous appearance of a stable and unstable negative solution branch, both of which have $S_{g}<0$ and $S_{c}<0$. According to the energy plot in Fig. 1(c), there is only one first-order phase transition between the weakly-ordered paranematic phase (solid blue line) and strongly ordered nematic solution (solid yellow line). Again, decreasing $D$ is analogous to increasing temperature in [17].

For the characteristic low temperature $A=-0.04$ (Fig. 1(d)), we still retain three solution branches and the global bulk energy minimum is attained by the positive solution branch with positive $S_{g}$ and positive $S_{c}$. We note that $S_{g}$ and $S_{c}$ increase with increasing $E$, along the solid yellow solution branch in Fig. 1(b) and (d).

In the 2D setting in Fig. 2, we set $B=0$ in the algebraic equations (8)-(9). Compared to the 3D setting in Fig. 1, we lose the weakly-ordered paranematic phase and the associated saddle-node bifurcations. Without an external field ($E=0$), there are three solution branches — the isotropic solution branch with $S_{g}=S_{c}=0$; the positive nematic solution branch with positive $S_{g}$ and $S_{c}$, and the negative solution branch with $S_{g}$ and $S_{c}<0$. Since $B=0$ and $E=0$, the positive and negative solution branches have equal free energies for $A=\pm 0.04$. For $A=0.04$, the isotropic solution is stable when $D$ is small enough. When $D$ increases further, the isotropic solution loses stability and bifurcates into stable positive and negative ordered solutions. Compared to the 3D setting in Fig. 1(a), we only retain the second-order phase transition in Fig. 2(a). There is no bifurcation for $A=-0.04$ and the isotropic solution branch is unstable for all values of $D$ under consideration (Fig. 2(b)).

With a nonzero external field ($E=1.6\times 10^{-4}$), the positive nematic solution branch is energetically preferred for high and low temperatures. The positive solution branch corresponds to the GS and ES molecules being aligned with the external field. The negative solution branch (with $S_{g},S_{c}<0$) corresponds to the GS and ES molecules being oriented orthogonal to the external field. For $A=0.04$ (Fig. 2(c)), there is no phase transition in the 2D setting, contrary to the 3D setting which exhibits a first-order paranematic-nematic phase transition with increasing $D$. The globally stable solution branch is always the ordered nematic phase (solid yellow line). There are no bifurcations for $A=-0.04$, analogous to the $E=0$ case but the symmetry between the positive and negative solution branches is broken in the energy plot in Fig. 2(d).

We conclude this section by demonstrating that the generalised bulk energy in (4) can only have isotropic or uniaxial critical points for $d=3$, i.e., bulk biaxiality is outside the scope of this relatively simple generalised bulk energy. The critical points of the generalised bulk energy density (with $E=0$)

[TABLE]

are solutions of the coupled systems of partial differential equations:

[TABLE]

Setting $\sigma=(M+N\textrm{tr}\mathbf{Q}_{c}^{2})/D$ and substituting Eq. (12) into Eq. (11) yields

[TABLE]

Let $\mathbf{Q}_{c}$ be a diagonal $3\times 3$ matrix with diagonal entries, $Q_{1},Q_{2},Q_{3}$ respectively and zero non-diagonal entries. The above equation holds for all matrix components of $\mathbf{Q}_{c}$ and repeating the arguments from [18], we deduce that there are at least two equal diagonal entries. This implies that $\mathbf{Q}_{c}$ has at least two equal eigenvalues and hence all critical points of the generalised bulk energy density are either isotropic or uniaxial pairs, $(\mathbf{Q}_{g},\mathbf{Q}_{c})$, as outlined above. Hence, the generalised bulk energy density does not admit biaxial critical points in $d=3$ and bulk biaxiality is outside the scope of this model.

IV Confinement Effects for BCN Spatial Equilibria on Square Domains

Next we consider a thin three-dimensional square well filled with a prototype BCN material:

[TABLE]

where $\Omega$ is a square domain and $h$ is the height of the well. In the $h\to 0$ limit and for certain choices of the surface energies, one can prove that conventional LdG energy minimisers have a fixed eigenvector in $\mathbf{z}$-direction with a fixed constant eigenvalue and all dependent variables are independent of the $\mathbf{z}$-coordinate [6], so that it suffices to consider the reduced LdG tensor— a symmetric, traceless $2\times 2$ matrix with only two degrees of freedom. In other words, the LdG $\mathbf{Q}$-tensor can be written as

[TABLE]

where $\mathbf{x}=(1,0,0)^{T}$, $\mathbf{y}=(0,1,0)^{T}$, $\mathbf{z}=(0,0,1)^{T}$ and $s_{3}$ is a constant in the $h\to 0$ limit. Hence, it suffices to study the reduced LdG tensor as given below:

[TABLE]

Here, we use similar ideas of dimension reduction to study confined BCN systems on a square domain $\Omega=[0,\lambda]^{2}$, where $\lambda$ is the square edge length. We describe the GS and ES-induced ordering by two reduced LdG tensors: $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ as described below:

[TABLE]

where $q_{1},q_{2},p_{1},p_{2}$ only depend on $x$ and $y$, and $\mathbf{e}$ is a constant vector. Recall that in (2) and (3), we have $\mathbf{Q}_{g}=\sqrt{2}S_{g}\left(\mathbf{n}_{g}\otimes\mathbf{n}_{g}-\frac{\mathbf{I}_{2}}{2}\right)$ and $\mathbf{n}_{g}=(\cos\theta,\sin\theta)$ for some director angle $\theta$ in the plane, and analogously $\mathbf{Q}_{c}=\sqrt{2}S_{c}\left(\mathbf{n}_{c}\otimes\mathbf{n}_{c}-\frac{\mathbf{I}_{2}}{2}\right)$ and $\mathbf{n}_{c}=(\cos\phi,\sin\phi)$ for some director angle $\phi$ in the plane. In other words,

[TABLE]

Using the re-scaling, $\tilde{x}=\frac{x}{\lambda}$ and the assumption $K_{g}=K_{c}$, the non-dimensionalised free energy can be written as (refer to (5))

[TABLE]

where $\tilde{\Omega}=[0,1]^{2}$ and $\tilde{\lambda}^{2}=C_{g}\lambda^{2}/(2K_{g})$. The critical points of (15) are solutions of the corresponding Euler-Lagrange equations:

[TABLE]

with Dirichlet boundary conditions

[TABLE]

Here $S_{g}^{*}$ is the positive solution of (8)-(9) with $B=0$, minimising bulk energy density $f_{b}$ in (7). We impose tangential boundary conditions on $\mathbf{Q}_{g}$, i.e., $\mathbf{n}_{g}=(\pm 1,0)$ on $y=0,1$ and $\mathbf{n}_{g}=(0,\pm 1)$ on $x=0,1$. These translate to conflicting boundary conditions for $q_{1}$ on the horizontal and vertical edges. To deal with this, there must be an interior nodal line with $q_{1}=0$, with either $S_{g}=0$ or with $\mathbf{n}_{g}=(1/\sqrt{2},\pm 1/\sqrt{2})$ (captured by $\theta=\pm\frac{\pi}{4}$). For $\mathbf{Q}_{c}$, we impose $S_{c}=0$ on the square edges as in [17, 21, 19], to exclude cybotactic clusters on the square edges.

In Figs. 3 and 4, we focus on the high temperature case, $A=0.04$, and study the effects of coupling parameter $D$ and external field $E$ on the energy minimisers in (15).

For the small domain case in Figure 3, we recover the WORS (Well Order Reconstruction Solution) [12] where $S_{g}=0$ on the two square diagonals, and the BD (Bent Director/Boundary Distortion) solution [27] where $S_{g}=0$ on two curves localised near parallel square edges. We find that the WORS is stabilised by increasing $D$ for positive $A$. As $D$ increases in Fig. 3(a), the average orientational order $\int S_{g}dA$ and $\int S_{c}dA$ increases (Fig. 3(c)), corresponding to the fact that $S_{g}$ and $S_{c}$ increase with increasing $D$ for bulk energy minimisers in Fig 2. In our previous work [22, 9], as the average orientational order increases, the nodal lines of WORS deform from the diagonal lines to shorter lines along opposite square edges (BD). However, we observe the opposite trend here. The boundary condition of $\mathbf{Q}_{g}$ depends on $S_{g}^{*}$ in (17), which is an increasing function of $D$, leading to stronger constraining effects which stabilise symmetric structures like WORS. The symmetry of a profile is measured by $\int q_{1}dA$ and $\int p_{1}dA$ in (Fig. 3(c)). The integral $\int q_{1}dA$ initially increases as $D$ increases due to enhanced ordering and then decreases, indicating the increasing symmetry of the $\mathbf{Q}_{g}$ profile. Since $\tilde{\lambda}^{2}$ is relatively small, $\mathbf{Q}_{c}$ is significantly affected by the isotropic boundary condition in (18), so that $S_{c}$ is around $10^{-4}$ and is much smaller than its value in the homogeneous case in Fig. 2. Hence, $\mathbf{Q}_{c}$ has negligible impact on $\mathbf{Q}_{g}$ and the system is effectively decoupled.

As $E$ increases in (Fig. 3(b)) for a fixed $D$, the effects of the boundary conditions and $S_{g}^{*}$ do not change significantly from left to right. Hence, consistent with our previous work, the interior ordering increases, the nodal lines of $\mathbf{Q}_{g}$ become further apart and the profile is closer to the BD profile. All the four measurements, the average orientational order $\int S_{g}dA$ and $\int S_{c}dA$ and the symmetry measurements $\int q_{1}dA$ and $\int p_{1}dA$, reflect the transition from a WORS to a BD-type profile with increasing $E$.

In Figure 4, we consider a large domain $\tilde{\lambda}^{2}=2500$. For large $\tilde{\lambda}$, we expect energy minimisers of (15) to converge to minimisers of the bulk energy uniformly, almost everywhere away from defects. In other words, we expect $S_{g}$ to converge almost uniformly to $S_{g}^{*}$ in the square interior, for sufficiently large $\tilde{\lambda}$. The bulk energy minimisers for $A=0.04$ are plotted as a function of $D$ in Figure 2, with $E=0$. We note that $D$ is very large in Figure 4 and the corresponding value of $S_{g}^{*}$ is also expected to be large for $D>10^{-2}$. Hence, nodal lines are expensive and the energy-minimising $\mathbf{Q}_{g}$ profile mediates between the conflicting boundary conditions by means of a diagonal profile with $\mathbf{n}_{g}=\left(1/\sqrt{2},\pm 1/\sqrt{2}\right)$ for $D=5\times 10^{-2}$ and $D=10^{-1}$. For $D=10^{-2}$, there is competition between the elastic distortion effects and the interior ordering promoted by large $\tilde{\lambda}$ and $D$. It seems that the energy-minimising $\mathbf{Q}_{g}$ profile prefers a WORS-type profile for $D=10^{-2}$ with diagonal lines of low order. The cybotactic clusters become pronounced in the square interior with large values of $|\mathbf{Q}_{c}|^{2}$ since the effects of the isotropic boundary conditions are dominated by the effects of the bulk energy minimiser in (7) (with $B=0$) in the interior, for large $\tilde{\lambda}$ and large $D$. We note that one cannot achieve an ordered diagonal solution with $D=0$, $\tilde{\lambda}^{2}=2500$ and $A=0.04$, but the diagonal solution can be stabilised by large $D$ at high temperatures.

In Figure 4(b), we fix $D$ and $A=0.04$ and the direction of the external field to be in the diagonal direction, and increase the strength of the external field from left to right. The effect of increasing $E$ is analogous to the effect of increasing $D$ for fixed $A=0.04$; $S_{g}^{*}$ is an increasing function of $E$ for fixed $D$ and the order parameters of the energy minimising $(\mathbf{Q}_{g},\mathbf{Q}_{c})$-profiles converge almost uniformly to the bulk energy minimisers $(S_{g}^{*},S_{c}^{*})$ in (7), for sufficiently large $\tilde{\lambda}$. Hence, we note a transition from a WORS-type profile to a diagonal $\mathbf{Q}_{g}$ solution with increasing $E$, for a fixed $D$ at $A=0.04$. Analogously, the cybotactic clustering is also enhanced in the interior with increasing $E$ or increasing $D$, for $A=0.04$, and this manifests in enhanced values of $|\mathbf{Q}_{c}|^{2}$ in the square interior. We expect the qualitative conclusions to carry over to arbitrary high temperatures captured by fixed $A>0$, while we recover only the isotropic phase at high enough temperatures for fixed $D$ or $E$.

In Figure 4(c), we fix $D$ and $E=0$, and increase the temperature $A$ from left to right. The effect of increasing $A$ is inverse to the effect of increasing $D$ or $E$. According to the results in [17], $S_{g}^{*}$ and $S_{c}^{*}$ are decreasing functions of $A$. Hence, we note a transition from a diagonal $\mathbf{Q}_{g}$ solution to a WORS-type profile with increasing $A$. The cybotactic clustering is weakened in the interior with increasing $A$, and this manifests in reduced values of $|\mathbf{Q}_{c}|^{2}$ in the square interior.

Figures 3 and 4(a-b) focus on the high-temperature regime $A=0.04$. In Fig. 5, we focus on computing solutions of the Euler-Lagrange equations (16) in the low temperature regime $A=-0.04$, where the global bulk energy minimiser is always the positive solution branch, $(S_{g}^{*},S_{c}^{*})$ in Fig 2, with and without an external field. We work with a large domain $\tilde{\lambda}^{2}=2500$ and relatively weak coupling $D=4\times 10^{-3}$ in Figure 5. We note that $\tilde{\lambda}^{2}=2500$ so that we expect energy minimisers to approach $(S_{g}^{*},S_{c}^{*})$ in the square interior (the bulk energy dominates the elastic energy for large $\tilde{\lambda}$). Further, the values $M=9.7\times 10^{-5}$, $N=1.78\times 10^{-3}$ are quite small and in this respect, we expect $\mathbf{Q}_{c}$ to be tailored by $\mathbf{Q}_{g}$, i.e., the coupling term promotes co-alignment of $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ and the coupling term dominates the bulk potential of $\mathbf{Q}_{c}$ (parameterised by $M$ and $N$). As such, we largely recover the $\mathbf{Q}_{g}$ profiles as reported in the reduced LdG study (with $D=0$) in [30]. The $\mathbf{n}_{c}$ profiles follow the $\mathbf{n}_{g}$ profiles with enhanced interior cybotactic clustering tailored by $S_{c}^{*}$. For example, we find the following solutions of (16) or critical points of (15): the T solution which is disordered along a single diagonal line; the H solution, featuring a +1/2 and a -1/2 defect; the J solution, which has lower order near the top boundary than R; R (rotated solution) for which $\mathbf{n}_{g}$ rotates by $\pi$ radians between a pair of opposite square edges (Fig. 5(a)). We follow the terminology used in [30]. The stability of a critical point can be measured in terms of its Morse index [20], i.e., the number of negative eigenvalues of the second variation of the generalised energy (15) evaluated at the critical point or equivalently, the number of unstable eigendirections for a critical point, $(\mathbf{Q}_{g},\mathbf{Q}_{c})$ of (15). The index of these solutions is almost the same as in the reduced LdG model with $D=0$ except that the index of the WORS is $5$ with $D=0$. When $D=0$, the WORS has an unstable eigenvector corresponding to a negative eigenvalue close to zero. This negative eigenvalue changes sign when $D>0$, making the WORS index-$4$ with weak coupling. The unstable critical points with positive Morse index are not experimentally observable but they play a crucial role in transition pathways between stable critical points and the selection of the stable solution for multistable systems. For example, index-$1$ unstable critical points are often referred to as transition states and can be experimentally observed during switching processes. Moreover, unstable critical points can be stabilised by external controls.For example, the unstable BD configuration for $D=E=0$ becomes stable with the constraint $q_{2}=p_{2}\equiv 0$ or under an applied electric field, as shown in Fig. 3.

We also use HiOSD dynamics [31] in Sec. V to construct pathways between the distinct critical points in Figure 5 with the WORS as the parent state. If we are able to find a critical point B by perturbing the critical point A, by following the HiOSD dynamics, we say that the critical points A and B are connected. These pathways/connections in Fig. 5(b) can provide guidance on how to effectively manipulate the defects, directors and cluster properties in confined BCNs.

In Fig. 5, the values of $M$ and $N$ are small so that the $\mathbf{Q}_{g}$ profiles tailor the $\mathbf{Q}_{c}$ profiles, i.e., $\mathbf{n}_{c}$ co-aligns with $\mathbf{n}_{g}$ in the solutions plotted in Figure 5. In Fig. 6, we select a negative value of $M$ (this further enhances the value of $S_{c}^{*}$) and take $M$ and $N$ to be $200$ times larger in magnitude than the values in Figure 5. We have no physical motivation for these choices but simply want to investigate the effects of $M$ and $N$ on the solution landscapes and how this can be used to yield greater autonomy to the $\mathbf{Q}_{c}$ profiles. This can necessarily deepen our understanding of the generalised model and its implications for BCN systems.

We have three notable observations. (i) As shown in Fig. 6(a), there are some $\mathbf{Q}_{g}$-solutions similar to those reported in [30]: index-$9$ WORS, index-$8$ C$\pm$, index-$6$ T, index-$5$ T*±*, index-$4$ H, index-$2$ BD, index-$1$ R, index-[math] D. We note that the increased values of $M$ and $N$ seem to increase the Morse index of the $\mathbf{Q}_{g}$ solutions compared to Figure 5. Interestingly, the rotated (R) solution—usually stable—is now destabilized in this regime. The $\mathbf{Q}_{c}$ profiles are more autonomous compared to Figure 5. Since the boundary condition for $\mathbf{Q}_{c}$ is set to zero, it is not constrained by the tangential anchoring condition. Therefore, compared with Fig. 5, the $\mathbf{Q}_{c}$ profile of the R state in Fig. 6 has less bending and closely resembles the $\mathbf{Q}_{c}$ profile of a BD solution, which is unstable in Figs. 5 and 6. We speculate this generates more unstable directions for the solution pair, $(\mathbf{Q}_{g},\mathbf{Q}_{c})$, and consequently, increases the corresponding Morse index. We also note that $\mathbf{n}_{c}$ and $\mathbf{n}_{g}$ are largely co-aligned for these examples.

(ii) There are at least three novel solution pairs in Fig. 6(b): the corresponding $\mathbf{Q}_{g}$ profiles are not attainable with $D=0$ or are not reported in the existing literature. Two novel solutions are connected to the WORS; the cross-shaped defect still exists but the centre is shifted along the diagonal or $x=0.5$. The third novel solution is connected to the T solution with an interior $-1/2$ defect. The solutions in Fig. 6(c) have familiar $\mathbf{Q}_{g}$ profiles but relatively novel $\mathbf{Q}_{c}$ profiles: the $\mathbf{n}_{c}$ profile tilts upwards to the right in the upper part, tilts downwards to the right in the lower part and is roughly coaligned with $\mathbf{n}_{g}$.

(iii) We surprisingly find a set of solutions labelled as index-$5$ BD*⟂, index-$3$ J⟂, index-$5$ J±⟂* and index-$1$ D*⟂*, for which the corresponding $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ are almost perpendicular to each other (Fig. 6(d)).

We take the WORS solution—the configuration with the highest Morse index in our manuscript—as an example to illustrate how $\mathbf{Q}_{c}$ generates more unstable directions. In Fig. 7(a), the first eigenvector drives the WORS solution towards the D solution, as the interior order is enhanced at the center with orientations along a diagonal. The second eigenvector drives WORS towards the R, BD, and H states, where the perturbation enhances the order at the center with directions parallel to one edge of the square domain. The third and fourth eigenvectors are degenerate, corresponding to the formation of the $T$ and $T_{\pm}$ states. Their linear combinations strengthen the order on both sides of the diagonal: one side is oriented nearly perpendicular to the diagonal while the other is almost parallel to the diagonal. The first four eigenvectors of the index‑9 WORS are similar to those of the index‑4 WORS shown in Fig. 7(b). Index-9 WORS has five more unstable directions than index-4 WORS. The combination of the seventh and eighth degenerate unstable directions leads to the index-6 and index-7 critical states in Fig. 6(b). The ninth unstable direction corresponds to the $C_{\pm}$ states, where the order is enhanced midway along the lines connecting the midpoint of the square to its four vertices, driving the formation of the $\pm$1 defects at the center. For the eigenvectors discussed above, the spatial distributions of $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$, $S_{g}$ and $S_{c}$ in the $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ profiles are consistent with each other. However, for the fifth and sixth eigenvectors, the orientations in the $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ profiles are almost perpendicular to each other, indicating the influence of $\mathbf{Q}_{c}$ on the instability of solutions. The fifth unstable direction corresponds to the $D_{\perp}$ states, while the sixth unstable direction corresponds to the $BD_{\perp}$, $J_{\perp}$, and $J_{\pm\perp}$ states in Fig. 6(d).

Here, we provide a heuristic explanation for the emergence of critical points with orthogonal $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$; see (Figs. 6 and 8). In the absence of defects with non-zero $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$, the variables $(q_{1},q_{2},p_{1},p_{2})$ can be written as:

[TABLE]

where $s=\sqrt{q_{1}^{2}+q_{2}^{2}}\geq 0$ and $p=\sqrt{p_{1}^{2}+p_{2}^{2}}\geq 0$. The relations between $S_{g},S_{c}$ in (14) and $s,p$ are $s=S_{g}/\sqrt{2}$ and $p=S_{c}/\sqrt{2}$. Substituting the expressions in (19) into (15), the energy can be written in terms of $(s,p,\theta,\phi)$ as given below:

[TABLE]

The corresponding Euler-Lagrange equations are:

[TABLE]

[TABLE]

[TABLE]

We define a defect-free region $\Omega_{\parallel\perp}$ (see Fig. 8) where $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ can be either parallel or perpendicular to each other, in paired critical points, i.e.,

[TABLE]

In the following, we take the diagonal states D (with co-aligned $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ in $\Omega_{\parallel\perp}$) and the paired critical state D*⟂* (with $\mathbf{n}_{g}\cdot\mathbf{n}_{c}=0$ in $\Omega_{\parallel\perp}$) in Fig. 8, the paired critical BD and BD*⟂* in Fig. 6 as examples.

We first note that, for both types of critical points—those with parallel and those with orthogonal $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$—the conditions (23) and (24) can be satisfied at least on $\Omega_{\parallel\perp}$. For D and D*⟂, the associated $\Omega_{\parallel\perp}$ is located around the square centre in Figure 8. For BD and BD⟂*, the corresponding $\Omega_{\parallel\perp}$ is also away from the square edges, the defect lines and localised near the square centre. Substituting $\sin(2(\theta-\phi))=0$ into the EL equations above, the equations (23)-(24) reduce to

[TABLE]

when $p>0$ and $s>0$. Due to the symmetry of D and D*⟂*,

[TABLE]

we have $\Delta\theta|_{(0.5,0.5)}=\Delta\phi|_{(0.5,0.5)}=0$ and we deduce that $\Delta\theta|_{\Omega_{\parallel\perp}}\approx 0$ and $\Delta\phi|_{\Omega_{\parallel\perp}}\approx 0$. Similarly, $\theta$ and $\phi$ are almost constant around the square centre for the paired critical points: BD and BD⟂ and hence, we deduce that $\Delta\theta|_{\Omega_{\parallel\perp}}\approx 0$ and $\Delta\phi|_{\Omega_{\parallel\perp}}\approx 0$.

As the re-scaled domain size $\tilde{\lambda}^{2}\to\infty$, the polynomial terms dominate in (21)-(22) and we require (to leading order) that

[TABLE]

Substituting $\phi=\theta+\frac{k\pi}{2}$ into the above equations, we have

[TABLE]

If (27) and (28) admit at least one positive solution pair with $s,p>0$, then we may have a critical point $(\mathbf{Q}_{g},\mathbf{Q}_{c})$ with approximately co-aligned $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$. If the algebraic equations (29) and (30) admit a solution with positive $s$ and $p$, then we may have a critical point with almost orthogonal GS and ES directors: $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$. Note that (27) and (28) are the same as (8) and (9) with $B=E=0$ and $S_{g}=\sqrt{2}s$, $S_{c}=\sqrt{2}p$. As in the previous discussion, we can derive a ninth-degree polynomial for $s$ or $p$ from (27)-(28) or (29)-(30). Since one of the solutions is $s=p=0$, the degree can be reduced to eight. Although, one can numerically solve the eighth-degree polynomial with specific parameter values and check for the existence of a positive solution with $s>0$ and $p>0$, it is difficult to analyse roots of an eighth-degree polynomial systematically.

Therefore, we assume that equations (27) and (29) ((28) and (30)) are cubic polynomials in terms of $s$ ($p$). Then the question can be partially answered by studying the number of positive roots of the following cubic equation:

[TABLE]

and how the number and sign of the roots depend on $a$, $c$ and $d$.

• •

Case 1: $a>0$, $c>0$

By using Cardano’s formula [13], since the discriminant of the cubic is given by

[TABLE]

we have one real solution and two complex solutions. The real solution is given by

[TABLE]

and

[TABLE]

• •

Case 2: $a>0$, $c<0$

If we choose $a>0$, $c<0$ and $d$ such that

[TABLE]

then we have three degenerate or non-degenerate real roots given by [23]

[TABLE]

where $\theta=arccos\left(\frac{3d}{2c}\sqrt{-\frac{c}{3a}}\right)$ and

[TABLE]

These arguments can be mapped to the parameter values in Figures 5 and 6. The choice of $A=-0.04$, $C=0.9$, $D=4\times 10^{-3}$ in (27) and (29) corresponds to Case 2 above. When $M=9.7\times 10^{-5}$, $N=1.78\times 10^{-3}$, $D=4\times 10^{-3}$, then the polynomials in (28) and (30) correspond to Case 1. Hence, we have a positive solution pair, $(s,p)$, only if $k$ (in $\phi=\theta+\frac{k\pi}{2}$) is even in (27)-(28), i.e., the directors $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ are almost co-aligned on $\Omega_{\parallel\perp}$. When $M=-200\times(9.7\times 10^{-5})$, $N=200\times(1.78\times 10^{-3})$, $D=4\times 10^{-3}$ as in Figure 6, the polynomials in (28) and (30) correspond to Case 2 and we always have at least one positive solution $s,p>0$, for both (27)-(28) and (29)-(30). Hence, the directors $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ can be either co-aligned or perpendicular on $\Omega_{\parallel\perp}$, yielding paired critical points $(\mathbf{Q}_{g},\mathbf{Q}_{c})$.

To summarise, we have investigated solution landscapes of prototype BCN systems within the generalised model in (5), i.e., we have studied critical points $(\mathbf{Q}_{g},\mathbf{Q}_{c})$ of the free energy in (5). We have focused on a two-dimensional example of BCNs on square domains, subject to experimentally relevant tangent boundary conditions for $\mathbf{Q}_{g}$. We impose isotropic boundary conditions on $\mathbf{Q}_{c}$ but other choices are possible and it is reasonable to assume that cybotactic clusters are confined to the domain interior. The essential novelty is introduced by the coupling between $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$. This coupling promotes interior ordering and stabilises ordered GS profiles even for high temperatures, when we only expect predominantly isotropic solutions at $D=0$ for high temperatures. In fact, for $A=0.04$, the GS-ES coupling and/or external fields can stabilise the diagonal solution, which is typically only found for low temperatures. The choice of boundary conditions can also stabilise symmetric WORS-type solutions, i.e., the coupling parameter $D$ can compete with the disordering effects of high temperature. For relatively small values of $M$ and $N$, the $\mathbf{Q}_{c}$ profiles are tailored by the $\mathbf{Q}_{g}$ profiles and there are few surprises for low temperatures. For larger values of $M$ and $N$, the $\mathbf{Q}_{c}$ profiles can be more independent of the $\mathbf{Q}_{g}$ profiles leading to non-trivial and novel effects. For example, we find novel solution pairs or novel critical points in Figure 6 which are not found with $D=0$. The Morse index of critical points of (5) typically increases with increasing values of $M$ and $N$, since there are new unstable directions determined by the relatively autonomous $\mathbf{Q}_{c}$ profiles. For example, the rotated $R$ critical point has index-[math] in Figure 5 and has index-$1$ in Figure 6. Importantly, we recover critical points $(\mathbf{Q}_{g},\mathbf{Q}_{c})$ with almost orthogonal $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ in Figure 6. Whilst we cannot make things precise, it is possible that the optical measurements are determined by an appropriate combination of $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ within the generalised model in (5). If $\mathbf{n}_{g}$ and $\mathbf{n}_{c}$ are orthogonal to each other in experimentally attainable BCN profiles, then this could lead to experimentally observable macroscopic biaxiality. Of course, such arguments need much more investigation than is possible in this work.

The domains with $\tilde{\lambda}^{2}=2500$ and $5000$ correspond to experimentally relevant micron-scale areas of $(0.26\mu m)^{2}$ and $(0.52\mu m)^{2}$, respectively. As the domain size $\lambda^{2}$ increases, the defects become smaller, requiring a finer computational grid. This, in turn, leads to a higher computational cost, lower numerical accuracy, and greater difficulty in achieving convergence. As discussed in the context of the LdG theory for nematics in [30] and [25], computations on small domains remain relevant for large domains. The solutions of the Euler-Lagrange equations on small domains serve as excellent initial conditions for numerical continuation schemes, for large values of $\lambda^{2}$. Most of the solution branches, for small domains, continue to exist on larger domains although the stability properties might change, i.e., some solution branches lose stability as domain size increases. However, some solution branches retain stability or gain stability, as domain size increases or $\lambda^{2}$ increases. The unstable solutions also play a crucial role in the system dynamics and the selection of stable solutions for multistable systems, and hence, our results for $\lambda^{2}=2500$ (although not large enough for experiments on the several-hundred-micrometer scale), can still provide direct/indirect guidance for future experimental studies.

V Numerical methods

The numerical results of Section IV, including the minimisers and critical points of the non-dimensionalised free energy in 2D setting $\tilde{F}$ in (15), are calculated by the numerical methods as follows. The study is implemented in Matlab R2021b .

V.1 Spatial discrete method

We use the finite difference method to estimate the spatial derivative on a 2D square $\tilde{\Omega}=[0,1]^{2}$ taking the nodes $(x_{i},y_{j})$, $i,j=0,1,\cdots,N_{0}$ with step length $h=\frac{1}{N_{0}}$, where

[TABLE]

The five-point stencil method is used to approximate the Laplacian $\nabla^{2}\hat{\mathbf{x}}$ of a function $\hat{\mathbf{x}}(x,y)$ at a point $(x_{i},y_{j})$ on the discrete grid with second-order accuracy:

[TABLE]

V.2 HiOSD method

We use the HiOSD method [31], which is a generalization of the optimization-based shrinking dimer method [32] to find stable and unstable solutions of the Euler-Lagrange equations in (16) on the energy landscape of $\tilde{F}$ in (15).

The stability of a critical point of the free energy can be measured by its Morse index [20]. For a non-degenerate index-$k$ saddle point $\hat{\mathbf{x}}=(\mathbf{Q}_{g},\mathbf{Q}_{c})$, the Hessian $\mathbb{H}(\hat{\mathbf{x}})$ has exactly $k$ negative eigenvalues $\hat{\lambda}_{1}\leqslant\cdots\leqslant\hat{\lambda}_{k}$ with corresponding unit eigenvectors $\hat{\mathbf{v}}_{1},\cdots,\hat{\mathbf{v}}_{k}$ satisfying $\big\langle\hat{\mathbf{v}}_{j},\hat{\mathbf{v}}_{i}\big\rangle=\delta_{ij}$, $1\leqslant i,j\leqslant k$. Define a $k$-dimensional subspace $\hat{\mathcal{V}}=\mathrm{span}\big\{\hat{\mathbf{v}}_{1},\cdots,\hat{\mathbf{v}}_{k}\big\}$, then $\hat{\mathbf{x}}$ is a local maximum on a $k$-dimensional linear manifold $\hat{\mathbf{x}}+\hat{\mathcal{V}}$ and a local minimum on $\hat{\mathbf{x}}+\hat{\mathcal{V}}^{\perp}$, where $\hat{\mathcal{V}}^{\perp}$ is the orthogonal complement space of $\hat{\mathcal{V}}$. In particular, a stable state is an index-[math] solution.

An index-$k$ critical point can be found by the $k$-HiOSD dynamics with a certain initial condition. A $k$-HiOSD dynamics is a transformed gradient flow of $\mathbf{x}$ coupling with the search for an orthonormal basis $\mathcal{V}$ which minimises $k$ Rayleigh quotients simultaneously.

[TABLE]

To avoid evaluating the Hessian of $\tilde{F}(\mathbf{x})$, we use central difference schemes for directional derivatives to approximate Hessians along $i$th dimer with length $2l$ centred at $\mathbf{x}$,

[TABLE]

We use the time-discrete Euler scheme with $\Delta t=10^{-4}$.

We obtain a critical point of the free energy if the HiOSD dynamics converges with error tolerance $|\nabla\tilde{F}(\hat{\mathbf{x}})|\leq 10^{-6}$. The index of the critical point is checked by the calculation of the smallest $k+1$ eigenvalues of the Hessian $\nabla^{2}\tilde{F}(\hat{\mathbf{x}})$ by using the simultaneous Rayleigh-quotient minimisation method [15].

V.3 Algorithm for building the connectivity of critical solutions

Following the HiOSD dynamics, we build the connectivity of critical points using two algorithms: a downward search that enables us to search for all connected lower-index saddles from an index-$k$ saddle; an upward search with a selected direction to find the higher-index saddles. The two algorithms drive the entire search to navigate up and down on the energy landscape. Please refer to [9] for detailed procedure of downward and upward search.

VI Conclusion

In this paper, we propose and study a generalised LdG-type model for BCN, building on the two-state model proposed by Madhusudana in [17]. Our model is based on two tensor order parameters: $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ instead of two scalar order parameters, $S_{g}$ and $S_{c}$ as in [17]. Our model has all the capabilities of Madhusudana’s model and more — it can account for the GS director and additional directors induced by the cybotactic clusters, spatial inhomogeneities, effects of boundary conditions, geometric frustration, defects etc. Notably, we find that the GS-ES coupling, measured in terms of $D$ in (5), can induce phase transitions at a fixed temperature. For example, we can reproduce the nematic-paranematic phase transition and the paranematic-isotropic phase transitions by reducing $D$, at a fixed high temperature. Madhusudana reports the same sequence of phase transitions with increasing temperature at fixed $D$, in 3D. We also study phase transitions in 2D although strictly 2D systems are not realistic. However, $\mathbf{Q}_{g}$ and $\mathbf{Q}_{c}$ are reduced to symmetric traceless $2\times 2$ matrices in quasi-2D settings as in Section IV and as such, the 2D phase diagrams in Section III can shed light into planar ordering and disordering as a function of $D$ in such settings. The 2D phase diagrams in Section III certainly help us understand structural phase transitions in confined quasi-2D systems as a function of $D$. Of course, we do not have any physical insight into how to tune or manipulate $D$, except that it is an intrinsic material parameter.

Our model can be embellished to capture emergent chirality, transitions to the twist-bend and splay-bend phases and macroscopic biaxiality. This could be done by introducing elastic anisotropy into the generalised free energy, i.e., more general and higher-order terms (cubic terms such as $\mathbf{Q}_{g}\nabla\mathbf{Q}_{g}\nabla\mathbf{Q}_{g}$) in the elastic energy in (5). Inspiration can be sought from existing work on the critical role of elastic constants in the emergence of the twist-bend and splay-bend phases in BCN; see [11]. However, we have limited insight into typical values for $K_{c}$ — the elastic constant associated with $\mathbf{Q}_{c}$. We have employed a one-constant approximation in this paper, with $K_{g}=K_{c}$, but the precise form of the elastic energy of $\mathbf{Q}_{c}$ is an open question. Octupolar order could be introduced by adding a higher-order tensor such as $T_{ijk}$ with an appropriate physical interpretation [11] but open questions remain with regards to the associated energetic contributions.

It is commendable that a simple two order parameter model in [17], further generalised in our paper, is capable of capturing some complex experimental trends. This is, in itself, a strong reason for studying such models since they offer tractable approaches to complex multi-physics problem. However, they might miss subtle underpinning physics about the cybotactic clusters and can only offer a coarse-grained homogenised perspective. For example, our model does not have additional positional order parameters to account for density modulations within cybotactic clusters. This model does not account for the spacing between the cybotactic clusters or the Smectic-type ordering within the cybotactic clusters.To describe these additional features, one may need to adopt a model for SmC‑type ordering, such as [29], in which a density order parameter is introduced in addition to $\mathbf{Q}_{c}$, and a coupling term between the layer normal (which can be represented through the first‑ or second‑order derivatives of the density field) and $\mathbf{Q}_{c}$ is included in the free energy. This leads to a more complex mathematical formulation, which will be explored in our future work. Further, Madhusudana’s model and our generalised model can offer no insight into how the opening angle of BCN molecules affects the macroscopic BCN phases. We employ a generic coupling term proportional to $\mathbf{Q}_{g}\cdot\mathbf{Q}_{c}$ as an immediate generalisation of the GS-ES coupling term in [17]. This is a canonical coupling term that arises in the dilute limit of suspensions, e.g., when we treat the cybotactic clusters as a dilute suspension of smectic nanoparticles or nanoclusters in a nematic GS medium, consistent with the assumptions in Madhusdana’s paper [1]. Similar coupling terms are also employed in phenomenological studies of ferronematics — dilute suspensions of magnetic nanoparticles in nematic media [8]. This coupling term can miss nonlinear phenomena. For example, as the global nematic order $\mathbf{Q}_{g}$ increases, $\mathbf{Q}_{c}$ may exhibit nonlinear behavior such as saturation characteristics. In bent‑core nematics, the molecules are intrinsically bent and typically respond differently along the bending plane and the perpendicular plane. Possible improvements include introducing higher‑order or anisotropic couplings. Further, it is often difficult to decide on the exact form of the elastic energy associated with $\mathbf{Q}_{c}$ or the exact origins of the Ginzburg-Landau type potential for $\mathbf{Q}_{c}$ — what is the underlying physics? Finally, how do we estimate the parameters $M,N$ and $D$ in the free energy and relate these phenomenological parameters to physical experiments on BCNs in a reliable way? Optical measurements probably only yield a mean scalar order parameter as suggested in [17] but one needs quantitative information about $S_{c}$ to infer meaningful information about $M$, $N$ and $D$. Open questions remain about the interpretability of this model but it is undeniable that this simple two order parameter model captures complex nonlinear phenomena qualitatively and certainly gives insight into cybotactic cluster formation, location and stability in prototype confined systems, making this a worthwhile and interesting study.

Acknowledgements.

YH gratefully thanks Prof. Lei Zhang and Beijing International Center for Mathematical Research of Peking University for hosting her as a Visiting Scholar. PR acknowledges support from Leverhulme Research Project Grant RPG-2021-401. AM is supported by the University of Strathclyde New Professors Fund, the Humboldt Foundation and a Leverhulme Research Project Grant RPG-2021-401.