# A diffuse-interface Landau-de Gennes model for free-boundary nematic liquid
crystals

Dawei Wu  Baoming Shi  Yucen Han  Pingwen Zhang  Apala Majumdar  Lei Zhang

###### Abstract

We introduce a diffuse-interface Landau-de Gennes free energy for free-
boundary nematic liquid crystals (NLC) in three dimensions submerged in
isotropic liquid, where a phase field is introduced to model the deformable
interface. The energy we propose consists of the original Landau-de Gennes
free energy, three penalty terms and a volume constraint. We prove the
existence and regularity of minimizers to the diffuse-interface energy
functional. We also prove a uniform maximum principle of the minimizer under
appropriate assumptions, together with a uniqueness result for small domains.
Then, we establish a sharp-interface limit where minimizers of the diffuse-
interface energy converge to a minimizer of a sharp-interface energy under the
framework of $\Gamma$-convergence. Finally, we conduct numerical experiments
with the diffuse-interface model, the findings of which are compared with
existing works.

#### Keywords.

nematic liquid crystals, phase separation, Landau-de Gennes,
$\Gamma$-convergence

#### AMS subject classifications.

76A15, 49J27, 49K20, 35B50, 49J45, 35F30

##  1 Introduction

The liquid crystal (LC) is an intermediate state between solids and liquids,
in which positional order of the crystal is partially or completely lost, but
molecular anisotropy is still present [9]. The simplest case of liquid
crystals is the nematic phase in which the constituent molecules have no
positional order, but prefer to align locally along certain directions known
as directors. Nematic liquid crystals (NLC) are endowed with direction-
dependent physical, optical and rheological properties [53]. NLC’s exhibit
intriguing morphology because of their molecular anisotropy, especially when a
droplet of undetermined shape interacts with a different material, about which
people have conducted many experiments. NLC dispersed in polymer (often called
polymer-dispersed liquid crystal in the community) forms droplets in various
shapes and sizes, leading to peculiar optoelectric properties when sandwiched
among planar electrodes [62, 13]. A special droplet morphology of interest is
the nematic tactoid, i.e. spindle-shaped droplets filled with NLC whose
directors align with the surface, which tends to emerge during the nematic-
isotropic phase transition induced by temperature [44, 57, 30] or around rod-
shaped bacteria as protective sheaths [56]. The myriad phenomena of free-
boundary NLC droplets have inspired widespread applications. For instance,
polymer-dispersed LC have been produced with simple and cost-effective phase
separation methods, leading to huge advances in the optoelectric display
industry [13]. Surface grating patterns can be carefully tailored to obtain
bistable nematic states, promising the production of low-energy, easily
switchable LC displays [31, 54, 67]. With the aid of NLC morphology, people
have also accurately simulated polymerized nematic elastic ribbons with
multiple stable shapes [48], biomimetic nematic droplets with non-equilibrium
oscillation [29] and the morphogenesis of living cells [35, 40]. These
examples demonstrate the vast potentials of free-boundary NLC.

To study free-boundary NLC mathematically, we need to construct a model in two
aspects: one about free boundary and the other about NLCs. For modelling NLCs,
there is the Onsager (molecular) model based on local director distributions,
the Ossen-Frank (vector) based on the unit director field, the Landau-de
Gennes (LdG, tensor) model based on moments of the director distribution and
many others [21, 60]. Each model has its strengths and weaknesses, but all of
them have been successful in predicting the behaviour of NLC, see e.g. [65]
(Onsager), [46] (Oseen-Frank) and [33] (LdG). For modelling free boundary, we
quote existing methods from the field of shape and topology optimization [12].
Among the most popular are the surface mesh method [3] that discretizes the
free boundary with a finite-element mesh known as the interface segregating
different phases, and the diffuse-interface method [61] that assumes a
continuous variation of the phase within the interfacial regions described by
a smooth field (space-dependent function) known as the phase field. In recent
decades, there has been a multitude of research in free-boundary NLCs with
fairly diverse modelling approaches, both in free boundary and in NLC [28, 66,
36, 16, 63, 30, 41, 11, 20, 2, 37]. The common practice is to design energy
functionals with respect to the shape and the NLC order parameter, and find
the stable configuration with energy minimization methods. For instance, in
[2] the authors describe a deformable two-dimensional NLC confinement with a
finite-element mesh, and their energy functional consists of the Landau-de
Gennes free energy and a penalty function for Dirichlet boundary values; they
discover that with a strong penalty, the shape of the confinement converged
into a tactoid through numerical simulation, which is consistent with
experimental findings. The authors of [11] study the same problem with the
finite element method (FEM) as well, but they choose the Oseen-Frank (vector)
free energy instead and further penalize the energy functional with the
perimeter of the confinement; they also find tactoids when the penalty factors
are large. In [66] the authors use a phase field $\phi$ to denote separation
of NLC and a Newtonian fluid with $\phi=1$ for NLC and $\phi=-1$ for fluid,
and their energy functional consists of a boundary anchoring with $\nabla\phi$
substituted for the boundary normal, a mixing energy of van der Waals-Cahn-
Hilliard form [4] to penalize the sharp phase separation ($\phi\approx\pm 1$)
and smoothness of $\phi$, and finally the Oseen-Frank free energy density
masked by $\frac{\phi+1}{2}$ so that it is only integrated over the NLC phase;
numerical tests of their energy minimization find multiple stable states,
including one radially symmetric state and one slightly elliptical bipolar
state reminiscent of the tactoid. In general, the application of
shape/topology optimization methods to free-boundary NLC has been highly
successful in explaining their complex anisotropy-induced morphology, and
promises an immense area of exploration both numerically and theoretically.

We propose a diffuse-interface model for a free-boundary nematic phase in
three dimensions (3D) submerged in the isotropic phase, using the Landau-de
Gennes free energy for the modelling of liquid crystals and a phase field for
deformable shapes. Our approach is mainly inspired by [66, 17, 8]. The model
is relatively new since NO work as far as we know combines LdG theory with
diffuse-interface models to study free-boundary NLCs. We choose the LdG model
to describe NLC because the $\mathbf{Q}$-tensor order parameter has fewer
degrees of freedom than the actual director distribution, so the computation
is less intensive than the Onsager model. The LdG model can also represent
$+1/2$ and $-1/2$ defects in two dimensions (a discontinuity of the director
field around which the director rotates $+1/2$ or $-1/2$ cycles, quite common
in real-life NLCs [32, 59]), where the director field is non-orientable and
thus impossible for the Oseen-Frank model. Theoretical properties such as
maximum principle and uniaxial symmetry are also available [39, 25, 6], which
mathematically guarantee the feasibility of the LdG model. With the LdG model,
one can also produce a plethora of complicated energy landscapes over
confinements of different shapes, see e.g. [64, 23, 22, 24, 51, 50].
Therefore, the LdG model possesses huge potentials in the free boundary
problem. We choose the diffuse-interface model to describe the deformable
shape because the phase field can be defined on a fixed grid [52, 61] while a
surface mesh requires a full finite-element mesh [3]. The degrees of freedom
are therefore considerably lower for the phase field, since it only involves a
scalar function, but the finite-element mesh requires the spatial coordinates
of all nodes and the adjacency relations between them. Moreover, the boundary
of the deformable shape is identified by variations of the phase field, so
topological changes can be implemented, e.g. one region splitting into two or
two regions merging into one [66, 35, 7]. On the contrary, the surface mesh
cannot describe topological changes because the node adjacency relations are
fixed.

Our work is presented in the following order. In Section 2, we define a
diffuse-interface energy functional $E_{\varepsilon}$ made up of the classic
LdG free energy in 3D defined by the $\mathbf{Q}$-tensor, penalty terms and a
volume constraint on the phase field $\phi$. There are three extra competing
energy terms besides the LdG energy. The volume constraint ensures the mass
conservation of the nematic phase. There is also a hyperparameter
$\varepsilon>0$ that signifies the width of the diffuse interface. In Section
3, we establish the solvability, i.e. existence of minimizers of
$E_{\varepsilon}$, with classical methods of calculus of variation [14]. In
Section 4, we discuss the problem of minimizing $E_{\varepsilon}$ on a fixed
phase field. We discuss the maximum principle with the general theory of
elliptic equations, and then derive uniqueness of the minimizer based on a
convexity argument from [34]. Our discoveries correspond to well-known results
of the classical LdG model on a fixed region. In Section 5, we present our
main contribution — the limit of $E_{\varepsilon}$ as $\varepsilon\to 0,$
known as the sharp-interface limit, with $\Gamma$-convergence [10] as the
primary tool. We prove under appropriate assumptions the convergence of
minimizers of $E_{\varepsilon}$ to the minimizer of a sharp-interface
functional $E_{0}$ defined over deformable regions that also models the
separation of nematic and isotropic phases. In Section 6, we conduct numerical
experiments on the minimization of a reduced version of $E_{\varepsilon}$ on a
two-dimensional domain. We adjust hyperparameters to observe changes in the
minimizer configuration, and compare them with past findings.

##  2 Diffuse-interface Landau-de Gennes model

Let $\Omega\subset\mathbb{R}^{3}$ be a region with Lipschitz boundary, where a
blob of nematic liquid crystals is surrounded by isotropic liquid. The
separating interface is called the nematic-isotropic (N-I) interface. We
propose a diffuse-interface LdG energy

(2.1) $\displaystyle E_{\varepsilon}=E^{\rm LdG}+\omega_{p}E^{\rm
mix}_{\varepsilon}+\omega_{a}E^{\rm anch}_{\varepsilon}+\omega_{v}E^{\rm
void},$

where $\omega_{p},\omega_{v},\omega_{a}$ are positive weights of competitive
energy terms and $\varepsilon>0$ is a positive parameter appearing in the
second and third terms. The units of $\omega_{p},\omega_{v},\omega_{a}$ are
$\mathrm{N\cdot m^{-1}}$, $\mathrm{N\cdot m^{-1}}$ and $\mathrm{N\cdot
m^{-2}}$ respectively. The diffuse-interface LdG energy is composed of four
types of free energy: (a) LdG free energy of nematic liquid crystals, (b)
mixing energy of the N-I interface, (c) tangential anchoring energy of the
liquid crystal molecules on the interface and (d) void energy (penalty) of the
isotropic phase.

#### (a) LdG energy

To describe the director and corresponding orientational order of the nematic
phase, we introduce the tensor variable $\mathbf{Q}=(Q^{ij})_{3\times 3}$, a
$3\times 3$ symmetric traceless matrix encoding information about nematic
director distribution [9]. $\mathbf{Q}=\mathbf{0}$ represent the isotropic
phase; if only two of the eigenvalues are the same, $\mathbf{Q}$ is uniaxial;
if the eigenvalues are distinct, $\mathbf{Q}$ is biaxial. A uniaxial state,
which can be represented by

$\mathbf{Q}=s\quantity(\mathbf{n}\otimes\mathbf{n}-\frac{1}{3}\mathbf{I}),$

with $s\in\mathbb{R}$ and $\mathbf{n}\in\mathbb{S}^{2}$, is the most
interesting because $\mathbf{n}$ recovers the nematic director, i.e. local
preferred molecular alignment. The term $E^{\rm LdG}$ in (2.1) is the
classical LdG free energy

(2.2) $E^{\rm
LdG}[\mathbf{Q}]=\int_{\Omega}(F_{el}+F_{b})\mathop{}\\!\mathrm{d}x,$

with elastic energy density

(2.3) $F_{el}=\frac{L}{2}|\gradient\mathbf{Q}|^{2},$

where $L>0$ is the material-dependent elastic constant,
$|\gradient\mathbf{Q}|^{2}\triangleq\sum_{i,j,k}|Q^{ij}_{x_{k}}|^{2}$, and
bulk energy density

(2.4)
$F_{b}(\mathbf{Q})=\frac{A}{2}\tr\mathbf{Q}^{2}-\frac{B}{3}\tr\mathbf{Q}^{3}+\frac{C}{4}\quantity(\tr\mathbf{Q}^{2})^{2},$

where $A<0$ is a rescaled temperature that can be tuned to achieve the
nematic-isotropic phase transition (we assume that $A<0$ to stabilize the
nematic phase), and $B,C>0$ are material-dependent constants. The units of
$A,B,C$ are $\mathrm{N\cdot m^{-2}}$, and the unit of $L$ is $\mathrm{N}$. As
a reminder, the elastic energy density (2.3) adopts the simplest one-constant
form, but more general forms exist [43, 21]. Under our setup, the minimizers
of the functional $F_{b}$ (2.4) are uniaxial states
$\mathbf{Q}^{*}=s_{+}\quantity(\mathbf{n}\otimes\mathbf{n}-\frac{1}{3}\mathbf{I})$,
with arbitrary unit vector $\mathbf{n}\in\mathbb{S}^{2}$ and

(2.5) $s_{+}=\frac{B+\sqrt{B^{2}-24AC}}{4C}.$

#### (b) Mixing energy

The term $E^{\rm mix}_{\varepsilon}$ in (2.1) is the mixing energy. It is the
well-known van der Waals-Cahn-Hilliard energy functional [4, 58] describing
binary mixture of two phases, which in this case are nematic and isotropic
liquid crystals.

(2.6) $E^{\rm
mix}_{\varepsilon}[\phi]=\int_{\Omega}\left[\varepsilon|\nabla\phi|^{2}+\varepsilon^{-1}\phi^{2}(1-\phi)^{2}\right]\mathop{}\\!\mathrm{d}x,$

where $\phi$ is a phase field. Informally speaking, the continuous double-well
function $\phi^{2}(1-\phi)^{2}$ prefers the separation of $\Omega$ into
subdomains with $\phi\approx 1$ for nematics and $\phi\approx 0$ for
isotropic, while the gradient term $|\nabla\phi|^{2}$ smoothens out $\phi$ to
form a diffuse interface, whose thickness is represented by $\varepsilon>0$,
the capillary width with unit $\mathrm{m}$ (see Figure 1). The van der Waals-
Cahn-Hilliard energy is known to approximate the perimeter of the interface in
the sense of $\Gamma$-convergence as $\varepsilon\to 0$ (the sharp interface
limit) [42], so it models the interfacial tension, as is common in numerous
physical models [27, 38].

$\phi=1$$\phi=0$$O(\varepsilon)$$\partial D$ Figure 1: Schematic figure of N-I
mixing. As colour change from bright to dark (greyish-red online), phase
changes from nematic ($\phi=1$) to isotropic ($\phi=0$) via diffuse N-I
interface with thickness $O(\varepsilon)$ surrounding the surface $\partial
D$.

#### (c) Anchoring energy

The term $E^{\rm anch}$ in (2.1), aimed to enforce tangential anchoring [9,
59] on the N-I interface, is given by

(2.7) $E^{\rm
anch}_{\varepsilon}[\mathbf{Q},\phi]=\int_{\Omega}\varepsilon\left|\left(\mathbf{Q}(x)+\frac{s_{+}}{3}\mathbf{I}\right)\nabla\phi\right|^{2}\mathop{}\\!\mathrm{d}x,$

with $s_{+}$ defined in (2.5). $E^{\rm anch}_{\varepsilon}$ is minimized when
$\mathbf{Q}\nabla\phi=-\frac{s_{+}}{3}\nabla\phi$, i.e. $\nabla\phi$ is an
eigenvector of $\mathbf{Q}$ corresponding to the negative eigenvalue
$-\frac{s_{+}}{3}$, indicating that the nematic director of $\mathbf{Q}$ is
perpendicular to $\nabla\phi$, which approximates the normal direction within
the diffuse interface. To keep this energy at $O(1)$ order in the sharp-
interface limit $\varepsilon\to 0$, the factor $\varepsilon$ is needed in
(2.7). One can verify this fact by checking from Figure 1 that
$|\nabla\phi|=O(\varepsilon^{-1})$ and that the size of its support is
$O(\varepsilon)$, making the integral without $\varepsilon$ in (2.7)
approximately $O(\varepsilon^{-2})\cdot O(\varepsilon)=O(\varepsilon^{-1})$.
This weak anchoring has been applied in works such as [26] and [49].

#### (d) Void energy

The term $E^{\rm void}$ in (2.1), aimed to penalize the isotropic phase
($\phi\approx 0$) surrounding the nematic phase, is given by

(2.8) $E^{\rm
void}[\mathbf{Q},\phi]=\int_{\Omega}\frac{1}{2}(1-\phi)^{2}|\mathbf{Q}|^{2}\mathop{}\\!\mathrm{d}x,$

where $|\mathbf{Q}|^{2}=\sum_{i,j}|Q^{ij}|^{2}$ is the Frobenius norm. It is
inspired by similar penalty terms in the studies of phase separation in fluids
[17, 7].

We impose the volume constraint

$\int_{\Omega}\phi\mathop{}\\!\mathrm{d}x=V_{0},$

to conserve the mass of NLCs with $0<V_{0}<|\Omega|$ ($|\Omega|$ is the volume
of $\Omega$), and impose homogeneous Dirichlet boundary conditions

$\mathbf{Q}|_{\partial\Omega}=\mathbf{0},\quad\phi|_{\partial\Omega}=0,$

to model the isotropic phase surrounding the nematic phase.

The energy (2.1) can be nondimensionalized. Let $\bar{x}=\lambda^{-1}x$ for
all $x\in\Omega$. After substituting $\bar{x}$ in the integrals and dividing
the entire energy by $\lambda L$, we obtain

(2.9) $\displaystyle\overline{E}_{\bar{\varepsilon}}[\mathbf{Q},\phi]$
$\displaystyle=\bar{E}^{\rm
LdG}+\bar{\omega}_{p}\bar{\lambda}\bar{E}_{\bar{\varepsilon}}^{\rm
mix}+\bar{\omega}_{a}\bar{\lambda}\bar{E}_{\bar{\varepsilon}}^{\rm
anch}+\bar{\omega}_{p}\bar{\lambda}^{2}\bar{E}^{\rm void}$
$\displaystyle=\int_{\Omega^{(0)}}\quantity[\frac{1}{2}|\gradient\mathbf{Q}|^{2}+\bar{\lambda}^{2}\bar{F}_{b}(\mathbf{Q})]\mathop{}\\!\mathrm{d}\bar{x}+\bar{\omega}_{p}\bar{\lambda}\int_{\Omega^{(0)}}\quantity[\bar{\varepsilon}|\nabla\phi|^{2}+\frac{1}{\bar{\varepsilon}}\phi^{2}(1-\phi)^{2}]\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\bar{\omega}_{a}\bar{\lambda}\int_{\Omega^{(0)}}\bar{\varepsilon}\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\nabla\phi\right|^{2}\mathop{}\\!\mathrm{d}\bar{x}+\bar{\omega}_{v}\bar{\lambda}^{2}\int_{\Omega^{(0)}}\frac{1}{2}(1-\phi)^{2}|\mathbf{Q}|^{2}\mathop{}\\!\mathrm{d}\bar{x}$

with nondimensional domain $\Omega^{(0)}$ and rescaled nondimensional
constants

$\bar{\lambda}^{2}\triangleq\frac{\lambda^{2}C}{L},\
\bar{\varepsilon}\triangleq\frac{\varepsilon}{\lambda},\
\bar{\omega}_{p}\triangleq\frac{\omega_{p}}{\sqrt{CL}},\
\bar{\omega}_{v}\triangleq\frac{\omega_{v}}{C},\
\bar{\omega}_{a}\triangleq\frac{\omega_{a}}{\sqrt{CL}},$

and rescaled LdG bulk energy $\bar{F}_{b}\triangleq\frac{1}{C}F_{b}$, subject
to the nondimensional volume constraint

(2.10) $\int_{\Omega^{(0)}}\phi\mathop{}\\!\mathrm{d}\bar{x}=\bar{V}_{0},$

with $\bar{V}_{0}=\lambda^{-3}V_{0}$ satisfying
$0<\bar{V}_{0}<|\Omega^{(0)}|$.

From this point onwards, we will focus on the nondimensional functional
$\overline{E}_{\bar{\varepsilon}}$ (2.9), which is almost identical to (2.1)
except for the additional $\bar{\lambda}$ in the leading coefficients and the
bars above. For notational simplicity, we rewrite $\Omega^{(0)}$ as $\Omega$
and drop bars over all constants from (2.9) (note: the bar over $\bar{F}_{b}$
is kept in order to differ from the classical bulk energy $F_{b}$ (2.4)). The
four terms of (2.9) are identified with the LdG energy $E^{\rm LdG}$ (2.2),
the mixing energy $E^{\rm mix}_{\varepsilon}$ (2.6), the anchoring energy
$E^{\rm anch}_{\varepsilon}$ (2.7) and the void energy $E^{\rm void}$ (2.8)
respectively. We also denote the admissible space of the order parameters
$(\mathbf{Q},\phi)$ by

(2.11) $\mathscr{A}=\quantity{(\mathbf{Q},\phi):\mathbf{Q}\in
H_{0}^{1}(\Omega;\mathcal{S}_{0}),\ \phi\in H_{0}^{1}(\Omega),\
\int_{\Omega}\phi\mathop{}\\!\mathrm{d}x=V_{0}}$

where we use $H_{0}^{1}$ to denote the Sobolev space [18, 1, 14]

$H_{0}^{1}(\Omega)=\quantity{u\in L^{2}(\Omega):\int_{\Omega}[u^{2}+|\nabla
u|^{2}]\mathop{}\\!\mathrm{d}x<\infty,\ u|_{\partial\Omega}=0},$

$H_{0}^{1}(\Omega;V)$ to denote the Sobolev space with function values in a
normed vector space $V$, and also

(2.12) $\mathcal{S}_{0}=\\{\mathbf{Q}\in\mathbb{R}^{3\times
3}:\mathbf{Q}^{T}=\mathbf{Q},\tr\mathbf{Q}=0\\}$

to denote the space of symmetric traceless $3\times 3$ matrices.

If $(\mathbf{Q},\phi)$ is a critical point of the functional (2.9), then they
are weak solutions to their corresponding Euler–Lagrange equations:

(2.13) $\displaystyle\triangle\mathbf{Q}$
$\displaystyle=\lambda^{2}\quantity[\quantity(\frac{A}{C}+\tr\mathbf{Q}^{2})\mathbf{Q}-\frac{B}{C}\mathscr{P}_{0}(\mathbf{Q}^{2})]+\omega_{v}\lambda^{2}(1-\phi)^{2}\mathbf{Q}$
$\displaystyle\quad+\omega_{a}\lambda\varepsilon\mathscr{P}_{0}\quantity[(\nabla\phi\otimes\nabla\phi)\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})+\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})(\nabla\phi\otimes\nabla\phi)],$

and

(2.14) $\triangle\phi+\frac{\omega_{a}}{\omega_{p}}\operatorname{\rm
div}\quantity[\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})^{2}\nabla\phi]=\frac{1}{\varepsilon^{2}}\phi(\phi-1)(2\phi-1)+\frac{\omega_{v}\lambda|\mathbf{Q}|^{2}}{2\omega_{p}\varepsilon}(\phi-1)+\xi,$

where $\nabla_{\mathcal{S}_{0}}$ is the projected gradient operator on
$\mathcal{S}_{0}$ (2.12), which maps the matrix gradient onto
$\mathcal{S}_{0}$ with the projection operator

(2.15)
$\mathscr{P}_{0}(\mathbf{A})=\mathbf{A}-\frac{\tr\mathbf{A}}{3}\mathbf{I}$

to maintain the traceless constraint on $\mathbf{Q}$, and $\xi\in\mathbb{R}$
is a Lagrangian multiplier originating from the volume constraint (2.10). For
a detailed calculation of (2.13) and (2.14), please refer to Appendix A.

##  3 Existence of minimizers

We use the direct method in calculus of variations [14] to prove the existence
of a minimizer to the diffuse-interface LdG functional.

######  Lemma 1.

The functional $E_{\varepsilon}$ (2.9) is weak lower semi-continuous
(w.l.s.c.) on the Banach space $(\mathbf{Q},\phi)\in
H_{0}^{1}(\Omega;\mathcal{S}_{0})\times H_{0}^{1}(\Omega)$, which means that
if $\mathbf{Q}_{k}\rightharpoonup\mathbf{Q},\phi_{k}\rightharpoonup\phi_{k}$
are weakly convergent sequences in $H^{1}$,

$\liminf_{k\to\infty}E_{\varepsilon}[\mathbf{Q}_{k},\phi_{k}]\geq
E_{\varepsilon}[\mathbf{Q},\phi].$

###### Proof.

As [14, §8.2,Theorem 1], for an integral energy functional of the form

$\int_{\Omega}F(\nabla\bm{u},\bm{u},x)\mathop{}\\!\mathrm{d}x$

where $\bm{u}:\Omega\to\mathbb{R}^{m}$ and $F:\mathbb{R}^{m\times
n}\times\mathbb{R}^{m}\times\Omega\to\mathbb{R}$, if $F$ is uniformly bounded
below and convex in its first argument the functional is w.l.s.c. with respect
to $\bm{u}\in H^{1}(\Omega;\mathbb{R}^{m})$. We look at individual terms.

Examine the form (2.9) where we take
$\bm{u}=(\mathbf{Q},\phi):\Omega\to\mathcal{S}_{0}\times\mathbb{R}$. We notice
that the energy density is bounded below because all terms are positive except
for the LdG bulk energy $\bar{F}_{b}=C^{-1}F_{b}$ (2.4), which satisfies that

(3.1) $\bar{F}_{b}(\mathbf{Q})\geq
m_{b}>-\infty,\quad\forall\mathbf{Q}\in\mathcal{S}_{0}$

for some number $m_{b}\in\mathbb{R}$ whose explicit expression can be found in
[43, eq. (41)]. Moreover, the terms involving the gradients, which read

$\frac{1}{2}|\gradient\mathbf{Q}|^{2}+\omega_{p}\lambda\varepsilon|\nabla\phi|^{2}+\omega_{a}\lambda\varepsilon\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\nabla\phi\right|^{2},$

are convex with respect to $(\gradient\mathbf{Q},\nabla\phi)$ since they are
positive-defintie quadratic polynomials. Therefore, the energy density of
$E_{\varepsilon}$ (2.9) satisfies the criterion, and is w.l.s.c. over the
space $H_{0}^{1}(\Omega;\mathcal{S}_{0})\times H_{0}^{1}(\Omega)$. ∎

######  Proposition 2.

The energy functional $E_{\varepsilon}$ (2.9) has at least one global
minimizer
$(\mathbf{Q}^{*}_{\varepsilon},\phi^{*}_{\varepsilon})\in\mathscr{A}$, with
the admissible set $\mathscr{A}$ defined as in (2.11).

###### Proof.

The admissible set $\mathscr{A}$ (2.11) is obviously nonempty since there
always exists smooth $\phi\in H_{0}^{1}$ such that $\int_{\Omega}\phi=V_{0}$.

Moreover, the fact that all terms in $E_{\varepsilon}$ are bounded below
(recall (3.1)) implies the following coercivity condition.

(3.2) $\displaystyle E_{\varepsilon}[\mathbf{Q},\phi]$
$\displaystyle\geq\int_{\Omega}\left[\frac{1}{2}|\gradient\mathbf{Q}|^{2}+\omega_{p}\varepsilon|\nabla\phi|^{2}+\lambda^{2}m_{b}\right]\mathop{}\\!\mathrm{d}x$
$\displaystyle=\frac{L}{2}\|\gradient\mathbf{Q}\|_{L^{2}}^{2}+\omega_{p}\varepsilon\|\nabla\phi\|_{L^{2}}^{2}+\lambda^{2}m_{b}|\Omega|.$

The weak lower semi-continuity of $E_{\varepsilon}$ over the feasible set is
guaranteed by Lemma 1 since the admissible set $\mathscr{A}$ is a closed
subspace of $H_{0}^{1}(\Omega;\mathcal{S}_{0})\times H_{0}^{1}(\Omega)$.
Therefore, the existence of a global minimizer follows from the direct method
in calculus of variations. ∎

######  Corollary 3.

For each fixed $\phi\in H_{0}^{1}$, the w.l.s.c. and coercivity argument still
applies with respect to the variable $\mathbf{Q}$. Hence, there exists a
$\mathbf{Q}$-minimizer $\mathbf{Q}_{\varepsilon}\in
H_{0}^{1}(\Omega;\mathcal{S}_{0})$ such that

$\mathbf{Q}_{\varepsilon}=\mathbf{Q}_{\varepsilon}(\phi)=\mathop{\arg\min}_{\mathbf{Q}\in
H_{0}^{1}(\Omega;\mathcal{S}_{0})}E_{\varepsilon}[\mathbf{Q},\phi].$

##  4 Maximum principle on a fixed domain

In this section, we prove that under a fixed $\phi$, the critical points of
$E_{\varepsilon}$ with respect to $\mathbf{Q}$ (referred to as
$\mathbf{Q}$-critical points) are bounded in $L^{\infty}$ uniformly as
$\varepsilon\to 0$. The uniqueness of $\mathbf{Q}$-critical points follows
immediately. The maximum principle and uniqueness results of the LdG model on
a fixed region have been studied thoroughly [39, 34], but results for the
diffuse-interface model (2.1) involve different techniques, especially when
dealing with the anchoring penalty (2.7).

###  4.1 Uniform maximum principle

Before introducing the uniform maximum principle, we improve the regularity of
$\mathbf{Q}_{\varepsilon}$ with the $L^{p}$ theory of second-order elliptic
equations [18]. We denote by

$W^{m,p}(\Omega)=\quantity{u\in
L^{p}(\Omega):\|u\|_{W^{m,p}}^{p}=\sum_{k=0}^{m}\int_{\Omega}|\nabla^{k}u|^{p}<\infty}$

the general Sobolev spaces.

######  Lemma 4.

If $\phi\in W^{1,\infty}(\Omega)$, then a $\mathbf{Q}$-critical point
$\mathbf{Q}_{\varepsilon}\in H_{0}^{1}(\Omega;\mathcal{S}_{0})$ of
$E_{\varepsilon}[\cdot,\phi]$ (2.9) is in $W^{2,p}$ for all $p<\infty$ and
$C^{1,\alpha}$ for all $\alpha<1$. Since $\mathbf{Q}_{\varepsilon}$ is a weak
solution, continuity means that $\mathbf{Q}_{\varepsilon}$ is equal to a
$C^{1,\alpha}$ function a.e..

###### Proof.

By definition, the critical point $\mathbf{Q}_{\varepsilon}$ satisfies the
equation (2.13). Examining its right-hand side, we find that it is a cubic
polynomial of $\mathbf{Q}$ (because of the bulk energy derivative), and its
coefficients, excluding constants, are related to $\phi,\nabla\phi$, which are
uniformly bounded by assumption.

If $\mathbf{Q}_{\varepsilon}\in L^{p}(\Omega),$ then the RHS (right-hand side)
of (2.13) as a whole is in $L^{p/3}(\Omega)$. Therefore, by the $L^{p}$ theory
of elliptic equations [18], $\mathbf{Q}_{\varepsilon}\in W^{2,\frac{p}{3}}$ if
$p>3$. This enables us to improve the regularity of $\mathbf{Q}_{\varepsilon}$
through a chain of Sobolev embeddings.

We have the a priori estimate $\mathbf{Q}_{\varepsilon}\in H_{0}^{1},$ which
indicates $\mathbf{Q}_{\varepsilon}\in L^{6}$ by the Sobolev embedding theorem
in 3D. Carrying out the chain of embeddings, we eventually reach
$\mathbf{Q}_{\varepsilon}\in W^{2,p}$ for all $p<\infty$, and with one extra
step of embedding we get $\mathbf{Q}_{\varepsilon}\in C^{1,\alpha}$ for all
$\alpha<1$. ∎

Under stronger smoothness assumptions, we prove that $\mathbf{Q}$-critical
points of the nondimensional energy (2.9) is uniformly bounded. The bound is
uniform in the sense that it is independent of the capillary width
$\varepsilon$.

######  Theorem 5.

Let $\phi_{\varepsilon}\in W^{1,\infty}(\Omega)$ be a phase field, with the
following extra assumptions.

  1. (i)

$\phi_{\varepsilon}\in[0,1]$ for all $x$.

  2. (ii)

Over the diffuse interface
$\Gamma_{\varepsilon}=\overline{\\{x\in\Omega:0<\phi(x)<1\\}},$ the diffuse
unit normal vector $\mathbf{n}_{\varepsilon}$, defined by

(4.1)
$\mathbf{n}_{\varepsilon}=\begin{cases}\dfrac{\nabla\phi_{\varepsilon}}{|\nabla\phi_{\varepsilon}|},&\nabla\phi_{\varepsilon}\neq
0,\\\ \bm{0},&\nabla\phi_{\varepsilon}=0,\end{cases}$

is part of a $C^{2}$-smooth vector field on $\Omega$ such that

$\begin{cases}|\mathbf{n}_{\varepsilon}|\leq 1,\\\
\|\mathbf{n}_{\varepsilon}\|_{C^{2}(\Omega)}\leq M_{N}<\infty\end{cases}$

for all $\varepsilon$. In other words, the normal vector in the diffuse
interface $\Gamma$ is part of a uniform $C^{2}$ vector field.

Suppose $\mathbf{Q}_{\varepsilon}$ is a $\mathbf{Q}$-critical point of
$E_{\varepsilon}[\cdot,\phi_{\varepsilon}]$ (2.9). Then with any fixed
$\lambda_{0}>0$, there exists a uniform constant $M$ determined by
$A,B,C,\omega_{v},\lambda_{0}$ and the bound $M_{N}$, but independent of
$\varepsilon$, such that

$\sup_{x\in\Omega}|\mathbf{Q}_{\varepsilon}|\leq M\lambda^{-\frac{2}{3}}.$

The bound is given in the a.e. sense.

###### Proof.

In the following proof we drop the subscript $\varepsilon$ for notational
simplicity.

Take the dot product of $\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}$ with
the Euler-Lagrange equation (2.13), where $\alpha$ is a number yet to be
determined. Recall the operator $\mathscr{P}_{0}$ as defined in (2.15).

(4.2)
$\displaystyle(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\triangle\mathbf{Q}$
$\displaystyle=\lambda^{2}\underbrace{(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\quantity[\quantity(\frac{A}{C}+\tr\mathbf{Q}^{2})\mathbf{Q}-\frac{B}{C}\mathscr{P}_{0}(\mathbf{Q}^{2})]}_{(I)}$
$\displaystyle\quad+\omega_{v}\lambda^{2}(1-\phi)^{2}\underbrace{(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\mathbf{Q}}_{(II)}$
$\displaystyle\quad+\omega_{a}\lambda\varepsilon\underbrace{(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\mathscr{P}_{0}\quantity[(\nabla\phi\otimes\nabla\phi)\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})+\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})(\nabla\phi\otimes\nabla\phi)]}_{(III)}.$

For the LHS (left-hand side) of (4.2), we compute

$\displaystyle(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\triangle\mathbf{Q}$
$\displaystyle=(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\triangle(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n})-(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\triangle(\alpha\mathbf{n}\otimes\mathbf{n})$
$\displaystyle=\frac{1}{2}\triangle(|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2})-|\bm{\nabla}(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n})|^{2}-\alpha(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\triangle(\mathbf{n}\otimes\mathbf{n}).$

We have used the elementary equality
$\mathbf{A}:\triangle\mathbf{A}=\frac{1}{2}\triangle(|\mathbf{A}|^{2})-|\gradient\mathbf{A}|^{2}$
($\mathbf{A}=\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}$). It is well-
defined in the weak sense because
$\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}\in C^{0}\cap H^{1}$ by the
regularity estimate from Lemma 4.

For the RHS of (4.2), we compute the three terms separately. The first term
equals

$\displaystyle(I)$
$\displaystyle=(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\quantity[\quantity(\frac{A}{C}+\tr\mathbf{Q}^{2})\mathbf{Q}-\frac{B}{C}\quantity(\mathbf{Q}^{2}-\frac{\tr\mathbf{Q}^{2}}{3}\mathbf{I})]$
$\displaystyle=\quantity(\frac{A}{C}+\tr\mathbf{Q}^{2})(\tr\mathbf{Q}^{2}+\alpha\mathbf{n}\cdot\mathbf{Qn})-\frac{B}{C}\quantity(\tr\mathbf{Q}^{3}+\alpha|\mathbf{Qn}|^{2}-\frac{\alpha\tr\mathbf{Q}^{2}}{3}).$

The second term equals

$(II)=\tr\mathbf{Q}^{2}+\alpha\tr(\mathbf{n}\otimes\mathbf{n}\mathbf{Q})=\tr\mathbf{Q}^{2}+\alpha\mathbf{n}\cdot\mathbf{Qn}.$

We have extensively used the equality $\tr\mathbf{AB}=\tr\mathbf{BA}$.

Suppose first that $x\in\Gamma$, so by the assumption (4.1) we have
$\nabla\phi(x)=|\nabla\phi(x)|\mathbf{n}(x)$. Then we write

$\displaystyle(III)$
$\displaystyle=|\nabla\phi|^{2}(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\mathscr{P}_{0}\quantity[(\mathbf{n}\otimes\mathbf{n})\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})+\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})(\mathbf{n}\otimes\mathbf{n})]$
$\displaystyle=2|\nabla\phi|^{2}(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\mathscr{P}_{0}\quantity[(\mathbf{n}\otimes\mathbf{n})\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})]$

using
$\mathbf{n}\otimes\mathbf{n}=\frac{1}{|\nabla\phi|^{2}}\nabla\phi\otimes\nabla\phi$
and that $\mathbf{Q}$ is symmetric. Notice that since the image of
$\mathscr{P}_{0}$ has zero trace, the product would be identical if one added
$\frac{s_{+}}{3}\mathbf{I}$ to $\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}$.
Denote by $\tilde{\mathbf{Q}}=\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I}$ (so
$\tr\tilde{\mathbf{Q}}=s_{+}$), and we get

$\displaystyle(III)$
$\displaystyle=2|\nabla\phi|^{2}(\tilde{\mathbf{Q}}+\alpha\mathbf{n}\otimes\mathbf{n}):\mathscr{P}_{0}((\mathbf{n}\otimes\mathbf{n})\tilde{\mathbf{Q}})$
$\displaystyle=2|\nabla\phi|^{2}(\tilde{\mathbf{Q}}+\alpha\mathbf{n}\otimes\mathbf{n}):\quantity[(\mathbf{n}\otimes\mathbf{n})\tilde{\mathbf{Q}}-\frac{\tr((\mathbf{n}\otimes\mathbf{n})\tilde{\mathbf{Q}})}{3}\mathbf{I}]$
$\displaystyle=2|\nabla\phi|^{2}\quantity[|\tilde{\mathbf{Q}}\mathbf{n}|^{2}+\frac{2\alpha-
s_{+}}{3}(\mathbf{n}\cdot\tilde{\mathbf{Q}}\mathbf{n})].$

Setting $\alpha=\frac{s_{+}}{2}$ eliminates the term $\frac{2\alpha-
s_{+}}{3}(\mathbf{n}\cdot\tilde{\mathbf{Q}}\mathbf{n})$, which leads us to

$(III)=2\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\nabla\phi\right|^{2},\quad
x\in\Gamma.$

If $x\notin\Gamma$, then $\phi$ is a constant in a neighbourhood of $x$, so
$\nabla\phi(x)\equiv 0$ and the third term equals zero. Therefore, it holds
for all $x\in\Omega$ that

$(III)\geq 0.$

Then, substituting the results above in (4.2), and removing all positive terms
on the RHS, we get the following inequality satisfied by the function
$|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2}$ over the entire $\Omega$,
known as a subsolution condition.

(4.3)
$-\triangle\quantity(|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2})+2T(\mathbf{Q})\leq
0.$

The function $T(\mathbf{Q})$ is defined as

(4.4) $\displaystyle T(\mathbf{Q})$
$\displaystyle=T^{(1)}(\mathbf{Q})+\lambda^{2}T^{(2)}(\mathbf{Q}),$
$\displaystyle T^{(1)}(\mathbf{Q})$
$\displaystyle=\alpha(\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}):\triangle(\mathbf{n}\otimes\mathbf{n})$
$\displaystyle T^{(2)}(\mathbf{Q})$
$\displaystyle=(\tr\mathbf{Q}^{2})^{2}-\frac{B}{C}\tr\mathbf{Q}^{3}+\alpha(\tr\mathbf{Q}^{2})(\mathbf{n}\cdot\mathbf{Qn})+\frac{3A+B\alpha}{3C}\tr\mathbf{Q}^{2}-\frac{B}{C}\alpha|\mathbf{Qn}|^{2}$
$\displaystyle\quad+\alpha\quantity(\frac{A}{C}+\omega_{v}(1-\phi)^{2})\mathbf{n}\cdot\mathbf{Qn}.$

Despite the complexity of its expression, $T(\mathbf{Q})$ is basically a
fourth order polynomial in terms of $\mathbf{Q}$, with uniformly bounded
coefficients and a strictly positive leading term
$\lambda^{2}(\tr\mathbf{Q}^{2})^{2}$. Therefore, there exists a number $R_{0}$
depending on $A,B,C,\lambda,\omega_{v}$ and the $C^{2}$ norm $M_{N}$ from
assumption (ii), such that

(4.5) $T(\mathbf{Q})\geq 0,\quad\forall|\mathbf{Q}|\geq R_{0}.$

Note that $\alpha=\frac{s_{+}}{2}$ is also determined by $A,B,C$ according to
(2.5).

Finally, we use the standard technique in the maximum principle estimate of
weak solutions. One can consult [34, App. 2] for an example. Define the cut-
off

$\psi_{k}=\max\quantity{|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2}-k,0},$

such that $\psi_{k}\in H_{0}^{1}(\Omega)$ for $k>\alpha^{2}$ because
$\mathbf{Q}|_{\partial\Omega}\equiv 0$ and that
$|\alpha\mathbf{n}\otimes\mathbf{n}|^{2}\leq\alpha^{2}$. By the chain rule, we
also have

(4.6)
$\nabla\psi_{k}=\begin{cases}\nabla(|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2}),&|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2}>k,\\\
0,&|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2}\leq k,\end{cases}$

Multiplying both sides of (4.3) with $\psi_{k}$, integrating by parts and
eliminating positive terms, we get an inequality

(4.7)
$\int_{\Omega}\nabla(|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2})\cdot\nabla\psi_{k}\mathop{}\\!\mathrm{d}x+2\int_{\Omega}\psi_{k}T(\mathbf{Q})\mathop{}\\!\mathrm{d}x\leq
0,\quad k>\alpha^{2}.$

By (4.6), the first term of the LHS of (4.7) is simply
$\int_{\Omega}|\nabla\psi_{k}|^{2}.$ The second term of (4.7) is supported on
the region where $\psi_{k}\neq 0$, i.e.
$|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|^{2}>k$. By the triangle
inequality of the Frobenius norm,

$|\mathbf{Q}|\geq|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|-|\alpha\mathbf{n}\otimes\mathbf{n}|=|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|-\alpha,$

so choosing $k>(R_{0}+\alpha)^{2}$ is sufficient to ensure
$|\mathbf{Q}|>R_{0}$ whenever $\psi_{k}$ is nonzero. Thus, according to (4.5),
$T(\mathbf{Q})$ is positive on the support of $\psi_{k}$, and we get

$\int_{\Omega}|\nabla\psi_{k}|^{2}\mathop{}\\!\mathrm{d}x\leq
0\Rightarrow\psi_{k}\equiv 0,$

which means that $|\mathbf{Q}+\alpha\mathbf{n}\otimes\mathbf{n}|$ is uniformly
bounded by $k$ for all $k>(R_{0}+\alpha)^{2}$. Therefore, using the triangle
inequality on $\mathbf{Q}$ and $\alpha\mathbf{n}\otimes\mathbf{n}$ again, we
get

(4.8) $\sup_{x\in\Omega}|\mathbf{Q}|\leq R_{0}+\alpha,$

which is the uniform bound we desired.

Evidently, the bound (4.8) depends on the choice of $R_{0}$, which ensures the
positivity of $T(\mathbf{Q})$ for all $|\mathbf{Q}|>R_{0}$. Hence, to complete
the proof, we also need to study the dependence of $R_{0}$ (4.5) on $\lambda$.
Closely observing the expressions of $T^{(1)}$ and $T^{(2)}$ (4.4), we find
that their coefficients are controlled by a uniform bound dependent on
$A,B,C,\omega_{v},M_{N}$ but not on $\lambda$. We then get

$\displaystyle T^{(1)}(\mathbf{Q})\geq-C_{1}(|\mathbf{Q}|+1),$ $\displaystyle
T^{(2)}(\mathbf{Q})\geq\frac{1}{2}|\mathbf{Q}|^{4}-C_{2},$

with constants $C_{1},C_{2}>0$ depending on $A,B,C,\omega_{v},M_{N}$, since
$T^{(1)}$ is affine in $\mathbf{Q}$ and $T^{(2)}$ is degree-4 in $\mathbf{Q}$
with leading term $(\tr\mathbf{Q}^{2})^{2}$. In order that $T(\mathbf{Q})\geq
0$, it suffices to force the lower bound to be positive, i.e.

$T^{(1)}+\lambda^{2}T^{(2)}\geq\lambda^{2}\quantity(\frac{1}{2}|\mathbf{Q}|^{4}-C_{2})-C_{1}(|\mathbf{Q}|+1)\geq
0.$

Multiplying both sides by $\lambda^{\frac{2}{3}}$, it reduces to checking that

$\frac{1}{2}|\lambda^{\frac{2}{3}}\mathbf{Q}|^{4}-C_{1}|\lambda^{\frac{2}{3}}\mathbf{Q}|\geq
C_{1}\lambda^{\frac{2}{3}}+C_{2}\lambda^{\frac{8}{3}}.$

Suppose that $\lambda_{0}>0$ is given, so the RHS is bounded for all
$\lambda<\lambda_{0}$, whose bound depends on $C_{1},C_{2},\lambda_{0}$.
Notice also that the LHS is a degree-4 polynomial of
$|\lambda^{\frac{2}{3}}\mathbf{Q}|$ as a whole. By the elementary property of
polynomial functions, there exists $C_{3}>0$ depending on
$C_{1},C_{2},\lambda_{0}$ such that the inequality holds whenever
$|\lambda^{\frac{2}{3}}\mathbf{Q}|>C_{3}$ i.e.
$|\mathbf{Q}|>C_{3}\lambda^{-\frac{2}{3}}$. Therefore, we take
$R_{0}=C_{3}\lambda^{-\frac{2}{3}}$ in the uniform bound (4.8) and get

$\sup_{x\in\Omega}|\mathbf{Q}|\leq C_{3}\lambda^{-\frac{2}{3}}+\alpha\leq
M\lambda^{-\frac{2}{3}},\quad\forall\lambda<\lambda_{0}$

with $M=C_{3}+\alpha\lambda_{0}^{\frac{2}{3}}$ dependent on
$A,B,C,L,\omega_{v},\lambda_{0},M_{N}$ as desired. ∎

###  4.2 Uniqueness on small domains

An immediate result of the maximum principle is the uniqueness of the critical
points when $\lambda$ is small.

######  Proposition 6.

Under the assumptions of Theorem 5, there exists a $\lambda_{1}>0$ depending
on $\Omega,A,B,C,L,\omega_{v}$ and the bound $M_{N}$, such that for all
$\lambda<\lambda_{1}$, the $\mathbf{Q}$-critical point of
$E_{\varepsilon}[\cdot,\phi]$ (2.9) is unique.

###### Proof.

Following the approach of Lemma 8.2 from [34, App. 2], we only need to prove
the strict convexity of the energy functional $E_{\varepsilon}$ on the closed
convex set

$X_{\lambda}=\left\\{\mathbf{Q}\in
H_{0}^{1}(\Omega;\mathcal{S}_{0}):\sup_{x\in\Omega}|\mathbf{Q}|\leq
M\lambda^{-\frac{2}{3}}\right\\},$

where $M=M(A,B,C,L,\omega_{v},\lambda_{0},M_{N})$ is the constant stated in
Theorem 5, with $\lambda_{0}=1$. By definition, $X_{\lambda}$ contains all
critical points of the energy. We highlight the dependence on $\lambda$ since
the bound changes as $\lambda\to 0$. It suffices to check that

$\frac{E_{\varepsilon}[\mathbf{Q}+\mathbf{h},\phi]+E_{\varepsilon}[\mathbf{Q}-\mathbf{h},\phi]}{2}\geq
E_{\varepsilon}\quantity[\mathbf{Q},\phi]$

if $\mathbf{Q}\pm\mathbf{h}\in X_{\lambda}$ and that the strict inequality
holds when $\mathbf{h}\neq 0$. Since the void energy and the anchoring energy
are all positive-definite quadratic polynomials with respect to $\mathbf{Q}$,
they are both convex, and we only need to check the LdG energy

$E^{\rm
LdG}[\mathbf{Q}]=\int_{\Omega}\quantity[\frac{1}{2}|\nabla\mathbf{Q}|^{2}+\lambda^{2}\bar{F}_{b}(\mathbf{Q})]\mathop{}\\!\mathrm{d}x.$

By the convexity of $X_{\lambda}$, $\mathbf{Q}\pm\mathbf{h}\in X_{\lambda}$
implies that $\mathbf{Q}\in X_{\lambda}$. Then, strong convexity of $E^{\rm
LdG}$ is equivalent to the positivity of the second-order central difference

$\Delta(\mathbf{Q},\mathbf{h})=E^{\rm LdG}[\mathbf{Q}+\mathbf{h}]+E^{\rm
LdG}[\mathbf{Q}-\mathbf{h}]-2E^{\rm LdG}[\mathbf{Q}]$

whenever $\mathbf{h}\neq\mathbf{0}$. From the expression of $E^{\rm LdG}$ we
compute

$\Delta(\mathbf{Q},\mathbf{h})=\int_{\Omega}|\gradient\mathbf{h}|^{2}+\lambda^{2}\int_{\Omega}\quantity[\frac{A}{C}\tr\mathbf{h}^{2}-\frac{2B}{C}\tr(\mathbf{Qh}^{2})+2(\tr\mathbf{Qh})^{2}+\tr\mathbf{Q}^{2}\tr\mathbf{h}^{2}+\frac{1}{2}(\tr\mathbf{h}^{2})^{2}].$

Since $|\mathbf{Q}|\leq M\lambda^{-\frac{2}{3}}$, we can estimate a lower
bound of $\Delta(\mathbf{Q},\mathbf{h})$ by removing all positive terms and
amplifying $\mathbf{Q}$ to its largest norm. We get

$\displaystyle\Delta(\mathbf{Q},\mathbf{h})$
$\displaystyle\geq\int_{\Omega}|\gradient\mathbf{h}|^{2}\mathop{}\\!\mathrm{d}x-\lambda^{2}\int_{\Omega}\quantity[\frac{|A|}{C}\tr\mathbf{h}^{2}+\frac{2B}{C}\cdot
M\lambda^{-\frac{2}{3}}\tr\mathbf{h}^{2}]\mathop{}\\!\mathrm{d}x$
$\displaystyle\geq\int_{\Omega}|\gradient\mathbf{h}|^{2}\mathop{}\\!\mathrm{d}x-M^{\prime}\lambda^{\frac{4}{3}}\int_{\Omega}|\mathbf{h}|^{2}\mathop{}\\!\mathrm{d}x,$

where the constant $M^{\prime}$ is scaled from $M$ by a constant related to
$\frac{|A|}{C},\frac{B}{C}$. We note that the factor $\lambda^{2}$ manages to
absorb the unbounded $\lambda^{-\frac{2}{3}}$ from the uniform maximum value
introduced in Theorem 5.

Poincaré’s inequality [18] implies that

$\int_{\Omega}|\gradient\mathbf{h}|^{2}\mathop{}\\!\mathrm{d}x\geq\frac{1}{C_{4}(\Omega)}\int_{\Omega}|\mathbf{h}|^{2}\mathop{}\\!\mathrm{d}x,$

so there exists a
$\lambda_{1}=\min\\{\lambda_{0},(C_{4}M^{\prime})^{-\frac{3}{4}}\\}$ such that
when $\lambda<\lambda_{1}$, the quantity $\Delta(\mathbf{Q},\mathbf{h})$ is
always strictly positive for all $\mathbf{h}\neq\mathbf{0}$. Therefore,
$\bar{E}^{LdG}$ is strictly convex on $X_{\lambda}$, and the uniqueness of
$\mathbf{Q}$-critical points follows naturally. ∎

##  5 Sharp-interface limit

The limiting behaviour of $E_{\varepsilon}[\mathbf{Q},\phi]$ (2.9) when the
capillary width $\varepsilon\to 0$ is known as the sharp-interface limit. A
famous result of the sharp-interface limit is about van der Waals-Cahn-
Hilliard functional (as used in our mixing energy (2.6)), which states that it
converges to the interface perimeter in the sense of $\Gamma$-limit [42, 55].
We are interested in the sharp-interface limit $\varepsilon\to 0$ because it
guarantees that solving the continuous diffuse-interface functional (e.g. the
van der Waals-Cahn-Hilliard functional) approximates difficult-to-compute
sharp-interface problems (e.g. interface perimeter).

For our diffuse-interface LdG energy $E_{\varepsilon}$ (2.9), we also assert
such a sharp-interface limit.

######  Theorem 7.

Let $(\mathbf{Q}_{\varepsilon}^{*},\phi_{\varepsilon}^{*})\in\mathscr{A}$ be a
minimizer of $E_{\varepsilon}$ (2.9) in its admissible space (2.11). Assume
further that as $\varepsilon\to 0$,

  1. (i)

$\\{\phi_{\varepsilon}^{*}\\}$ converges to the indicator function
$I_{D_{0}}=\begin{cases}1,&x\in D_{0},\\\ 0&x\notin D_{0}\end{cases}$ in
$L^{1}$, where $D_{0}\subset\subset\Omega$ (abbreviation for
$D_{0}\subset\Omega$ and $\partial D_{0}\subset\Omega$, used throughout the
rest of the text) such that $|D_{0}|=V_{0}$;

  2. (ii)

$\\{\mathbf{Q}_{\varepsilon}^{*}\\}$ converges uniformly to a continuous
function $\mathbf{Q}_{0}\in C(\Omega;\mathcal{S}_{0})$.

Then, $(\mathbf{Q}_{0},D_{0})$ minimizes the energy functional

(5.1) $\displaystyle E_{0}[\mathbf{Q},D]$
$\displaystyle=\int_{\Omega}\quantity[\frac{L}{2}|\gradient\mathbf{Q}|^{2}+\lambda^{2}\bar{F}_{b}(\mathbf{Q})]\mathop{}\\!\mathrm{d}x+\omega_{v}\lambda^{2}\int_{D^{c}}\frac{1}{2}|\mathbf{Q}|^{2}\mathop{}\\!\mathrm{d}x$
$\displaystyle\quad+\frac{\omega_{p}\lambda}{3}\int_{\partial
D}\quantity[1+\frac{\omega_{a}}{\omega_{p}}\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\bm{\nu}\right|^{2}]^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x)$

among all matrix functions $\mathbf{Q}\in C^{0}\cap H^{1}$ and regions
$D\subset\subset\Omega$ with smooth boundary satisfying $|D|=V_{0}$, where
$\bm{\nu}$ represents the unit normal vector of $\partial D$.

######  Remark 1.

The sharp-interface limit (5.1) also comprises the LdG free energy and penalty
factors on the void, the perimeter and the weak anchoring on $\partial D$. The
first integral is identical to $E^{\rm LdG}[\mathbf{Q}]$, and the second is
identical to $E^{\rm void}[\mathbf{Q},I_{D}]$ (2.8) where $I_{D}$ is the
indicator function of $D$. With an elementary inequality
$\frac{1+|x|}{\sqrt{2}}\leq\sqrt{1+|x|^{2}}\leq 1+|x|,$ we find that

$1\leq\frac{\displaystyle\mathcal{P}(D)+\sqrt{\frac{\omega_{a}}{\omega_{p}}}\int_{\partial
D}\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\bm{\nu}\right|}{\displaystyle\int_{\partial
D}\quantity[1+\frac{\omega_{a}}{\omega_{p}}\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\bm{\nu}\right|^{2}]^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x)}\leq\sqrt{2}.$

That is, the boundary integral is equivalent to the sum of the perimeter and
an $L^{1}$ anchoring penalty.

###  5.1 Proof by $\Gamma$-convergence

Theorem 7 is proven in the framework of $\Gamma$-convergence, which is
composed of an upper bound estimate and a lower bound estimate [10].

We start by rewriting the diffuse-interface functional (2.9) as

(5.2) $E_{\varepsilon}[\mathbf{Q},\phi]=E^{\rm
LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q},\phi]+\omega_{p}\lambda\int_{\Omega}\quantity[\varepsilon
a_{\mathbf{Q}}(\nabla\phi)+\varepsilon^{-1}W(\phi)]\mathop{}\\!\mathrm{d}x$

with $W(s)=s^{2}(1-s)^{2}$ denoting the double-well function, and the sharp-
interface functional (5.1) as

(5.3) $E_{0}[\mathbf{Q},D]=E^{\rm LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q},I_{D}]+\frac{\omega_{p}\lambda}{3}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x).$

Note that the expressions (5.2) and (5.3) share the LdG and void energy, as
well as a quadratic form $a_{\mathbf{Q}}(\cdot):\mathbb{R}^{3}\to\mathbb{R}$
defined by

(5.4)
$a_{\mathbf{Q}}(\xi)=|\xi|^{2}+\frac{\omega_{a}}{\omega_{p}}\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\xi\right|^{2},\quad\xi\in\mathbb{R}^{3},$

which strongly implies the connection between $E_{\varepsilon}$ and $E_{0}$.

Then, we state the $\Gamma$-convergence of $E_{\varepsilon}$ to $E_{0}$.

######  Proposition 8.

Work under the same assumptions as Theorem 7. Then, the following properties
hold for all $\mathbf{Q}\in C^{0}\cap H^{1}$ and smooth region
$D\subset\subset\Omega$ satisfying $|D|=V_{0}$.

  1. (a)

For all sequences $\phi_{\varepsilon}\to I_{D}$ ($I_{D}$ is the indicator
function of $D$) in $L^{1}$, $\mathbf{Q}_{\varepsilon}\to\mathbf{Q}$
uniformly,

(5.5) $\liminf_{\varepsilon\to
0}E_{\varepsilon}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}]\geq
E_{0}[\mathbf{Q},D];$

  2. (b)

There exists a sequence $\phi_{\varepsilon}\to I_{D}$ in $L^{1}$ satisfying
the volume constraint (2.10), $\mathbf{Q}_{\varepsilon}\to\mathbf{Q}$
uniformly, such that

(5.6) $\limsup_{\varepsilon\to
0}E_{\varepsilon}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}]\leq
E_{0}[\mathbf{Q},D].$

Proposition 8 relies on the following claims. We work under the same
assumptions as Proposition 8 throughout the claims.

  * •

Claim (I): For all sequences $\phi_{\varepsilon}\to I_{D}$ in $L^{1}$ and
$\mathbf{Q}_{\varepsilon}\to\mathbf{Q}$ uniformly,

(5.7) $\liminf_{\varepsilon\to 0}\quantity(E^{\rm
LdG}[\mathbf{Q}_{\varepsilon}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}])\geq E^{\rm
LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm void}[\mathbf{Q},I_{D}]$

  * •

Claim (II): For all $\phi_{\varepsilon}\to I_{D}$ a.e. with
$0\leq\phi_{\varepsilon}\leq 1$, it holds that

(5.8) $\lim_{\varepsilon\to 0}\quantity(E^{\rm
LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q},\phi_{\varepsilon}])=E^{\rm
LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm void}[\mathbf{Q},I_{D}]$

  * •

Claim (III): For all sequences $\phi_{\varepsilon}\to I_{D}$ in $L^{1}$,
$\mathbf{Q}_{\varepsilon}\to\mathbf{Q}$ uniformly,

(5.9) $\liminf_{\varepsilon\to 0}\int_{\Omega}\quantity[\varepsilon
a_{\mathbf{Q}_{\varepsilon}}(\nabla\phi_{\varepsilon})+\varepsilon^{-1}W(\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x\geq\frac{1}{3}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x).$

  * •

Claim (IV): There exists a sequence $\phi_{\varepsilon}\to I_{D}$ a.e. with
$0\leq\phi_{\varepsilon}\leq 1$ and the volume constraint (2.10), such that

(5.10) $\limsup_{\varepsilon\to 0}\int_{\Omega}\quantity[\varepsilon
a_{\mathbf{Q}}(\nabla\phi_{\varepsilon})+\varepsilon^{-1}W(\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x\leq\frac{1}{3}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x).$

Obviously, Claims (I) and (II) concern the convergence of the first two terms
of (5.2) and (5.3), while Claims (III) and (IV) are about the final term. For
technical reasons, we note that the sequences in Claim (II) and (IV) have
$\mathbf{Q}$ fixed, and that the sequences in Claim (IV) satisfy Claim (II).

###  5.2 Detailed proof

We start by proving Claims (I)–(IV). Claims (I) and (II) can be obtained
through basic calculus.

###### Proof of Claim (I).

Assume w.l.o.g. that

$\liminf_{\varepsilon\to
0}\int_{\Omega}|\gradient\mathbf{Q}_{\varepsilon}|^{2}\mathop{}\\!\mathrm{d}x<\infty,$

for otherwise the conclusion is trivial. Due to the weak compactness of
$H^{1}$, we extract a subsequence of $\\{\mathbf{Q}_{\varepsilon}\\}$, denoted
by $\\{\mathbf{Q}_{\varepsilon_{k}}\\}$ such that

$\lim_{k\to\infty}\int_{\Omega}|\gradient\mathbf{Q}_{\varepsilon_{k}}|^{2}\mathop{}\\!\mathrm{d}x=\liminf_{\varepsilon\to
0}\int_{\Omega}|\gradient\mathbf{Q}_{\varepsilon}|^{2}\mathop{}\\!\mathrm{d}x\text{
and }\mathbf{Q}_{\varepsilon_{k}}\rightharpoonup\mathbf{Q}\text{ in }H^{1}.$

Hence, we have by the w.l.s.c. property [14, sect. 8.2, Theorem I] that

$\liminf_{\varepsilon\to
0}\int_{\Omega}|\gradient\mathbf{Q}_{\varepsilon}|^{2}\mathop{}\\!\mathrm{d}x=\lim_{k\to\infty}\int_{\Omega}|\gradient\mathbf{Q}_{\varepsilon_{k}}|^{2}\mathop{}\\!\mathrm{d}x\geq\int_{\Omega}|\gradient\mathbf{Q}|^{2}\mathop{}\\!\mathrm{d}x.$

The uniform convergence $\mathbf{Q}_{\varepsilon}\to\mathbf{Q}$ implies that
$\bar{F}_{b}(\mathbf{Q}_{\varepsilon})\to\bar{F}_{b}(\mathbf{Q})$ uniformly as
well, since the bulk energy $\bar{F}_{b}$ (2.4) is a continuous polynomial
function. Thus,

$\lim_{\varepsilon\to
0}\int_{\Omega}\bar{F}_{b}(\mathbf{Q}_{\varepsilon})\mathop{}\\!\mathrm{d}x=\int_{\Omega}\bar{F}_{b}(\mathbf{Q})\mathop{}\\!\mathrm{d}x.$

For the void penalty $E^{\rm void}$ (2.8), we once again extract a subsequence
$\\{\phi_{\varepsilon_{k}}\\}\to I_{D}$ (w.l.o.g. the same subsequence as
before) such that

$\lim_{k\to\infty}E^{\rm
void}[\mathbf{Q}_{\varepsilon_{k}},\phi_{\varepsilon_{k}}]=\liminf_{\varepsilon\to
0}E^{\rm void}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}]\text{ and
}\phi_{\varepsilon_{k}}\to I_{D},\text{a.e.}.$

Recall the expression (2.8) of $E^{\rm void}$. It follows from Fatou’s lemma
[47] that

$\displaystyle\liminf_{\varepsilon\to 0}E^{\rm
void}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}]$
$\displaystyle=\lim_{k\to\infty}\int_{\Omega}\frac{1}{2}(1-\phi_{\varepsilon_{k}})^{2}|\mathbf{Q}_{\varepsilon_{k}}|^{2}\mathop{}\\!\mathrm{d}x$
$\displaystyle\geq\int_{\Omega}\liminf_{k\to\infty}\frac{1}{2}(1-\phi_{\varepsilon_{k}})^{2}|\mathbf{Q}_{\varepsilon_{k}}|^{2}\mathop{}\\!\mathrm{d}x=E^{\rm
void}[\mathbf{Q},I_{D}].$

Summing up the formulae above gives (5.7). ∎

###### Proof of Claim (II).

The LdG energy is the same on both sides of (5.8), so we only have to check
the void energy. Since $\phi_{\varepsilon}$ and $\mathbf{Q}\in C^{0}$ are
uniformly bounded and since $\phi_{\varepsilon}\to I_{D}$ a.e., it follows
from the dominated convergence theorem that

$\lim_{\varepsilon\to
0}\int_{\Omega}\frac{1}{2}(1-\phi_{\varepsilon})^{2}|\mathbf{Q}|^{2}\mathop{}\\!\mathrm{d}x=\int_{\Omega}\frac{1}{2}(1-I_{D})^{2}|\mathbf{Q}|^{2}\mathop{}\\!\mathrm{d}x,$

which leads to (5.8). ∎

The proof of Claims (III) and (IV) involves discussion on a generalized form
of the van der Waals-Cahn-Hilliard functional (2.6). We briefly digress from
the main thread and present a general result, which are of critical importance
to the proof. The notations used in the following two Lemmas are independent
of the rest of the paper.

######  Lemma 9.

Work on a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{n}$ in general
dimensions $n$. Define the generalized van der Waals-Cahn-Hilliard functional
$L_{\varepsilon}^{a}:H^{1}(\Omega)\to\mathbb{R}\cup\\{+\infty\\}$ as

(5.11)
$L_{\varepsilon}^{a}(\phi)=\int_{\Omega}\quantity[\varepsilon^{-1}W(\phi)+\varepsilon
a(x,\nabla\phi)]\mathop{}\\!\mathrm{d}x,$

with the following assumptions.

  1. (i)

$a(x,p)=a^{ij}(x)p_{i}p_{j}$ with Einstein’s summation convention applied. The
symmetric-matrix-valued function
$\mathbf{A}=((a^{ij})):\overline{\Omega}\to\mathbb{R}^{n\times n}$ is
continuous, and satisfies the uniform elliptic condition

(5.12)
$\lambda\mathbf{I}\preceq\mathbf{A}(x)\preceq\Lambda\mathbf{I},\quad\forall
x.$

The notation $\mathbf{A}\preceq\mathbf{B}$ means that $\mathbf{B}-\mathbf{A}$
is positive definite, where $\mathbf{A},\mathbf{B}$ are symmetric matrices.

  2. (ii)

$W(s)$ is a nonnegative continuous double-well potential with strict minima at
$s=\alpha,\beta$ ($\alpha<\beta$), and $W(\alpha)=W(\beta)=0$.

  3. (iii)

$D\subset\subset\Omega$ is a region with $C^{2}$ boundary.

Let $\phi_{0}=\alpha I_{D}+\beta I_{D^{c}}$ be the separation of the phases
$\alpha$ and $\beta$. Also, denote the perimeter of $D$ by
$\mathcal{P}(D)=\mathcal{H}_{n-1}(\partial D)$ (the $(n-1)$-dimensional
Hausdorff measure [15]) and introduce the constant

$c_{0}=\int_{\alpha}^{\beta}\sqrt{W(s)}\mathop{}\\!\mathrm{d}s.$

Then, for any sequence $\phi_{\varepsilon}\to\phi_{0}$ in $L^{1}(\Omega)$, it
holds that

(5.13) $\liminf_{\varepsilon\to 0}L_{\varepsilon}^{a}(\phi_{\varepsilon})\geq
2c_{0}\int_{\partial
D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x),$

where $\bm{\nu}:\partial D\to\mathbb{S}^{n-1}$ is the unit outward normal
vector.

######  Lemma 10.

Under the same assumptions as Lemma 9, suppose further that $\phi_{0}=\alpha
I_{D}+\beta I_{D^{c}}$ satisfies the volume constraint

(5.14) $\int_{\Omega}\phi\mathop{}\\!\mathrm{d}x=m$

where $\alpha|\Omega|<m<\beta|\Omega|$, or equivalently
$|D|=\frac{\beta|\Omega|-m}{\beta-\alpha}.$ Then, there exists a sequence
$\phi_{\varepsilon}\in W^{1,\infty}(\Omega)$ satisfying the constraint (5.14),
$\alpha\leq\phi_{\varepsilon}\leq\beta$, $\phi_{\varepsilon}\to\phi_{0}$ a.e.,
such that

(5.15) $\limsup_{\varepsilon\to 0}L_{\varepsilon}^{a}(\phi_{\varepsilon})\leq
2c_{0}\int_{\partial
D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x).$

######  Remark 2.

Remarks on the results of Lemmas 9 and 10:

  * •

These lemmas are stated in the exact same fashion as [42, Prop. 1, 2], where
the same properties are derived for the degenerate case
$a^{ij}(x,p)\equiv|p|^{2}$ and (5.11) is reduced to the classical van der
Waals-Cahn-Hilliard energy.

  * •

The theory of Lemmas 9 and 10 can apply to our model with
$\alpha=1,\beta=0,W(s)=s^{2}(1-s)^{2}$, and $\phi_{0}=I_{D}$. The reason is
that when $\alpha>\beta$, one can set $\tilde{W}(s)=W(-s)$ and apply the
theory to $\\{-\phi_{\varepsilon}\\}$ with $\tilde{W}$ substituted for $W$.

For the sake of continuity, the proof of Lemmas 9 and 10 will be left for
Appendices B and C, respectively. Assuming their correctness for the time
being, we can proceed to prove Claims (III) and (IV).

###### Proof of Claim (III).

Evidently, the quadratic form $a_{\mathbf{Q}}(x)$ (5.4) satisfies the uniform
elliptic condition (5.12) as long as $\mathbf{Q}:\Omega\to\mathcal{S}_{0}$ is
continuous. Since $\mathbf{Q}_{\varepsilon}\to\mathbf{Q}$ uniformly,
$(\mathbf{Q}_{\varepsilon}+\frac{s_{+}}{3}\mathbf{I})^{2}\to(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})^{2}$
uniformly as well. Therefore, for any $\delta>0$, there exists
$\varepsilon_{0}$ such that when $\varepsilon<\varepsilon_{0}$,

$a_{\mathbf{Q}_{\varepsilon}}(p)\geq
a_{\mathbf{Q}}(p)-\delta|p|^{2},\quad\forall x\in\Omega,p\in\mathbb{R}^{3}.$

As $a_{\mathbf{Q}}(p)\geq|p|^{2}$ by definition (see (5.4)), the quadratic
form $a_{\mathbf{Q}}(p)-\delta|p|^{2}$ is also uniformly elliptic when
$\delta$ is sufficiently small (e.g. $\delta<\frac{1}{2}$). Therefore,
applying Lemma 9 to the sequence $\\{\phi_{\varepsilon}\\}$ with quadratic
form $a(x,p)=a_{\mathbf{Q}(x)}(p)-\delta|p|^{2}$, phases $\alpha=1,\beta=0$
and region $D$ shows that

LHS of (5.9) $\displaystyle\geq\liminf_{\varepsilon\to
0}\int_{\Omega}\quantity[\varepsilon(a_{\mathbf{Q}}(\nabla\phi_{\varepsilon})-\delta|\nabla\phi_{\varepsilon}|^{2})+\varepsilon^{-1}W(\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x$
$\displaystyle\geq 2c_{0}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})-\delta|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S.$

Since $W(s)=s^{2}(1-s)^{2}$, the constant
$c_{0}=\int_{0}^{1}\sqrt{W(s)}\mathop{}\\!\mathrm{d}s=\int_{0}^{1}s(1-s)\mathop{}\\!\mathrm{d}s=\frac{1}{6}$,
so the coefficient is $2c_{0}=\frac{1}{3}$. Taking the limit $\delta\to 0$ the
boundary integral, we immediately obtain

$\text{LHS of \eqref{gamlim-3-4-inf}}\geq\frac{1}{3}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S=\text{RHS of
\eqref{gamlim-3-4-inf}}$

as desired. ∎

###### Proof of Claim (IV).

Let $\phi_{\varepsilon}$ be the sequence guaranteed by Lemma 10, with
$a(x,p)=a_{\mathbf{Q}(x)}(p)$ and $\alpha=1,\beta=0$, which satisfies
$\phi_{\varepsilon}\to I_{D}$ a.e., $0\leq\phi\leq 1$ and that the volume
constraint (2.10) holds. Then, the conclusion of the lemma shows that

LHS of (5.10) $\displaystyle=\limsup_{\varepsilon\to
0}\int_{\Omega}\quantity[\varepsilon
a_{\mathbf{Q}}(\nabla\phi_{\varepsilon})+\varepsilon^{-1}W(\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq 2c_{0}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S=\text{RHS of
\eqref{gamlim-3-4-sup}}$

as desired, since $c_{0}=\frac{1}{6}$. ∎

Having established the claims, we proceed to prove Proposition 8.

###### Proof of Proposition 8.

i. Lower bound proof. For the lower bound (5.5), we add up (5.7) and (5.9) for
all sequences $\phi_{\varepsilon}\to I_{D}$ in $L^{1}$ and
$\mathbf{Q}_{\varepsilon}\to\mathbf{Q}$ uniformly. Using Claims (I) and (III),
we conduct the following computation to obtain (5.5).

$\displaystyle\liminf_{\varepsilon\to
0}E_{\varepsilon}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}]$
$\displaystyle\geq\liminf_{\varepsilon\to 0}\quantity(E^{\rm
LdG}[\mathbf{Q}_{\varepsilon}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}])$
$\displaystyle\quad+\omega_{p}\lambda\liminf_{\varepsilon\to
0}\int_{\Omega}\quantity[\varepsilon
a_{\mathbf{Q}_{\varepsilon}}(\nabla\phi_{\varepsilon})+\varepsilon^{-1}W(\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x$
$\displaystyle\geq E^{\rm LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q},I_{D}]+\frac{\omega_{p}\lambda}{3}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x)$
$\displaystyle=E_{0}[\mathbf{Q},D].$

ii. Upper bound proof. Take the sequence guaranteed by the Claim (IV), where
$\phi_{\varepsilon}\to I_{D}$ a.e. with $0\leq\phi_{\varepsilon}\leq 1$ (so
$\\{\phi_{\varepsilon}\\}$ also satisfies Claim (II)) with the volume
constraint (2.10). We also let $\mathbf{Q}_{\varepsilon}\equiv\mathbf{Q}$ be
constant. Therefore, (5.8) and (5.10) both hold for the sequence
$\\{(\phi_{\varepsilon},\mathbf{Q})\\}$. Adding them together leads to (5.6).

$\displaystyle\limsup_{\varepsilon\to
0}E_{\varepsilon}[\mathbf{Q},\phi_{\varepsilon}]$
$\displaystyle=\lim_{\varepsilon\to 0}\quantity(E^{\rm
LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q},\phi_{\varepsilon}])$
$\displaystyle\quad+\omega_{p}\lambda\limsup_{\varepsilon\to
0}\int_{\Omega}\quantity[\varepsilon
a_{\mathbf{Q}}(\nabla\phi_{\varepsilon})+\varepsilon^{-1}W(\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq E^{\rm LdG}[\mathbf{Q}]+\omega_{v}\lambda^{2}E^{\rm
void}[\mathbf{Q},I_{D}]+\frac{\omega_{p}\lambda}{3}\int_{\partial
D}|a_{\mathbf{Q}}(\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x)$
$\displaystyle=E_{0}[\mathbf{Q},D].$

We have proven the assertions of Proposition 8 as desired. ∎

Finally, the proof of Theorem 7 follows from Proposition 8 by a standard
argument.

###### Proof of Theorem 7.

For any $(\mathbf{Q},D)$ satisfying $\mathbf{Q}\in C^{0}\cap H^{1}$ and
$D\subset\subset\Omega$, take the sequence
$\\{(\mathbf{Q}_{\varepsilon},\phi_{\varepsilon})\\}$ guaranteed by
Proposition 8(b) to get

$\displaystyle E_{0}[\mathbf{Q},D]$ $\displaystyle\geq\limsup_{\varepsilon\to
0}E_{\varepsilon}[\mathbf{Q}_{\varepsilon},\phi_{\varepsilon}]$
$\displaystyle\geq\liminf_{\varepsilon\to
0}E_{\varepsilon}[\mathbf{Q}_{\varepsilon}^{*},\phi_{\varepsilon}^{*}]\geq
E_{0}[\mathbf{Q}_{0},D_{0}].$

The second inequality comes from the assumption that
$(\mathbf{Q}_{\varepsilon}^{*},\phi_{\varepsilon}^{*})$ minimizes
$E_{\varepsilon}$, and the third inequality utilizes Proposition 8(a) at
$(\mathbf{Q}_{0},D_{0})$. Hence, by definition $(\mathbf{Q}_{0},D_{0})$
minimizes $E_{0}$ among all $(\mathbf{Q},D)$ of concern. ∎

##  6 Numerical results

In this section, we study the numerical behaviour of the diffuse-interface LdG
model.

To simplify implementation and demonstration, and to facilitate further
comparison with the numerical results in literature, we use the 2D reduced LdG
model on 2D $\tilde{\Omega}$. The reduced LdG model is the thin film limit of
the LdG model in 3D domain $\Omega=\tilde{\Omega}\times[0,h]$, as $h\to 0$,
under certain boundary conditions [5]. Under the special reduced temperature
$A=-\frac{B^{2}}{3C}$, the critical state $\mathbf{Q}$ of the LdG free energy
(2.2) can be related to the reduced $2\times 2$ symmetric traceless tensor
$\mathbf{P}$ as follows [23]

(6.1)
$\mathbf{Q}\bumpeq\begin{bmatrix}\mathbf{P}+\frac{s_{+}}{6}\mathbf{I}_{2}&0\\\
0&-\frac{s_{+}}{3}\end{bmatrix}.$

From now on, we return to the energy $\overline{E}_{\bar{\varepsilon}}$ (2.9)
and recommence the discrimination between quantities with and without bars.
The corresponding reduced energy is then obtained through shifting
$\overline{E}_{\bar{\varepsilon}}$ (2.9) by an additive constant.

(6.2) $\displaystyle\widetilde{E}_{\bar{\varepsilon}}[\mathbf{P},\phi]$
$\displaystyle=\int_{\tilde{\Omega}}\quantity[\frac{1}{2}|\gradient\mathbf{P}|^{2}+\bar{\lambda}^{2}\quantity(-\frac{B^{2}}{4C^{2}}\tr\mathbf{P}^{2}+\frac{1}{4}(\tr\mathbf{P}^{2})^{2})]\mathop{}\\!\mathrm{d}x$
$\displaystyle\quad+\bar{\omega}_{p}\bar{\lambda}\int_{\tilde{\Omega}}\quantity[\bar{\varepsilon}|\nabla\phi|^{2}+\bar{\varepsilon}^{-1}\phi^{2}(1-\phi)^{2}]\mathop{}\\!\mathrm{d}x+\bar{\omega}_{a}\bar{\lambda}\int_{\tilde{\Omega}}\bar{\varepsilon}\left|(\mathbf{P}+s_{+}/2)\nabla\phi\right|^{2}\mathop{}\\!\mathrm{d}x$
$\displaystyle\quad+\bar{\omega}_{v}\bar{\lambda}^{2}\int_{\tilde{\Omega}}\frac{1}{2}(1-\phi)^{2}|\mathbf{P}|^{2}\mathop{}\\!\mathrm{d}x,$
s.t.
$\displaystyle\int_{\tilde{\Omega}}\phi\mathop{}\\!\mathrm{d}x=\bar{V}_{0}.$

Recall that $\bar{\varepsilon}=\varepsilon/\lambda$ is the relative capillary
width, and the rescaled constants are defined by

$\bar{\lambda}^{2}\triangleq\frac{\lambda^{2}C}{L},\
\bar{\omega}_{p}=\frac{\omega_{p}}{\sqrt{CL}},\
\bar{\omega}_{v}=\frac{\omega_{v}}{C},\
\bar{\omega}_{a}=\frac{\omega_{a}}{\sqrt{CL}},$

so the coefficients of the trailing integrals in (6.2) are respectively
$\frac{\omega_{p}\lambda}{L},\frac{\omega_{a}\lambda}{L},\frac{\omega_{v}\lambda^{2}}{L}$.
We take $B=0.64\times 10^{4}\,\mathrm{N\cdot m^{-2}}$, $C=0.35\times
10^{4}\,\mathrm{N\cdot m^{-2}}$, and $L=4\times 10^{-11}\,\mathrm{N}$ in the
following numerical simulations, which are the physical constants of the
commonly used nematic liquid crystal material MBBA [39]. We also fix the
following hyperparameters

$\bar{\varepsilon}=\frac{\varepsilon}{\lambda}=0.005,\ \bar{V}_{0}=0.09,\
\frac{\omega_{p}}{L}=3\times 10^{7}\,\mathrm{m^{-1}},\
\frac{\omega_{v}}{L}=6\times 10^{14}\,\mathrm{m^{-2}}.$

Using the finite difference method, we discretize the 2D domain
$\tilde{\Omega}=[0,1]\times[0,1]$ with an $N\times N$ grid where $N=128$, and
minimize $E_{\varepsilon}$ numerically following the gradient flow with
respect with $\mathbf{Q}$ and $\phi$. We investigate the minimizers of the
diffuse-interface LdG energy under various domain sizes $\lambda$ and
anchoring strengths $\omega_{a}$ in Figure 2.

The radial and polar states on a 2D disc with tangential boundary condition
have been studied in [26]. The radial state is radially symmetric with a
central $+1$ defect, and the polar state has two $+1/2$ defects along the
diameter. People have also found the tactoid, which is a spindle-shaped
nematic droplet within which directors are almost parallel to each other, and
with either defects near the pointy ends of the spindle or none at all [11,
2].

In Figure 2, when $\lambda$ is small enough and $\omega_{a}$ is small, the
minimizer is a radial state with circular N-I interface and a $+1$ defect in
the centre of the nematic phase. As $\omega_{a}$ increases, the minimizer
becomes a polar state with an ellipsoidal interface and the two $+1/2$ defects
along the long axis, and larger $w_{a}$, i.e. stronger tangential anchoring,
causes defects to move towards the centre. When $\lambda$ is large, the
minimizer is tactoidal state where nematic directors are almost parallel to
each other and there are no defects.

$\omega_{a}/L$$\lambda$ Figure 2: The plots of stable states under various
values of $\lambda$ in rows from top to bottom: $0.8,1,2,5,7.5$ (unit:
$10^{-6}\,\mathrm{m}$); values of $\omega_{a}/L$ in columns from left to
right: $1,3,6,9,15,30$ (unit: $10^{7}m^{-1}$), with $\bar{\varepsilon}=0.005$,
$\bar{V}_{0}=0.09$, $\omega_{p}/L=3\times 10^{7}\,\mathrm{m^{-1}}$,
$\omega_{v}/L=6\times 10^{14}\,\mathrm{m^{-2}}$. The colour represents the
order parameter $|\mathbf{P}|^{2}$. The white contour along N-I interface
denotes the level set $\\{\phi=0.5\\}$. For the green line, the length
represents the positive eigenvalue of $\mathbf{P}$ and the direction
represents the eigenvector of $\mathbf{P}$ corresponding to the positive
eigenvalue.

##  7 Conclusion and discussions

In this paper, we propose a diffuse-interface LdG energy functional
$E_{\varepsilon}$ (2.1) for the free boundary between a nematic droplet and an
isotropic surrounding. The model incorporates a LdG $\mathbf{Q}$-tensor and a
phase field $\phi$ as order parameters, with an extra hyperparameter
$\varepsilon$ introduced to describe the width of interface.

We establish the existence (Proposition 2) of global minimizers of
$E_{\varepsilon}$. With fixed $\phi$, we have also shown that the
$\mathbf{Q}$-critical points satisfy the maximum principle (Theorem 5) and
uniqueness (Proposition 6) under certain conditions. The existence, uniqueness
and regularity properties as a whole ensure the solvability of the diffuse-
interface functional $E_{\varepsilon}$. The model is therefore applicable in
the computation of free-boundary NLCs. We also prove the sharp-interface limit
of $E_{\varepsilon}$ when $\varepsilon\to 0$ in the sense of
$\Gamma$-convergence (Theorem 7 and Proposition 8), i.e. convergence of
minimizers to minimizers. The sharp-interface limit is a functional with
respect to the nematic region $D$ and the $\mathbf{Q}$-tensor order parameter
that is also made up of LdG, void, perimeter and anchoring energies.

The numerical tests we conduct produce physically interpretable results, which
are also consistent with existing observations in other free-boundary NLC
models. In particular, we observe the transition form radial to polar and
tactoidal states when adjusting the scale $\lambda$ under different
hyperparameters.

The paper also motivates research in the future. For example, there is some
room for improvement for our theoretical results such as the maximum principle
and uniqueness, which is only proven in the fixed-$\phi$ case under smoothness
assumptions. It could be generalized to global minimizers under weaker
conditions. The numerical results in Figure 2 have demonstrated exceedingly
complex structures within the minimizers of $E_{\varepsilon}$ as parameters
change. Therefore, it is necessary to further investigate into the structure
of minimizers under different parameters. One possible approach is to study
its energy landscape, i.e. multiple minimizers and the links between them,
using the method of saddle dynamics [64, 22, 50]. The mesh width of the finite
difference method has also limited our investigation into the sharp-interface
limit, so more advanced numerical schemes, such as the spectral method, should
also be implemented in the future. They will also be more useful for 3D
problems.

##  Appendix A Calculation of Euler-Lagrange equation

Consider a perturbation $\delta\mathbf{Q}$ to the function $\mathbf{Q}\in
H_{0}^{1}(\Omega;\mathcal{S}_{0}),$ and the variational derivative
$\frac{\delta E_{\varepsilon}}{\delta\mathbf{Q}}$ is an
$\mathcal{S}_{0}$-valued function that satisfies

$\delta E_{\varepsilon}=\int_{\Omega}\frac{\delta
E_{\varepsilon}}{\delta\mathbf{Q}}:\delta\mathbf{Q}\mathop{}\\!\mathrm{d}x+o(\|\delta\mathbf{Q}\|_{L^{2}}),$

where $\mathbf{A}:\mathbf{B}$ is the matrix inner product. Taking the
variational derivative of $E_{\varepsilon}$ with respect to $\mathbf{Q}$
results in the following Euler-Lagrange equation.

$\frac{\delta E^{\rm LdG}}{\delta\mathbf{Q}}+\omega_{v}\lambda^{2}\frac{\delta
E^{\rm void}}{\delta\mathbf{Q}}+\omega_{a}\lambda^{2}\frac{\delta
E_{\varepsilon}^{\rm anch}}{\delta\mathbf{Q}}=0.$

The mixing energy (2.6) is ignored because it does not involve $\mathbf{Q}$.

The variational derivative of the LdG energy (2.2) has been computed in [39,
Eq. (3.24)] as follows (notation adapted to our form).

(A.1) $\frac{\delta E^{\rm
LdG}}{\delta\mathbf{Q}}=-\triangle\mathbf{Q}+\lambda^{2}\nabla_{\mathcal{S}_{0}}\bar{F}_{b}(\mathbf{Q}),$

where the Laplacian operator $\triangle$ is defined by $\triangle
f=\sum_{i=1}^{3}f_{x_{i}x_{i}},$ and

$\nabla_{\mathcal{S}_{0}}\bar{F}_{b}(\mathbf{Q})=\mathscr{P}_{0}\quantity[\quantity(\frac{A}{C}+\tr\mathbf{Q}^{2})\mathbf{Q}-\frac{B}{C}\mathbf{Q}^{2}]=\quantity(\frac{A}{C}+\tr\mathbf{Q}^{2})\mathbf{Q}-\frac{B}{C}\quantity(\mathbf{Q}^{2}-\frac{\tr\mathbf{Q}^{2}}{3}\mathbf{I}).$

In the formulae above, $\nabla_{\mathcal{S}_{0}}$ is the projected gradient
operator on $\mathcal{S}_{0}$, which maps the matrix gradient with the
projection operator $\mathscr{P}_{0}$ (2.15). The projection is necessary
because of the traceless constraint on $\mathbf{Q}$.

With the help of $\nabla_{\mathcal{S}_{0}}$, we get that the variational
derivative of the void penalty (2.8) equals

$\frac{\delta E^{\rm
void}}{\delta\mathbf{Q}}=\omega_{v}(1-\phi)^{2}\nabla_{\mathcal{S}_{0}}\quantity(\frac{1}{2}|\mathbf{Q}|^{2}),$

and that the variational derivative of the anchoring penalty (2.7) equals

$\frac{\delta E_{\varepsilon}^{\rm
anch}}{\delta\mathbf{Q}}=\omega_{a}\varepsilon\nabla_{\mathcal{S}_{0}}\quantity[\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\nabla\phi\right|^{2}].$

We can then compute

(A.2) $\frac{\delta E^{\rm
void}}{\delta\mathbf{Q}}=\omega_{v}(1-\phi)^{2}\mathbf{Q}$

and

(A.3) $\frac{\delta E_{\varepsilon}^{\rm
anch}}{\delta\mathbf{Q}}=\omega_{a}\varepsilon\mathscr{P}_{0}\quantity[(\nabla\phi\otimes\nabla\phi)\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})+\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})(\nabla\phi\otimes\nabla\phi)],$

where the matrices
$\nabla\phi\otimes\nabla\phi,\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I}$ are
multiplied with the (non-commutative) matrix multiplication. During the
computation of (A.3), we have used the transcription of
$|(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\nabla\phi|^{2}$ into a tensor
contraction

$\left|\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})\nabla\phi\right|^{2}=(\nabla\phi\otimes\nabla\phi):\left(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I}\right)^{2},$

and an obvious result on the matrix gradient of polynomial functions.

######  Lemma A.1.

If $f(\mathbf{Q})=\mathbf{A}:\mathbf{Q}^{k}$ where $k\in\mathbb{N}$, then the
matrix gradient of $f$ is

$\nabla
f(\mathbf{Q})=\sum_{j=0}^{k-1}\mathbf{Q}^{j}\mathbf{A}\mathbf{Q}^{k-1-j}.$

Summing up (A.1), (A.2) and (A.3), we get the equation (2.13).

Taking the first variation of $E_{\varepsilon}$ with respect to $\phi$, we get

$\frac{\delta E_{\varepsilon}^{\rm
mix}}{\delta\phi}+\frac{\omega_{a}}{\omega_{p}}\frac{\delta
E_{\varepsilon}^{\rm
anch}}{\delta\phi}+\frac{\omega_{v}\lambda}{\omega_{p}}\frac{\delta E^{\rm
void}}{\delta\phi}+\xi=0,$

where $\xi$ is a Lagrangian multiplier from the volume constraint (2.10). The
LdG energy $E^{\rm LdG}$ (2.2) is ignored because it does not involve $\phi$.
We can then directly compute that

$\displaystyle\frac{\delta E_{\varepsilon}^{\rm mix}}{\delta\phi}$
$\displaystyle=-2\varepsilon\triangle\phi+2\varepsilon^{-1}\phi(\phi-1)(2\phi-1),$
$\displaystyle\frac{\delta E_{\varepsilon}^{\rm anch}}{\delta\phi}$
$\displaystyle=-2\varepsilon\operatorname{\rm
div}\quantity[\quantity(\mathbf{Q}+\frac{s_{+}}{3}\mathbf{I})^{2}\nabla\phi],$
$\displaystyle\frac{\delta E^{\rm void}}{\delta\phi}$
$\displaystyle=(\phi-1)|\mathbf{Q}|^{2}.$

Summing up these formulae, and dividing the entire equation by $2\varepsilon$,
we get the equation (2.14). Note that the Lagrangian multiplier $\xi$ has been
rescaled as well, but we still write it as $\xi$ to simplify notation.

##  Appendix B Proof of Lemma 9

The proof of the lower bound in Lemma 9 follows exactly [42, Prop. 1], which
employs the elementary inequality $a^{2}+b^{2}\geq 2ab$ and the lower semi-
continuity of the total variation. In fact, we prove that the lower bound
holds for any $L^{1}$-convergent sequence $\\{\phi_{\varepsilon}\\}$
regardless of the volume constraint (5.14).

First, we assume w.l.o.g. that $\alpha\leq\phi_{\varepsilon}\leq\beta$, for
otherwise we can define the truncated function
$\tilde{\phi}_{\varepsilon}(x)=\tau_{1}(\phi_{\varepsilon})$, where

$\tau_{1}(x)=\max\\{\min\\{x,\beta\\},\alpha\\}=\begin{cases}\alpha,&x<\alpha,\\\
x,&\alpha\leq x\leq\beta,\\\ \beta,&x>\beta.\end{cases}$

It is obvious that
$|\tilde{\phi}_{\varepsilon}-\phi_{0}|\leq|\phi_{\varepsilon}-\phi_{0}|,$ so
$\|\tilde{\phi}_{\varepsilon}-\phi_{0}\|_{L^{1}}\to 0$, i.e.
$\\{\tilde{\phi}_{\varepsilon}\\}$ converges in $L^{1}$ as well. In addition,
it follows from the property of $W$ and the chain rule that

$W(\tilde{\phi}_{\varepsilon})\leq
W(\phi_{\varepsilon}),\quad\nabla\tilde{\phi}_{\varepsilon}=I_{\\{\alpha<\phi_{\varepsilon}<\beta\\}}\nabla\phi_{\varepsilon},$

so

$\int_{\Omega}[\varepsilon^{-1}W(\tilde{\phi}_{\varepsilon})+\varepsilon
a(x,\nabla\tilde{\phi}_{\varepsilon})]\mathop{}\\!\mathrm{d}x\leq\int_{\Omega}[\varepsilon^{-1}W(\phi_{\varepsilon})+\varepsilon
a(x,\nabla\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x,$

i.e. $L^{a}_{\varepsilon}(\tilde{\phi}_{\varepsilon})\leq
L^{a}_{\varepsilon}(\phi_{\varepsilon})$. Therefore, in order to prove Lemma
9, it suffices to check that

$\liminf_{\varepsilon\to 0}L^{a}_{\varepsilon}(\tilde{\phi}_{\varepsilon})\geq
2c_{0}\int_{\partial D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S,$

where the sequence $\\{\tilde{\phi}_{\varepsilon}\\}$ also has the $L^{1}$
convergence, and comes with the extra assumption that
$\alpha\leq\tilde{\phi}_{\varepsilon}\leq\beta$.

Assume that $\alpha\leq\phi_{\varepsilon}\leq\beta$ and discard the volume
constraint. We can estimate the LHS of Lemma 9 with the elementary inequality
$a^{2}+b^{2}\geq 2ab.$ This step cancels out $\varepsilon$.

$\int_{\Omega}\quantity[\varepsilon^{-1}W(\phi_{\varepsilon})+\varepsilon
a(x,\nabla\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x\geq
2\int_{\Omega}\quantity[W(\phi_{\varepsilon})a(x,\nabla\phi_{\varepsilon})]^{\frac{1}{2}}\mathop{}\\!\mathrm{d}x.$

Define $\Phi(s)=\int_{\alpha}^{s}\sqrt{W(s)}\mathop{}\\!\mathrm{d}s$, so
$\Phi(\alpha)=0,\Phi(\beta)=c_{0}$ and $\Phi$ is strictly increasing on
$(\alpha,\beta)$. Then by the chain rule we have
$\nabla\Phi(\phi_{\varepsilon})=\Phi^{\prime}(\phi_{\varepsilon})\nabla\phi_{\varepsilon}=\sqrt{W(\phi_{\varepsilon})}\nabla\phi_{\varepsilon},$
which enables us to rewrite

$\displaystyle\int_{\Omega}\quantity[W(\phi_{\varepsilon})a(x,\nabla\phi_{\varepsilon})]^{\frac{1}{2}}\mathop{}\\!\mathrm{d}x$
$\displaystyle=\int_{\Omega}a\quantity(x,\sqrt{W(\phi_{\varepsilon})}\nabla\phi_{\varepsilon})^{\frac{1}{2}}\mathop{}\\!\mathrm{d}x$
$\displaystyle=\int_{\Omega}a(x,\nabla\Phi(\phi_{\varepsilon}))^{\frac{1}{2}}\mathop{}\\!\mathrm{d}x.$

Hence, we get the inequality

(B.1)
$\int_{\Omega}\quantity[\varepsilon^{-1}W(\phi_{\varepsilon})+\varepsilon
a(x,\nabla\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x\geq
2\int_{\Omega}a(x,\nabla\Phi(\phi_{\varepsilon}))^{\frac{1}{2}}\mathop{}\\!\mathrm{d}x.$

Assume for the moment that $a^{ij}\in C^{1}(\Omega)$. We define the following
functional inspired by the total variation [15].

(B.2)
$V(\varphi)\triangleq\sup\left\\{\int_{D}\phi\mathrm{D}_{i}(a^{ij}g_{j})\mathop{}\\!\mathrm{d}x:\bm{g}=[g_{j}]\in
C_{0}^{1}(\Omega;\mathbb{R}^{n}),\sup_{x}a(x,\bm{g})\leq 1\right\\},$

where Einstein’s convention is used for $\mathrm{D}_{i}(a^{ij}g_{j})$.
Obviously the set on the RHS is non-empty as $\bm{g}\equiv 0$ is an element.
Thus, $V(\varphi)$ is nonnegative (equals $+\infty$ if the set is unbounded).

We consider the following properties obvious, as the proof is identical to
corresponding properties of the total variation.

######  Lemma B.1.

We claim the following properties regarding the functional
$V:L^{1}(\Omega)\to\mathbb{R}_{\geq 0}\cup\\{+\infty\\}$.

  1. (a)

$V(\varphi)$ is sublinear and homogeneous, i.e.
$V(a\varphi+\psi)\leq|a|V(\varphi)+V(\psi)$ for all
$a\in\mathbb{R},\varphi,\psi\in L^{1}(\Omega)$.

  2. (b)

$V(\varphi)=0$ if $\varphi\equiv\text{const.}$, and as a consequence,
$V(\varphi+C)=V(\varphi)$ for all $\varphi\in L^{1}$ and constant $C$.

  3. (c)

$V(\varphi)$ is weakly lower semi-continuous on $L^{1}(\Omega)$, i.e. for any
sequence $\varphi_{k}\to\varphi$ in $L^{1}$,

$\liminf_{k\to\infty}V(\varphi_{k})\geq V(\varphi).$

  4. (d)

$V(\varphi)=\int_{\Omega}|a(x,\nabla\varphi)|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}x$
whenever $\varphi\in W^{1,1}(\Omega)$.

  5. (e)

$V(\varphi)=\int_{\partial
D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S$ whenever
$\varphi=I_{D}$ where $D\subset\subset D$ is Lipschitz.

Since $\Phi$ is Lipschitz continuous on $[\alpha,\beta]$ (by
$\Phi^{\prime}=\sqrt{W}$) and that $\phi_{\varepsilon}\to\phi_{0}$ in $L^{1}$,
$\Phi(\phi_{\varepsilon})\to\Phi(\phi_{0})$ in $L^{1}$ as well. As
$\phi_{0}=\alpha I_{D}+\beta I_{D^{c}}$ by definition, we get that
$\Phi(\phi_{0})=c_{0}(1-I_{D}).$ We then make the following deduction based on
Lemma B.1.

(B.3) $\displaystyle\liminf_{\varepsilon\to
0}\int_{\Omega}|a(x,\nabla\Phi(\phi_{\varepsilon}))|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}x$
$\displaystyle\overset{\rm(d)}{=}\liminf_{\varepsilon\to
0}V(\Phi(\phi_{\varepsilon}))$
$\displaystyle\overset{\rm(c)}{\geq}V(c_{0}-c_{0}I_{D})\overset{\rm(b)}{=}V(-c_{0}I_{D})$
$\displaystyle\overset{\rm(a)}{=}c_{0}V(I_{D})\overset{\rm(e)}{=}c_{0}\int_{\partial
D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S.$

Explanation of (B.3): the first equality uses claim (d) on differentiable
$\varphi$; the second inequality uses claim (c) regarding w.l.s.c. property;
the third equality uses claim (b) by shifting with a constant; the fourth
equality uses the homogeneity from claim (a); the final equality uses claim
(e) on indicator functions.

The result of Lemma 9 then follows by connecting the inequalities (B.1) and
(B.3).

$\liminf_{\varepsilon\to
0}\int_{\Omega}\quantity[\varepsilon^{-1}W(\phi_{\varepsilon})+\varepsilon
a(x,\nabla\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x\geq
2\liminf_{\varepsilon\to 0}V(\Phi(\phi_{\varepsilon}))\geq
2c_{0}\int_{\partial D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S.$

We have assumed $\mathbf{D}(x)\in C^{1}$ thus far. As for more general cases
where $\mathbf{A}(x)\in C^{0}$, we consider a uniform approximation
$\mathbf{A}_{k}(x)\to\mathbf{A}(x)$ “from below”, that is, satisfying
$\mathbf{A}_{k}\preceq\mathbf{A}$ everywhere. Then by the property of
positive-definite matrices,

$a_{k}(x,p)\leq a(x,p),\quad\forall x\in\Omega,p\in\mathbb{R}^{n},$

where $a_{k}(x,p)=a_{k}^{ij}(x)p_{i}p_{j}$ is the quadratic form defined by
$\mathbf{A}_{k}(x)$. Because of the uniform ellipticity (5.12), construction
of such an approximation is simple and hence omitted.

Applying the previous argument to the sequence $\mathbf{A}_{k}$, we then get

$\displaystyle\liminf_{\varepsilon\to
0}\int_{\Omega}\quantity[\varepsilon^{-1}W(\phi_{\varepsilon})+\varepsilon
a(x,\nabla\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x$
$\displaystyle\geq\liminf_{\varepsilon\to
0}\int_{\Omega}\quantity[\varepsilon^{-1}W(\phi_{\varepsilon})+\varepsilon
a_{k}(x,\nabla\phi_{\varepsilon})]\mathop{}\\!\mathrm{d}x$ $\displaystyle\geq
2c_{0}\int_{\partial D}|a_{k}(x,\bm{\nu})|\mathop{}\\!\mathrm{d}S(x).$

Let $k\to\infty$ in the inequality above, and the boundary integral on the RHS
converges uniformly to $\int_{\partial D}|a(x,\bm{\nu})|^{\frac{1}{2}},$
concluding the proof.

##  Appendix C Proof of Lemma 10

The proof of the upper bound follows [42, Prop. 2], which uses an explicit
construction of $\\{\phi_{\varepsilon}\\}$ aided by the signed distance
function (SDF) [45, Chap. 2].

If the upper limit were to hold, then the limit of
$L_{\varepsilon}^{a}(\phi_{\varepsilon})$ should be exactly equal to the RHS
because of the lower bound estimates by Lemma 9. Therefore, inequalities in
the proof of the lower bound must approximate equalities. Especially, for the
elementary inequality in (B.1) to become an equality, it must hold almost
everywhere that

$\varepsilon^{-1}W(\phi_{\varepsilon})\approx\varepsilon
a(x,\nabla\phi_{\varepsilon}).$

As a result, the integrations of the two terms are approximately equal,

$\int_{\Omega}\varepsilon^{-1}W(\phi_{\varepsilon})\approx\int_{\Omega}a(x,\nabla\phi_{\varepsilon}).$

This property is known as the equipartition of energy [8].

For the generalized model (5.11), we propose the following motivations for our
proof.

  1. i.

By constructing a function $\chi_{\varepsilon}:\mathbb{R}\to[\alpha,\beta]$
with
$\chi_{\varepsilon}^{\prime}\approx\varepsilon^{-1}\sqrt{W(\chi_{\varepsilon})}$
and a function $h^{a}\in W^{1,\infty}$, one obtains with the chain rule that

$\varepsilon
a(x,\nabla\chi_{\varepsilon}(h^{a}))\approx\varepsilon\cdot\varepsilon^{-2}W(\chi_{\varepsilon}(h^{a}))a(x,\nabla
h^{a})=\varepsilon^{-1}W(\chi_{\varepsilon}(h^{a})),$

if $a(x,\nabla h^{a})\equiv 1$.

  2. ii.

In order to fulfil $a(x,\nabla h^{a})\equiv 1$, we solve a nonlinear first-
order partial differential equation with the method of characteristics. The
solution we get is a generalized SDF.

  3. iii.

The level sets of $h^{a}$ are locally smooth mappings of the boundary
$\partial D$ near $h^{a}=0$, so we evaluate
$L_{\varepsilon}^{a}(\phi_{\varepsilon})$ with the coarea formula, and
complete the proof.

###  C.1 Approximate equipartition of energy

We briefly reproduce the derivation of the function $\chi_{\varepsilon}$ by
Modica [42]. The function $\chi_{\varepsilon}$ is intended to approximate the
relation
$\chi_{\varepsilon}^{\prime}\approx\varepsilon^{-1}\sqrt{W(\chi_{\varepsilon})},$
to fulfil the equipartition condition.

We propose for $\chi_{\varepsilon}:\mathbb{R}\to\mathbb{R}$ as the solution to
the following differential equation.

(C.1)
$\begin{cases}\varepsilon\chi_{\varepsilon}^{\prime}(t)=\sqrt{W(\chi_{\varepsilon}(t))+\varepsilon},\\\
\chi_{\varepsilon}(0)=\alpha.\end{cases}$

$\chi_{\varepsilon}$ is intended to approximate the standing-wave solution
$q$, which happens to be the minimizer of the one-dimensional van der Waals-
Cahn-Hilliard energy satisfying
$\varepsilon^{2}q^{\prime\prime}=\frac{1}{2}W(q)$ [66].

The following properties of $\chi_{\varepsilon}$ are evident.

######  Lemma C.1.

The function $\chi_{\varepsilon}$ satisfies that:

  1. (a)

$\chi_{\varepsilon}$ is strictly increasing.

  2. (b)

There exists an interval $[0,\eta_{\varepsilon}]$, such that
$\chi_{\varepsilon}(0)=\alpha,\chi_{\varepsilon}(\eta_{\varepsilon})=\beta$.

  3. (c)

$\eta_{\varepsilon}\to 0$ as $\varepsilon\to 0$.

###  C.2 Method of characteristics

The next step is to construct (locally) a Lipschitz function $h^{a}$
satisfying

$a(x,\nabla h^{a})=1,\quad\text{a.e.}.$

$h^{a}$ is essentially the solution to the following first-order nonlinear
partial differential equation (PDE)

(C.2) $\begin{cases}a^{ij}(x)u_{x_{i}}u_{x_{j}}=1,&x\in\Omega\\\
u=0,&x\in\partial D,\end{cases}$

which can be solved with the method of characteristics. Assume for the moment
that $a^{ij}(x)\in C^{\infty}(\overline{\Omega})$. Then by [14, sect. 3.2]
(with $F(x,u,p)=a(x,p)-1$), the solution to (C.2) is given by the
characteristic curves $(\bm{X}(s),U(s),\bm{P}(s))$ satisfying the dynamical
system

(C.3a)
$\begin{cases}\dot{X}_{i}=-\partial_{p_{i}}F(\bm{X},U,\bm{P})=2a^{ij}(\bm{X})P_{j},\\\
\dot{U}=\nabla_{p}F(\bm{X},U,\bm{P})\cdot\bm{P}=2a^{ij}(\bm{X})P_{i}P_{j},\\\
\dot{P}_{i}=-\partial_{x_{i}}F(\bm{X},U,\bm{P})-P_{i}\partial_{u}F(\bm{X},U,\bm{P})=-a^{jk}_{x_{i}}(\bm{X})P_{j}P_{k},\end{cases}\quad\quantity(\dot{(*)}=\frac{\partial(*)}{\partial
s})$ (Einstein’s summation convention is used) with initial values (C.3b)
$\begin{cases}\bm{X}(0)=y\in\partial D,\\\ U(0)=0,\\\
\bm{P}(0)=|a(y,\bm{\nu})|^{-\frac{1}{2}}\bm{\nu}.\end{cases}$ The curve
$\bm{X}(s)$ originates from the manifold of boundary conditions, which in this
case is $\partial D$, and the functions $U(s),\bm{P}(s)$ picture respectively
the evolution of $h^{a}(x)$ and $\nabla h^{a}(x)$ along the curves. The
initial values (C.3b) are designed to be compatible with the boundary
condition $h^{a}=0$ and $a(x,\nabla h^{a})=1$.

With the general theory of first-order nonlinear PDE’s, we assert that the
system (C.3a) solves (C.2) correctly in a neighbourhood of $\partial D$, with
$U(t),\bm{P}(t)$ standing for the respective value of $h^{a},\nabla h^{a}$ at
the point $\bm{X}(t)$.

######  Lemma C.2.

Assume that $\partial D$ is compact and $C^{2}$ smooth. Then, there exists an
open interval $I$ containing 0 such that:

  1. (a)

For each $y\in\partial D$ and $s\in I$, a unique $C^{1}$ characteristic curve
$(\bm{X}(s),U(s),\bm{P}(s))$ exists with the initial conditions (C.3b).

  2. (b)

The mapping

$\displaystyle I\times\partial D$
$\displaystyle\to\Omega\times\mathbb{R}\times\mathbb{R}^{n}$
$\displaystyle(s,y)$ $\displaystyle\mapsto(\bm{X},U,\bm{P})$

is $C^{1}$ with respect to $s\in I$ and $y\in\partial D$.

  3. (c)

The mapping $(s,y)\mapsto\bm{X}(s;y)$ for $s\in I,y\in\partial D$ is
bijective, and its inverse is also $C^{1}$, i.e. there exists $C^{1}$ inverse
mappings $s=s(x),y=y(x)$.

  4. (d)

On the open region containing $D$ defined by

$V=\left\\{\bm{X}(s;y):s\in I,y\in\partial D\right\\},$

the function

$h^{a}(x)=U(s(x);y(x))$

is a $C^{1}$ classical solution to the equation (C.2), with $\nabla
h^{a}(x)=\bm{P}(s(x);y(x))$.

We can also make the following observations on the solution to (C.3a) on the
neighbourhood $V$. It states that the sign of $h^{a}(x)$ determines which side
of $\partial D$ the point $x$ is on, a property shared by the classical SDF.

######  Lemma C.3.

Assume the same as Lemma C.2, and let the solution to (C.2) defined by
characteristic curves be $u(x)$. Then for any $x\in V$, the relative position
of $x$ to $D$ is determined by the sign of $h^{a}(x)$. Specifically,

$\begin{cases}x\in D,&h^{a}(x)<0,\\\ x\in\partial D,&h^{a}(x)=0,\\\
x\in(\overline{D})^{c},&h^{a}(x)>0.\end{cases}$

###### Proof.

Since $U(s)=h^{a}(\bm{X}(s)),\bm{P}(s)=\nabla h^{a}(\bm{X}(s))$ solves (C.2)
by Lemma C.2,

$a(\bm{X}(s),\bm{P}(s))\equiv 1$

for all $s\in I$. Therefore, we notice from the characteristic equation (C.3a)
that $\dot{U}(s)=2a(\bm{X},\bm{P})=2,$ which indicates

(C.4) $h^{a}(x)=U(s(x))=2s.$

By (C.3b), the initial values satisfy
$\dot{\bm{X}}(0)=2a^{ij}(y)P_{j}(0),\bm{P}(0)=c\bm{\nu}$ ($c>0$), so

$\dot{\bm{X}}(0)\cdot\bm{\nu}=c2a^{ij}(y)\nu_{i}\nu_{j}>0$

since $a^{ij}$ satisfies the elliptic condition (5.12). As $\bm{\nu}$ is the
exterior normal vector, the characteristic curve points outside the region $D$
as $s$ increases. Therefore, $s>0$ corresponds to the outside of $D$ and $s<0$
corresponds to the inside of $D$. By (C.4), the signs of $h^{a}$ and $s(x)$
are identical, and so we yield Lemma C.3. ∎

Level sets of $h^{a}$ are local continuous isomorphisms of the boundary
$\partial D$.

######  Corollary C.4.

Use the same assumptions and notations as Lemma C.2. The level sets
$\Sigma_{t}=\\{h^{a}=t\\}$ of $h^{a}(x)$ satisfies the following properties.

  1. (a)

For all $t$ with $\frac{t}{2}\in I$, $\Sigma_{t}$ is a $C^{1}$ hypersurface.

  2. (b)

For any continuous function $g(x):\Omega\to\mathbb{R}$, the integral

$\int_{\Sigma_{t}}g(x)\mathop{}\\!\mathrm{d}S(x)$

is continuous with respect to $t$. Specifically, when $g(x)\equiv 1$, we get
the continuity of $\mathcal{H}_{n-1}(\Sigma_{t})$ with respect to $t$.

###### Proof.

Recall from (C.4) that $U(s)=2s$, so the level set $h^{a}=t$ is expressed
equivalently as

$h^{a}(x)=t\Leftrightarrow s(x)=\frac{t}{2}.$

That is, $\Sigma_{t}$ is a section of the $C^{1}$ bijective mapping
$(s,y)\mapsto\bm{X}(s;y)$ (by Lemma C.2(c)) at $s=\frac{t}{2}.$

$\Sigma_{t}=\\{\bm{X}(t/2;y):y\in\partial D\\},\quad\forall\frac{t}{2}\in I.$

Therefore, $s=\frac{t}{2}$ is a $C^{1}$ parametrization of the surface family
$\\{\Sigma_{t}\\}$, and the assertions of Corollary C.4 follow immediately. ∎

###  C.3 Integration with coarea formula

Work on the interval $I$ determined by Lemma C.2. W.l.o.g., assume that
$I=(-\eta_{*},\eta_{*})$.

Since $\eta_{\varepsilon}\to 0$ as $\varepsilon\to 0$ by Lemma C.1, we choose
$\varepsilon$ small enough such that $\eta_{\varepsilon}<\eta_{*}$. Then
according to Corollary C.4, for all
$t\in[-\eta_{\varepsilon},\eta_{\varepsilon}]$, the level set $\Sigma_{t}$
resides within the neighbourhood $V$. We define on $V$

(C.5)
$\phi_{\varepsilon}(x)=\begin{cases}\chi_{\varepsilon}(h^{a}(x)+\delta_{\varepsilon}),&-\delta_{\varepsilon}\leq
h^{a}(x)\leq\eta_{\varepsilon}-\delta_{\varepsilon},\\\
\alpha,&h^{a}<-\delta_{\varepsilon},\\\
\beta,&h^{a}>\eta_{\varepsilon}-\delta_{\varepsilon},\end{cases}$

where $\delta_{\varepsilon}\in[0,\eta_{\varepsilon}]$ is yet to be determined.
Effectively, we have applied a truncated function
$\phi_{\varepsilon}=\tilde{\chi}_{\varepsilon}(h^{a}+\delta_{\varepsilon})$,
where

$\tilde{\chi}_{\varepsilon}(t)=\min\left\\{\beta,\max\\{\alpha,\chi_{\varepsilon}(t)\\}\right\\}=\begin{cases}\alpha,&t<0,\\\
\chi_{\varepsilon}(t),&t\in[0,\eta_{\varepsilon}],\\\
\beta,&t>\eta_{\varepsilon}\end{cases}$

is an increasing and Lipschitz continuous function from
$\mathbb{R}\to[\alpha,\beta]$, so $\phi_{\varepsilon}$ is also Lipschitz
continuous on $V$.

$t$$\tilde{\chi}_{\varepsilon}(t)$$\chi_{0}(t)$$\eta_{\varepsilon}$$0$ Figure
3: Sketch of $\tilde{\chi}_{\varepsilon}(t)$ compared with $\chi_{0}$

We choose $\delta_{\varepsilon}\in[0,\eta_{\varepsilon}]$ to meet the volume
constraint (5.14).

######  Lemma C.5.

There exists $\delta_{\varepsilon}\in[0,\eta_{\varepsilon}]$ such that

$\int_{V}\phi_{\varepsilon}(x)\mathop{}\\!\mathrm{d}x=\int_{V}\phi_{0}(x)\mathop{}\\!\mathrm{d}x,$

where $\phi_{0}=\alpha I_{D}+\beta I_{D^{c}}$.

###### Proof.

Define $\chi_{0}=\begin{cases}\alpha,&t<0\\\ \beta,&t>0\end{cases}$, which
satisfies

$\tilde{\chi}_{\varepsilon}(t)\leq\chi_{0}(t)\leq\tilde{\chi}_{\varepsilon}(t+\eta_{\varepsilon}),\quad\forall
t\in\mathbb{R}$

according to Figure 3. According to Lemma C.3(c), the sign of
$h^{a}(x)=U(s)=2s$ is determined by $x\in D^{c}$ or $x\in D$, so

$\chi_{0}(h^{a}(x))=\begin{cases}\alpha,&x\in D\\\ \beta,&x\in
D^{c}\end{cases}=\phi_{0}(x),$

where $\phi_{0}=\alpha I_{D}+\beta I_{D^{c}}$ is the sharp-interface limit in
Lemmas 9 and 10. Since the integral

$\int_{V}\phi_{\varepsilon}(x)\mathop{}\\!\mathrm{d}x=\int_{V}\tilde{\chi}_{\varepsilon}(h^{a}(x)+\delta)\mathop{}\\!\mathrm{d}x$

is continuous and monotonic with respect to $\delta$, and since

$\int_{V}\tilde{\chi}_{\varepsilon}(h^{a}(x))\mathop{}\\!\mathrm{d}x\leq\int_{V}\chi_{0}(h^{a}(x))\mathop{}\\!\mathrm{d}x=m\leq\int_{V}\chi_{\varepsilon}(h^{a}(x)+\eta_{\varepsilon})\mathop{}\\!\mathrm{d}x,$

there exists $\delta_{\varepsilon}\in[0,\eta_{\varepsilon}]$ such that
$\int_{V}\tilde{\chi}_{\varepsilon}(h^{a}+\delta_{\varepsilon})\mathop{}\\!\mathrm{d}x=m$
by the mean value theorem. ∎

We then extend $\phi_{\varepsilon}$ to the entire $\Omega$. Define

(C.6) $\phi_{\varepsilon}(x)=\begin{cases}\alpha,&x\in D\cap V^{c},\\\
\beta,&x\in D^{c}\cap V^{c}.\end{cases}$

We claim that:

  1. (i)

The extension (C.6) is a Lipschitz function over $\Omega$.

  2. (ii)

$\phi_{\varepsilon}$ satisfies the volume constraint (2.10).

  3. (iii)

$\phi_{\varepsilon}\to\phi_{0}$ a.e., where $\phi_{0}=\alpha I_{D}+\beta
I_{D^{c}}$.

To prove claim (i), we just need to prove the differentiability near the
boundary $\partial V$. Since $V$ is defined by images of the $C^{1}$ mapping
$\bm{X}(s;y)$ with $s\in(-\eta_{*},\eta_{*})$, the boundary $\partial V$
equals

$\partial V=\\{\bm{X}(s;y):s=\pm\eta_{*},y\in\partial D\\},$

which is composed of two patches of level set surfaces $\Sigma_{\pm
2\eta_{*}}=\\{h^{a}=\pm 2\eta_{*}\\}$. It holds on either of the patches that
$|h^{a}(s)|>\eta_{\varepsilon}$ in a neighbourhood, so for all $x\in V$ in
that neighbourhood, $\phi_{\varepsilon}(x)=\tilde{\chi}_{\varepsilon}(h^{a})$
is constantly $\alpha$ or $\beta$ according as $h^{a}<0$ or $h^{a}>0$. On the
other hand, the extension (C.6) on $x\in V^{c}$ also equals $\alpha$ or
$\beta$ according as $x\in D$ or $x\in D^{c}$, which is equivalent to $x\in D$
or $x\in D^{c}$ according to Lemma C.3. Therefore, the extended
$\phi_{\varepsilon}$ is constant near $\partial V$, and is obviously Lipschitz
continuous.

To prove claim (ii), we notice that $\phi_{0}$ is also constantly $\alpha$ or
$\beta$ according as $x\in D$ or $x\in D^{c}$ satisfies the relations in
(C.6), which ensures

$\int_{V^{c}}\phi_{\varepsilon}(x)\mathop{}\\!\mathrm{d}x=\int_{V^{c}}\phi_{0}(x)\mathop{}\\!\mathrm{d}x.$

Meanwhile, Lemma C.5 ensures that

$\int_{V}\phi_{\varepsilon}(x)\mathop{}\\!\mathrm{d}x=\int_{V}\phi_{0}(x)\mathop{}\\!\mathrm{d}x.$

Adding them together, we get the volume constraint (2.10) on the entire
$\Omega$.

To prove claim (iii), we observe that $\phi_{\varepsilon}=\phi_{0}$ whenever
$h^{a}\notin(-\delta_{\varepsilon},\eta_{\varepsilon}-\delta_{\varepsilon}).$
Therefore, for each $x\notin\partial D$, either $x\notin
V\Rightarrow\phi_{\varepsilon}\equiv\phi_{0}$ or $h^{a}(x)\neq 0\Rightarrow$
when $\varepsilon$ is sufficiently small,
$h^{a}\notin(-\delta_{\varepsilon},\eta_{\varepsilon}-\delta_{\varepsilon})\Rightarrow\phi_{\varepsilon}(x)=\phi_{0}(x)$.
This means that $\phi_{\varepsilon}\to\phi_{0}$ a.e..

Finally, we estimate the functional (5.11) with the aid of
$\phi_{\varepsilon}$ and the following coarea formula for a Lipschitz function
$h\in W^{1,\infty}(\Omega)$ and an integrable function $f\in L^{1}(\Omega)$
[19, 14, 15], which translates volume integrals to surface integrals.

(C.7) $\int_{\Omega}f(x)|\nabla
h|\mathop{}\\!\mathrm{d}x=\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}r\int_{\\{h=r\\}\cap\Omega}f(x)\mathop{}\\!\mathrm{d}S(x).$

The following steps are identical to those of Modica [42]. When
$h^{a}\geq\eta_{\varepsilon}-\delta_{\varepsilon}$ or
$h^{a}\leq-\delta_{\varepsilon}$, $\phi_{\varepsilon}$ is constantly $\alpha$
or $\beta$, so both $W(\phi_{\varepsilon})$ and $\nabla\phi_{\varepsilon}$ are
zero. Therefore, the integral (5.11) is supported on the “strip”
$\\{-\delta_{\varepsilon}<h^{a}<\eta_{\varepsilon}-\delta_{\varepsilon}\\},$
or equivalently, $\\{\alpha<\phi_{\varepsilon}<\beta\\}.$ Using the
differential equation (C.1) and the chain rule, we get

$\varepsilon
a(x,\nabla\phi_{\varepsilon})=\varepsilon^{-1}(W(\phi_{\varepsilon})+\varepsilon)a(x,\nabla
h^{a})=\varepsilon^{-1}(W(\phi_{\varepsilon})+\varepsilon).$

Thus, we transcribe the integral into

$\displaystyle L_{\varepsilon}^{a}(\phi_{\varepsilon})$
$\displaystyle=\int_{\\{\alpha<\phi_{\varepsilon}<\beta\\}}\varepsilon^{-1}(2W(\phi_{\varepsilon})+\varepsilon)\mathop{}\\!\mathrm{d}x$
$\displaystyle=\int_{-\delta_{\varepsilon}}^{\eta_{\varepsilon}-\delta_{\varepsilon}}\mathop{}\\!\mathrm{d}t\int_{\Sigma_{t}}\frac{\varepsilon^{-1}}{|\nabla
h^{a}|}(2W(\phi_{\varepsilon})+\varepsilon)\mathop{}\\!\mathrm{d}S(x),$

where the coarea formula (C.7) is applied to $h^{a}$. Note that because of
(5.12), $1=|a(x,\nabla h^{a})|\leq\Lambda|\nabla h^{a}|^{2}\Rightarrow|\nabla
h^{a}|\geq\Lambda^{-\frac{1}{2}}$, and thus it can be placed on the
denominator. As $\phi_{\varepsilon}=\tilde{\chi}_{\varepsilon}(h^{a})$, we can
extract the factor $W(\phi_{\varepsilon})+\varepsilon$ from the surface
integral.

$\displaystyle L_{\varepsilon}^{a}(\phi_{\varepsilon})$
$\displaystyle=\int_{-\delta_{\varepsilon}}^{\eta_{\varepsilon}-\delta_{\varepsilon}}\frac{2W(\tilde{\chi}_{\varepsilon}(t+\delta_{\varepsilon}))+\varepsilon}{\varepsilon}\mathop{}\\!\mathrm{d}t\int_{\Sigma_{t}}\frac{1}{|\nabla
h^{a}|}\mathop{}\\!\mathrm{d}S(x)$
$\displaystyle=\int_{0}^{\eta_{\varepsilon}}\frac{2W(\chi_{\varepsilon}(t))+\varepsilon}{\varepsilon}\mathop{}\\!\mathrm{d}t\int_{\Sigma_{t-\delta_{\varepsilon}}}\frac{1}{|\nabla
h^{a}|}\mathop{}\\!\mathrm{d}S(x).$

As $\chi_{\varepsilon}(t)$ is strictly increasing by Lemma C.1(a), we make the
change of variable $r=\chi_{\varepsilon}(t)\in(\alpha,\beta)$, whose
substitution formula is
$\mathop{}\\!\mathrm{d}r=\chi_{\varepsilon}^{\prime}(t)\mathop{}\\!\mathrm{d}t=\varepsilon^{-1}\sqrt{W(r)+\varepsilon}\mathop{}\\!\mathrm{d}t.$
Then,

(C.8)
$L_{\varepsilon}^{a}(\phi_{\varepsilon})=\int_{\alpha}^{\beta}\frac{2W(r)+\varepsilon}{\sqrt{W(r)+\varepsilon}}\mathop{}\\!\mathrm{d}r\int_{\Sigma_{\chi_{\varepsilon}^{-1}(r)-\delta_{\varepsilon}}}\frac{1}{|\nabla
h^{a}|}\mathop{}\\!\mathrm{d}S(x),$

where
$\chi_{\varepsilon}^{-1}(r)\in(-\delta_{\varepsilon},\eta_{\varepsilon}-\delta_{\varepsilon})$
is the inverse function of $\chi_{\varepsilon}$. By Lemma C.1,
$\eta_{\varepsilon}\to 0$ and $\delta_{\varepsilon}\to 0$ as $\varepsilon\to
0$, so $\chi_{\varepsilon}^{-1}(r)-\delta_{\varepsilon}\to 0$ uniformly for
all $r\in(\alpha,\beta)$. $\frac{1}{|\nabla h^{a}|}$ is also continuous over
$V$ because $h^{a}\in C^{1}$ and $|\nabla h^{a}|\geq\Lambda^{-\frac{1}{2}}$.
Thus, by Lemma C.4(b) we can assert that

$\lim_{\varepsilon\to
0}\int_{\Sigma_{\chi_{\varepsilon}^{-1}(r)-\delta_{\varepsilon}}}\frac{1}{|\nabla
h^{a}|}\mathop{}\\!\mathrm{d}S(x)=\int_{\Sigma_{0}}\frac{1}{|\nabla
h^{a}|}\mathop{}\\!\mathrm{d}S(x)$

uniformly for all $r\in(\alpha,\beta)$. Moreover,
$\frac{2W(r)+\varepsilon}{\sqrt{W(r)+\varepsilon}}\leq
2\sqrt{W(r)+\varepsilon}$ is uniformly bounded above by a continuous function
for all $r\in[\alpha,\beta]$. Therefore, we apply the dominated convergence
theorem on (C.8) by taking the limit $\varepsilon\to 0$ under the integral
sign, and get

$\lim_{\varepsilon\to
0}L_{\varepsilon}^{a}(\phi_{\varepsilon})=\int_{\alpha}^{\beta}2\sqrt{W(r)}\mathop{}\\!\mathrm{d}r\int_{\Sigma_{0}}\frac{1}{|\nabla
h^{a}|}\mathop{}\\!\mathrm{d}S(x).$

Since $\Sigma_{0}$ is just $\partial D$, and since $|\nabla
h^{a}(x)|=|\bm{P}(x;0)|=|a(x,\bm{\nu})|^{-\frac{1}{2}}$ ($x\in\partial D$)
according to (C.3b),

$\lim_{\varepsilon\to
0}L_{\varepsilon}^{a}(\phi_{\varepsilon})=2c_{0}\int_{\partial
D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x),$

matching the form in Lemma 10.

All previous derivations are based on the assumption that $a^{ij}(x)\in
C^{\infty}$. For the general case $a^{ij}\in C^{0}$, we only need to prove the
existence of a subsequence, i.e. $\\{\varepsilon_{k}\\}\to 0$ and
$\phi_{k}\to\phi_{0}$ a.e. such that

$\limsup_{k\to\infty}L_{\varepsilon_{k}}^{a}(\phi_{k})\leq
2c_{0}\int_{\partial
D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x).$

Choose a sequence of smooth matrix functions $\mathbf{A}_{k}(x)\in C^{\infty}$
such that $\mathbf{A}_{k}\to\mathbf{A}$ uniformly and that
$\mathbf{A}_{k}(x)\succeq\mathbf{A}(x),\forall x$. Let
$\\{\phi^{k}_{\varepsilon}\\}$ be the approximation sequence satisfying Lemma
10 with respect to $\mathbf{A}_{k}$. For each $k$, choose
$\varepsilon_{k}<\min\\{\frac{1}{k},\varepsilon_{k-1}\\}$ such that

$L_{\varepsilon}^{a_{k}}(\phi^{k}_{\varepsilon_{k}})\leq 2c_{0}\int_{\partial
D}|a_{k}(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x)+\frac{1}{k},$

and denote by $\phi_{k}=\phi^{k}_{\varepsilon_{k}}$. It is then obvious that

$\limsup_{k\to\infty}L_{\varepsilon_{k}}^{a}(\phi_{k})\leq\limsup_{k\to\infty}2c_{0}\int_{\partial
D}|a_{k}(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x)=2c_{0}\int_{\partial
D}|a(x,\bm{\nu})|^{\frac{1}{2}}\mathop{}\\!\mathrm{d}S(x),$

as desired.

## Acknowledgements

Dawei Wu and Baoming Shi appreciate the University of Strathclyde for their
support and hospitality when work on this paper was undertaken.

## References

  * [1] Robert A. Adams and John J. F. Fournier.  Sobolev spaces.  Number v. 140 in Pure and applied mathematics. Academic Press, Amsterdam Boston, 2nd edition, 2003.
  * [2] James H. Adler, Anca S. Andrei, and Timothy J. Atherton.  Nonlinear Methods for Shape Optimization Problems in Liquid Crystal Tactoids, October 2023.  arXiv:2310.04022 [cs, math].
  * [3] Kenneth A. Brakke.  The Surface Evolver.  Experimental Mathematics, 1(2):141–165, January 1992.
  * [4] John W. Cahn and John E. Hilliard.  Free Energy of a Nonuniform System. I. Interfacial Free Energy.  The Journal of Chemical Physics, 28(2):258–267, February 1958.
  * [5] Giacomo Canevari, Joseph Harris, Apala Majumdar, and Yiwei Wang.  The well order reconstruction solution for three-dimensional wells, in the Landau–de Gennes theory.  International Journal of Non-Linear Mechanics, 119:103342, March 2020.
  * [6] Giacomo Canevari, Apala Majumdar, and Amy Spicer.  Order Reconstruction for Nematics on Squares and Hexagons: A Landau–de Gennes Study.  SIAM Journal on Applied Mathematics, 77(1):267–293, January 2017.
  * [7] Huangxin Chen, Haitao Leng, Dong Wang, and Xiao-Ping Wang.  An efficient threshold dynamics method for topology optimization for fluids.  CSIAM Transactions on Applied Mathematics, 3(1):26–56, June 2022.
  * [8] Shibin Dai, Bo Li, and Jianfeng Lu.  Convergence of Phase-Field Free Energy and Boundary Force for Molecular Solvation.  Archive for Rational Mechanics and Analysis, 227(1):105–147, January 2018.
  * [9] Pierre Gilles de Gennes and Jacques Prost.  The Physics of Liquid Crystals.  Number 83 in Oxford science publications. Clarendon Press; Oxford University Press, Oxford : New York, 1st edition, 1974.
  * [10] Ennio de Giorgi.  Convergence problems for functionals and operators.  In Proceedings of the International Meeting on Recent Methods in Nonlinear Analysis, pages 131–188. Pitagora Ed., Bologna, 1978.
  * [11] Andrew DeBenedictis and Timothy J. Atherton.  Shape minimisation problems in liquid crystals.  Liquid Crystals, 43(13-15):2352–2362, December 2016.
  * [12] M. C. Delfour and J.-P. Zolésio.  Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization, Second Edition.  Society for Industrial and Applied Mathematics, 2nd edition, January 2011.
  * [13] J. W. Doane, A. Golemme, J. L. West, J. B. Whitehead, and B.-G. Wu.  Polymer Dispersed Liquid Crystals for Display Application.  Molecular Crystals and Liquid Crystals Incorporating Nonlinear Optics, 165(1):511–532, December 1988.
  * [14] Lawrence C. Evans.  Partial differential equations.  Number 19 in Graduate studies in mathematics. American Mathematical Society, Providence, Rhode Island, 2nd edition, 2010.
  * [15] Lawrence C. Evans and Ronald F. Gariepy.  Measure theory and fine properties of functions.  Textbooks in mathematics. CRC Press, Boca Raton, Fla., revised edition, 2015.
  * [16] M. Gregory Forest, Qi Wang, and Xiaofeng Yang.  LCP droplet dispersions: a two-phase, diffuse-interface kinetic theory and global droplet defect predictions.  Soft Matter, 8(37):9642, 2012.
  * [17] Harald Garcke, Claudia Hecht, Michael Hinze, and Christian Kahle.  Numerical approximation of phase field based shape and topology optimization for fluids.  SIAM Journal on Scientific Computing, 37(4):A1846–A1871, January 2015.
  * [18] David Gilbarg and Neil S. Trudinger.  Elliptic Partial Differential Equations of Second Order, volume 224 of Grundlehren der mathematischen Wissenschaften.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1977.
  * [19] Enrico Giusti.  Minimal Surfaces and Functions of Bounded Variation.  Birkhäuser Boston, Boston, MA, 1984.
  * [20] Dmitry Golovaty, Peter Sternberg, and Raghavendra Venkatraman.  A Ginzburg–Landau-Type Problem for Highly Anisotropic Nematic Liquid Crystals.  SIAM Journal on Mathematical Analysis, 51(1):276–320, January 2019.
  * [21] Jiequn Han, Yi Luo, Wei Wang, Pingwen Zhang, and Zhifei Zhang.  From microscopic theory to macroscopic theory: A systematic study on modeling for liquid crystals.  Archive for Rational Mechanics and Analysis, 215(3):741–809, March 2015.
  * [22] Yucen Han and Apala Majumdar.  Multistability for a reduced nematic liquid crystal model in the exterior of 2D polygons, October 2022.  arXiv:2112.05511 [math-ph].
  * [23] Yucen Han, Apala Majumdar, and Lei Zhang.  A reduced study for nematic equilibria on two-dimensional polygons.  SIAM Journal on Applied Mathematics, 80(4):1678–1703, January 2020.
  * [24] Yucen Han, Jianyuan Yin, Pingwen Zhang, Apala Majumdar, and Lei Zhang.  Solution landscape of a reduced Landau–de Gennes model on a hexagon.  Nonlinearity, 34(4):2048–2069, April 2021.
  * [25] Duvan Henao and Apala Majumdar.  Symmetry of Uniaxial Global Landau–de Gennes Minimizers in the Theory of Nematic Liquid Crystals.  SIAM Journal on Mathematical Analysis, 44(5):3217–3241, January 2012.
  * [26] Yucheng Hu, Yang Qu, and Pingwen Zhang.  On the disclination lines of nematic liquid crystals.  Communications in Computational Physics, 19(2):354–379, February 2016.
  * [27] Frank Jülicher and Reinhard Lipowsky.  Shape transformations of vesicles with intramembrane domains.  Physical Review E, 53(3):2670–2683, March 1996.
  * [28] A. V. Kaznacheev, M. M. Bogdanov, and A. S. Sonin.  The influence of anchoring energy on the prolate shape of tactoids in lyotropic inorganic liquid crystals.  Journal of Experimental and Theoretical Physics, 97(6):1159–1167, December 2003.
  * [29] Felix C. Keber, Etienne Loiseau, Tim Sanchez, Stephen J. DeCamp, Luca Giomi, Mark J. Bowick, M. Cristina Marchetti, Zvonimir Dogic, and Andreas R. Bausch.  Topology and dynamics of active nematic vesicles.  Science, 345(6201):1135–1139, September 2014.
  * [30] Young-Ki Kim, Sergij V Shiyanovskii, and Oleg D Lavrentovich.  Morphogenesis of defects and tactoids during isotropic–nematic phase transition in self-assembled lyotropic chromonic liquid crystals.  Journal of Physics: Condensed Matter, 25(40):404202, October 2013.
  * [31] Stephen Kitson and Adrian Geisow.  Controllable alignment of nematic liquid crystals around microscopic posts: Stabilization of multiple states.  Applied Physics Letters, 80(19):3635–3637, May 2002.
  * [32] M Kleman.  Defects in liquid crystals.  Reports on Progress in Physics, 52(5):555–654, May 1989.
  * [33] Samo Kralj and Apala Majumdar.  Order reconstruction patterns in nematic liquid crystal wells.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 470(2169):20140276, September 2014.
  * [34] Xavier Lamy.  Bifurcation Analysis in a Frustrated Nematic Cell.  Journal of Nonlinear Science, 24(6):1197–1230, December 2014.
  * [35] Marco Leoni, Oksana V. Manyuhina, Mark J. Bowick, and M. Cristina Marchetti.  Defect driven shapes in nematic droplets: analogies with cell division.  Soft Matter, 13(6):1257–1266, 2017.
  * [36] Jun Li, M.G. Forest, Qi Wang, and R. Zhou.  A kinetic theory and benchmark predictions for polymer-dispersed, semi-flexible macromolecular rods or platelets.  Physica D: Nonlinear Phenomena, 240(2):114–130, January 2011.
  * [37] Fanghua Lin and Changyou Wang.  Isotropic‐nematic phase transition and liquid crystal droplets.  Communications on Pure and Applied Mathematics, 76(9):1728–1792, September 2023.
  * [38] John S. Lowengrub, Andreas Rätz, and Axel Voigt.  Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission.  Physical Review E, 79(3):031926, March 2009.
  * [39] Apala Majumdar.  Equilibrium order parameters of nematic liquid crystals in the Landau-de Gennes theory.  European Journal of Applied Mathematics, 21(2):181–203, April 2010.
  * [40] Yonit Maroudas-Sacks, Liora Garion, Lital Shani-Zerbib, Anton Livshits, Erez Braun, and Kinneret Keren.  Topological defects in the nematic order of actin fibres as organization centres of Hydra morphogenesis.  Nature Physics, 17(2):251–259, February 2021.
  * [41] Matthew Mata, Carlos J. García-Cervera, and Hector D. Ceniceros.  Ordering kinetics of a conserved binary mixture with a nematic liquid crystal component.  Journal of Non-Newtonian Fluid Mechanics, 212:18–27, October 2014.
  * [42] Luciano Modica.  The gradient theory of phase transitions and the minimal interface criterion.  Archive for Rational Mechanics and Analysis, 98(2):123–142, June 1987.
  * [43] Nigel J. Mottram and Christopher J. P. Newton.  Introduction to Q-tensor theory, September 2014.  arXiv:1409.3542 [cond-mat].
  * [44] Yu. A. Nastishin, H. Liu, T. Schneider, V. Nazarenko, R. Vasyuta, S. V. Shiyanovskii, and O. D. Lavrentovich.  Optical characterization of the nematic lyotropic chromonic liquid crystals: Light absorption, birefringence, and scalar order parameter.  Physical Review E, 72(4):041711, October 2005.
  * [45] Stanley Osher and Ronald Fedkiw.  Level Set Methods and Dynamic Implicit Surfaces, volume 153 of Applied Mathematical Sciences.  Springer New York, New York, NY, 2003.
  * [46] Silvia Paparini and Epifanio G. Virga.  Spiralling defect cores in chromonic hedgehogs.  Liquid Crystals, 50(7-10):1498–1516, August 2023.
  * [47] Walter Rudin.  Real and complex analysis.  McGraw-Hill, New York, 3rd ed edition, 1987.
  * [48] Yoshiki Sawa, Kenji Urayama, Toshikazu Takigawa, Vianney Gimenez-Pinto, Badel L. Mbanga, Fangfu Ye, Jonathan V. Selinger, and Robin L. B. Selinger.  Shape and chirality transitions in off-axis twist nematic elastomer ribbons.  Physical Review E, 88(2):022502, August 2013.
  * [49] Baoming Shi, Yucen Han, Apala Majumdar, and Lei Zhang.  Multistability for Nematic Liquid Crystals in Cuboids with Degenerate Planar Boundary Conditions.  SIAM Journal on Applied Mathematics, 84(2):756–781, April 2024.
  * [50] Baoming Shi, Yucen Han, Jianyuan Yin, Apala Majumdar, and Lei Zhang.  Hierarchies of critical points of a Landau-de Gennes free energy on three-dimensional cuboids.  Nonlinearity, 36(5):2631–2654, May 2023.
  * [51] Baoming Shi, Yucen Han, and Lei Zhang.  Nematic Liquid Crystals in a Rectangular Confinement: Solution Landscape, and Bifurcation.  SIAM Journal on Applied Mathematics, 82(5):1808–1828, December 2022.
  * [52] I Singer-Loginova and H M Singer.  The phase field technique for modeling multiphase materials.  Reports on Progress in Physics, 71(10):106501, October 2008.
  * [53] Andrei A. Sonin.  Pierre-Gilles de Gennes and physics of liquid crystals.  Liquid Crystals Reviews, 6(2):109–128, July 2018.
  * [54] T. J. Spencer, C. M. Care, R. M. Amos, and J. C. Jones.  Zenithal bistable device: Comparison of modeling and experiment.  Physical Review E, 82(2):021702, August 2010.
  * [55] Peter Sternberg.  The effect of a singular perturbation on nonconvex variational problems.  Archive for Rational Mechanics and Analysis, 101(3):209–260, September 1988.
  * [56] Abul K. Tarafder, Andriko Von Kügelgen, Adam J. Mellul, Ulrike Schulze, Dirk G. A. L. Aarts, and Tanmay A. M. Bharat.  Phage liquid crystalline droplets form occlusive sheaths that encapsulate and protect infectious rod-shaped bacteria.  Proceedings of the National Academy of Sciences, 117(9):4724–4731, March 2020.
  * [57] Luana Tortora, Heung-Shik Park, Shin-Woong Kang, Victoria Savaryn, Seung-Ho Hong, Konstantine Kaznatcheev, Daniele Finotello, Samuel Sprunt, Satyendra Kumar, and Oleg D. Lavrentovich.  Self-assembly, condensation, and order in aqueous lyotropic chromonic liquid crystals crowded with additives.  Soft Matter, 6(17):4157, 2010.
  * [58] J. D. Van Der Waals.  The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density.  Journal of Statistical Physics, 20(2):200–244, February 1979.
  * [59] Epifanio G. Virga.  Variational Theories for Liquid Crystals.  Number 8 in Applied mathematics and mathematical computation. Chapman & Hall, London Weinheim, 1st edition, 1994.
  * [60] Wei Wang, Lei Zhang, and Pingwen Zhang.  Modelling and computation of liquid crystals.  Acta Numerica, 30:765–851, May 2021.
  * [61] Yunzhi Wang and Ju Li.  Phase field modeling of defects and deformation.  Acta Materialia, 58(4):1212–1235, February 2010.
  * [62] John L. West.  Phase separation of liquid crystals in polymers.  Molecular Crystals and Liquid Crystals Incorporating Nonlinear Optics, 157(1):427–441, April 1988.
  * [63] Xiangjun Xing, Homin Shin, Mark J. Bowick, Zhenwei Yao, Lin Jia, and Min-Hui Li.  Morphology of nematic and smectic vesicles.  Proceedings of the National Academy of Sciences, 109(14):5202–5206, April 2012.
  * [64] Jianyuan Yin, Yiwei Wang, Jeff Z. Y. Chen, Pingwen Zhang, and Lei Zhang.  Construction of a pathway map on a complicated energy landscape.  Physical Review Letters, 124(9):090601, March 2020.  arXiv:1909.13182 [cond-mat, physics:physics].
  * [65] Jianyuan Yin, Lei Zhang, and Pingwen Zhang.  Solution landscape of the Onsager model identifies non-axisymmetric critical points.  Physica D: Nonlinear Phenomena, 430:133081, February 2022.
  * [66] Pengtao Yue, James J. Feng, Chun Liu, and Jie Shen.  A diffuse-interface method for simulating two-phase flows of complex fluids.  Journal of Fluid Mechanics, 515:293–317, September 2004.
  * [67] Dimitrios C. Zografopoulos, Romeo Beccherelli, and Emmanouil E. Kriezis.  Beam-splitter switches based on zenithal bistable liquid-crystal gratings.  Physical Review E, 90(4):042503, October 2014.