remarkRemark hypothesisHypothesis claimClaim assumptionAssumption Optimal
Control of the Landau-de Gennes ModelT. M. Surowiec and S. W. Walker

# Optimal Control of the Landau-de Gennes Model
of Nematic Liquid Crystals††thanks: Submitted to the editors DATE.

Thomas M. Surowiec  Simula Research Laboratory, Department of Numerical
Analysis and Scientific Computing, Kristian Augusts gate 23, 0164 Oslo, Norway
<EMAIL_ADDRESS>Shawn W. Walker  Department of Mathematics, Center for
Computation and Technology (CCT), Louisiana State University, Baton Rouge, LA
70803, Tel.: +1-225-578-1603, Fax: +1-225-578-4276, S. W. Walker acknowledges
financial support by the NSF DMS-2111474<EMAIL_ADDRESS>
###### Abstract

We present an analysis and numerical study of an optimal control problem for
the Landau-de Gennes (LdG) model of nematic liquid crystals (LCs), which is a
crucial component in modern technology. They exhibit long range orientational
order in their nematic phase, which is represented by a tensor-valued
(spatial) order parameter $Q=Q(x)$. Equilibrium LC states correspond to $Q$
functions that (locally) minimize an LdG energy functional. Thus, we consider
an $L^{2}$-gradient flow of the LdG energy that allows for finding local
minimizers and leads to a semi-linear parabolic PDE, for which we develop an
optimal control framework. We then derive several a priori estimates for the
forward problem, including continuity in space-time, that allow us to prove
existence of optimal boundary and external “force” controls and to derive
optimality conditions through the use of an adjoint equation. Next, we present
a simple finite element scheme for the LdG model and a straightforward
optimization algorithm. We illustrate optimization of LC states through
numerical experiments in two and three dimensions that seek to place LC
defects (where $Q(x)=0$) in desired locations, which is desirable in
applications.

###### keywords:

nematic liquid crystals, defects, finite element method, adjoint equation

49M25, 35K91, 65N30

##  1 Introduction

Liquid crystals (LCs) are a critical material for emerging technologies [20,
37]. Their response to optical [8, 28, 53, 34, 18], electric/magnetic [14, 1,
51], and mechanical actuation [63, 21, 7, 48] has already yielded various
devices, e.g. electronic shutters [32], novel types of lasers [35, 16],
dynamic shape control of elastic bodies [15, 58], and others [46, 40, 55, 59,
57]. In fact, [60] demonstrate that LCs can enable logic operations within
soft matter, which can lead to creating autonomous active materials with the
capability to make decisions. Thus, optimization of LC devices in these
applications is of obvious interest.

LCs are considered a _meso-phase_ of matter in which its ordered macroscopic
state is _between_ a spatially disordered liquid, and a fully crystalline
solid [56]. In their nematic phase, in which long ranged orientational order
exists, the Landau-de Gennes (LdG) theory introduces a _tensor-valued_
function $Q$ to describe local order in the LC material. In particular, the
eigenframe of $Q$ yields information about the statistics of the distribution
of LC molecule orientations; see [56, Sec. 1.3] for an excellent derivation.
The energy functional for $Q$, which is minimized at equilibrium, involves
both a bulk potential, of “double-well” type, and an elastic contribution
involving the derivatives of $Q$. Often, an $L^{2}$-type gradient flow is used
to compute (local) minimizers of the LdG energy functional.

The goals of this paper are to formulate an optimal control problem for the
$L^{2}$-gradient flow of the LdG energy, derive several analytical results,
and demonstrate the ability to optimize LC behavior with numerical
simulations. To the best of our knowledge, a fully fledged, PDE-based, optimal
control formulation of the LdG model of LCs has not been done before.
Utilizing both boundary controls and external “force” controls, we prove
existence of optimal controls for the LdG model. In addition, we show several
numerical experiments, of tracking control type, that seek to place LC
_defects_ in desired locations. Defects correspond to sudden spatial changes
in $Q$ and are discussed more thoroughly in Section 7; also see [39] for an
introduction to defects in mathematical models of liquid crystals. Our method
should be useful for optimizing LC devices in a variety of applications.

The paper is organized as follows. Section 2 explains the LdG model and the
associated optimal control problem, as well as discuss related work on the
Allen-Cahn equation. The well-posedness of the parabolic PDE coming from the
$L^{2}$-gradient flow of the LdG energy is established in Section 3 along with
several analytical results. Existence of optimal controls is shown in Section
4 and first order optimality conditions are established in Section 5. Section
6 describes our finite element method for approximating the forward and
adjoint problems; see [2, 19, 47, 5, 9, 38, 61] for other numerical methods
for models related to LdG. We illustrate our method with numerical experiments
in Section 7 and close with some remarks in Section 8.

##  2 Liquid Crystal Theory

This section reviews the Landau-de Gennes (LdG) theory for a nematic phase. We
start with the following clarifications.

First, we note that the minimization of the standard free energy of the LdG
model gives rise to a semilinear elliptic partial differential equation (PDE),
which admits multiple solutions. Theoretically, this semilinear equation could
be included as a constraint in an optimization problem. However, there is no
guarantee that the second derivative of the LdG free energy functional is
surjective. This would significantly complicate the derivation of optimality
conditions and severely limit the convergence theory of numerical optimization
algorithms. To remedy these issues, we will consider an evolution equation,
which amounts to an $L^{2}$-gradient flow of the LdG free energy. This time-
dependent control strategy is analyzed in subsequent sections.

The second point of clarification involves the bulk potential used to model
the nematic-to-isotropic phase transition $\widetilde{\psi}$. In the
discussion below, we will first introduce a traditional double well function
and derive an associated evolution equation. For mathematical reasons, we then
modify this term beyond physically meaningful values of $|Q|$.

##### Notation

We typically denote scalars and vectors with lowercase letters, while tensors
are denoted with uppercase letters. Boldface capital letters typically denote
vector spaces or function spaces. Standard notation is used for Sobolev spaces
and inner products.

###  2.1 Landau-de Gennes Model

Let $\mathbf{S}$ be the space of symmetric, $d\times d$ tensors, and
$\mathbf{S}_{0}$ the set of symmetric, traceless $d\times d$ tensors, where
$d=2$ or $3$. The order parameter of the LdG theory is given by
$Q\in\mathbf{S}_{0}$, which represents the statistical distribution (i.e. a
covariance matrix) of LC molecules at a given point in space [56]. This means
that the eigenvalues $\lambda_{i}(Q)$ should satisfy the following bound:
$-1/d\leq\lambda_{i}(Q)\leq(d-1)/d$ for $i=1,...,d$. In the standard model
discussed below, the eigenvalue bounds are not explicitly enforced, though
they are usually satisfied through the effect of the double-well (c.f. [41]).

We mainly focus on the $d=3$ case, and represent the state of the LC material
by a tensor-valued function $Q:\Omega\to\mathbf{S}_{0}$, where
$\Omega\subset\mathbb{R}^{3}$ is the physical domain of interest. Moreover, we
take $\Omega$ to be an open, bounded, Lipschitz domain with boundary $\Gamma$,
and outward unit normal $\bm{\nu}$; the normal derivative is denoted
$\partial_{\bm{\nu}}$. The standard free energy of the LdG model is defined as
[43, 44]:

(1) $\begin{split}\mathcal{E}[Q]&:=\int_{\Omega}W(Q,\nabla
Q)\,dx+\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\widetilde{\psi}(Q)\,dx+\eta_{\Gamma}\int_{\Gamma}f_{\Gamma}(Q)\,dS(x)-\int_{\Omega}\chi(Q)\,dx,\end{split}$

where $W(Q,\nabla Q)$ is the elastic energy density [43, 44]. Since optimal
control of the Landau-de Gennes model has yet to be developed, we simplify
$W(Q,\nabla Q)$ to the one-constant model, i.e. $W(Q,\nabla
Q)=\frac{1}{2}|\nabla Q|^{2}$. Future extensions of this work will consider
more general elastic energies.

The bulk potential $\widetilde{\psi}$ models the nematic-to-isotropic phase
transition. It is a (non-symmetric) double-well type of function that is
typically given by

(2)
$\begin{split}\widetilde{\psi}(Q)=a_{0}-\frac{a_{2}}{2}\mathrm{tr}(Q^{2})-\frac{a_{3}}{3}\mathrm{tr}(Q^{3})+\frac{a_{4}}{4}\left(\mathrm{tr}(Q^{2})\right)^{2},\quad\widetilde{\psi}\geq
0,\end{split}$

where $a_{2}$, $a_{3}$, $a_{4}$ are material parameters. The choice of
constants affects the stability of the nematic phase. Since we are only
interested in the nematic phase, $a_{2}$, $a_{3}$, $a_{4}$ are positive, and
$a_{0}$ is a convenient positive constant to ensure $\widetilde{\psi}\geq 0$.
Stationary points of $\widetilde{\psi}$ are either _uniaxial_ or _isotropic_
[41], i.e. if $Q$ is uniaxial, then it corresponds to

(3) $\begin{split}Q_{ij}=s_{0}\left(n_{i}n_{j}-\delta_{ij}/3\right),\text{ for
}1\leq i,j\leq 3,\end{split}$

where $s_{0}>0$ depends on $\widetilde{\psi}$, $[n_{i}]_{i=1}^{3}\equiv
n\in\mathbb{R}^{3}$ is a unit vector, and $\delta_{ij}$ is the Kronecker
delta; $Q=0$ is the isotropic state. The parameter $\eta_{\mathrm{dw}}>0$,
appearing in Eq. 1, is known as the _nematic correlation length_ , and usually
satisfies $\eta_{\mathrm{dw}}\ll 1$.

The surface energy $f_{\Gamma}(Q)$, with parameter $\eta_{\Gamma}>0$, accounts
for _weak anchoring_ of the LC material at the boundary, i.e. it imposes an
energetic penalty on the boundary conditions for $Q$. In this paper, we use a
Rapini-Papoular type anchoring energy [4]:

(4)
$\begin{split}f_{\Gamma}(Q)=\frac{1}{2}\mathrm{tr}\left(Q-U_{\Gamma}\right)^{2}\equiv\frac{1}{2}|Q-U_{\Gamma}|^{2},\end{split}$

where $U_{\Gamma}:\Gamma\to\mathbf{S}_{0}$ and we take $U_{\Gamma}$ to be one
of the control variables in the optimal control problem stated in Section 2.2.
The function $\chi(\cdot)$ is used to (approximately) model interactions of
the LC material with external fields, e.g. an electric field. In this paper,
we take $\chi(Q)=Q:U_{\Omega}$, where $U_{\Omega}:\Omega\to\mathbf{S}_{0}$ is
also a control variable.

Local minimizers of the energy $\mathcal{E}[Q]$ can be found through an
$L^{2}$ gradient flow, which can be thought of as a simple damped,
evolutionary LdG model. This leads to the following parabolic equation for $Q$
in strong form:

(5a) $\displaystyle Q_{t}-\Delta
Q+\frac{1}{\eta_{\mathrm{dw}}^{2}}\widehat{\widetilde{\psi}^{\prime}(Q)}$
$\displaystyle=\eta_{\Omega}U_{\Omega},\text{ in
}\Omega\times(0,t_{\mathrm{f}}),$ (5b)
$\displaystyle\partial_{\bm{\nu}}Q+\eta_{\Gamma}Q$
$\displaystyle=\eta_{\Gamma}U_{\Gamma},\text{ on
}\Gamma\times(0,t_{\mathrm{f}}),$ (5c) $\displaystyle Q(\cdot,0)$
$\displaystyle=Q_{0},\quad~{}\text{ in }\Omega,$

where $\widehat{M}$ denotes the traceless part of a symmetric tensor $M$,
$Q_{0}:\Omega\to\mathbf{S}_{0}$ is the initial condition, and
$t_{\mathrm{f}}>0$ is a given final time. The system Eq. 5 can be viewed as a
tensor-valued analog of the Allen-Cahn equation with Robin boundary
conditions. Taking $\eta_{\Gamma}\to\infty$ recovers _strong anchoring_ , i.e.
the Dirichlet condition $Q=U_{\Gamma}$ on $\Gamma$. In this paper,
$\eta_{\Gamma}>0$ is fixed and finite.

The first and second derivatives of $\widetilde{\psi}$ are a 2-tensor and
4-tensor, respectively:

(6)
$\begin{split}\widetilde{\psi}^{\prime}(Q)&=-a_{2}Q-a_{3}Q^{2}+a_{4}\,\mathrm{tr}(Q^{2})Q,\\\
[\widetilde{\psi}^{\prime\prime}(Q)]_{ijkl}&=-a_{2}\delta_{ik}\delta_{jl}-2a_{3}Q_{jk}\delta_{li}+a_{4}\left(2Q_{ij}Q_{kl}+\mathrm{tr}(Q^{2})\delta_{ik}\delta_{jl}\right),\end{split}$

where the second derivative is written with indices for clarity. In addition,
for analytical purposes, we modify the bulk potential to have quadratic growth
as $|Q|\to\infty$. For instance, let $\rho:[0,\infty)\to\mathbb{R}_{+}$ be a
$C^{\infty}$ cut-off function such that

$\begin{split}\rho(r)=1,\text{ if
}r<b_{1},~{}~{}\rho(r)=\text{monotone},\text{ if }b_{1}\leq r\leq
b_{2},~{}~{}\rho(r)=0,\text{ if }r>b_{2},\end{split}$

where $1\leq b_{1}<b_{2}$ are two fixed constants. Then, the modified
potential is given by

(7)
$\begin{split}\psi(Q)=\widetilde{\psi}(Q)\rho(\mathrm{tr}(Q^{2}))+a_{4}^{2}\mathrm{tr}(Q^{2})\left[1-\rho(\mathrm{tr}(Q^{2}))\right],\end{split}$

which will be used throughout the remainder of the paper. Clearly, there exist
uniform constants $c_{0}$, $c_{1}$, $c_{2}$, depending on $a_{0}$, $a_{2}$,
$a_{3}$, $a_{4}$, $b_{1}$, $b_{2}$, such that

(8) $\begin{split}|\psi(Q)|\leq a_{0}+c_{0}|Q|^{2},\quad|\psi^{\prime}(Q)|\leq
c_{1}|Q|,\quad|\psi^{\prime\prime}(Q)|\leq c_{2},~{}\text{ for all
}Q\in\mathbf{S}_{0},\end{split}$

noting that $|Q|^{2}\equiv\mathrm{tr}(Q^{2})$. For typical choices of the
physical constants in Eq. 2, choosing $b_{1}=1$ and $b_{2}=2$ is effective
since this modification does not change the location of the local minimizers.
In addition to Eq. 8 we observe that $\psi^{\prime\prime}$ is globally
Lipschitz with uniformly bounded derivative $\psi^{\prime\prime\prime}(Q)$ for
all $Q\in\mathbf{S}_{0}$.

In addition, we will often make use of the convex splitting

(9)
$\begin{split}\psi(Q)&\equiv\left[(d_{2}/2)\mathrm{tr}(Q^{2})+\psi(Q)\right]-(d_{2}/2)\mathrm{tr}(Q^{2})=:\psi_{c}(Q)-\psi_{e}(Q),\end{split}$

where $\psi_{c},\psi_{e}$ are non-negative convex functions with $d_{2}>0$
(i.e. a “stabilization” constant) chosen sufficiently large to ensure a convex
split. In particular, we note that $\psi_{c}^{\prime}(Q)$ and
$\psi_{e}^{\prime}(Q)$ are monotone functions and there is a constant
$0<\tilde{d}_{2}<d_{2}$ such that, for $d_{2}$ sufficiently large,
$\psi_{c}(Q)$ satisfies the lower bounds:

(10) $\begin{split}\psi_{c}(Q)\geq
a_{0}+\frac{\tilde{d}_{2}}{2}|Q|^{2},\quad\psi_{c}^{\prime}(Q):Q\geq\tilde{d}_{2}|Q|^{2},\quad
P:\psi_{c}^{\prime\prime}(Q):P\geq 3\tilde{d}_{2}|P|^{2}.\end{split}$

Henceforth, we take all constants to be non-dimensional; see [24] for a
detailed treatment of how the LdG model is non-dimensionalized.

###  2.2 Optimal control problem

We formulate an optimal tracking control problem for the LdG model. The
following Sobolev spaces are used throughout:

$\begin{split}\mathbf{V}:=H^{1}(\Omega;\mathbf{S}_{0}),\quad\mathbf{H}:=L^{2}(\Omega;\mathbf{S}_{0}),\quad\mathbf{H}_{\Gamma}:=L^{2}(\Gamma;\mathbf{S}_{0}),\end{split}$

where each space is endowed with its respective natural norm. For every
$t\in(0,t_{\mathrm{f}}]$, the space-time cylinder and boundary are defined as

(11)
$\begin{split}\mathcal{C}_{t}&:=\Omega\times(0,t),~{}~{}\mathcal{C}\equiv\mathcal{C}_{t_{\mathrm{f}}},~{}~{}\mathcal{G}_{t}:=\Gamma\times(0,t),\text{
and }\mathcal{G}\equiv\mathcal{G}_{t_{\mathrm{f}}}.\end{split}$

Next, we introduce target functions:

(12) $\begin{split}Z_{\mathcal{C}}&\in L^{2}(\mathcal{C}),\quad
Z_{\mathcal{G}}\in L^{2}(\mathcal{G}),\quad Z_{t_{\mathrm{f}}}\in
H^{1}(\Omega).\end{split}$

In contrast to optimal control problems with scalar or vector-valued controls,
the bound constraints used to defined the set of admissible controls
$\mathcal{U}_{\rm{ad}}$ are slightly more complicated. These are discussed
below.

We now define the optimal control problem: minimize the functional

(13)
$\begin{split}J(Q,U_{\Omega},U_{\Gamma})&:=\frac{\beta_{\mathcal{C}}}{2}\|Q-Z_{\mathcal{C}}\|^{2}_{L^{2}(\mathcal{C})}+\frac{\beta_{\mathcal{G}}}{2}\|Q-Z_{\mathcal{C}}\|^{2}_{L^{2}(\mathcal{G})}\\\
+\frac{\beta_{t_{\mathrm{f}}}}{2}&\|Q(\cdot,t_{\mathrm{f}})-Z_{t_{\mathrm{f}}}\|^{2}_{\mathbf{H}}+\frac{\alpha_{\mathcal{C}}}{2}\|U_{\Omega}\|^{2}_{L^{2}(\mathcal{C})}+\frac{\alpha_{\mathcal{G}}}{2}\|U_{\Gamma}\|^{2}_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})},\end{split}$

over a set of admissible controls
$(U_{\Omega},U_{\Gamma})\in\mathcal{U}_{\rm{ad}}\subset
L^{2}(\mathcal{C})\times L^{2}(\mathcal{G})$, subject to $Q$ satisfying the
PDE constraint Eq. 5. The coefficients satisfy
$\beta_{\mathcal{C}},\beta_{\mathcal{G}},\beta_{t_{\mathrm{f}}},\alpha_{\mathcal{C}},\alpha_{\mathcal{G}}\geq
0$, where at least one of
$\beta_{\mathcal{C}},\beta_{\mathcal{G}},\beta_{t_{\mathrm{f}}}$ is nonzero,
and $\alpha_{\mathcal{C}},\alpha_{\mathcal{G}}>0$. Tracking objectives as in
Eq. 13 are ubiquitous in optimal control. For our application, the final two
summands are quadratic cost functionals that force a certain regularity. The
first two terms represent the desire to track a transition of targeted
textures $Z_{\mathcal{C}}$, whereas the third summand $Z_{t_{\mathrm{f}}}$ is
a desired stationary texture. Nematic textures correspond to director
orientations associated with $Q$ at each point $x$ in the domain. Therefore, a
desired texture $Z_{t_{\mathrm{f}}}$ could be one in which all directors are
oriented in the same direction, e.g., parallel to a surface, or in which the
directors follow a particular pattern.

In practice, the available control mechanisms may be technically limited,
e.g., finite dimensional, stationary or only on the boundary. We include
distributed controls in the bulk and boundary here for more generality.
Restrictions to the cases just mentioned would not change the core of the
analysis. The control set $\mathcal{U}_{\rm{ad}}$ is always taken to be
nonempty, closed, and convex. An example of such a set is

(14) $\begin{split}\mathcal{U}_{\rm{ad}}:=&\left\\{P\in
L^{2}(\mathcal{C})\mid|P|\leq{u}_{\Omega},\text{ a.e. in
}\mathcal{C}\right\\}\times\left\\{P\in
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})\mid|P|\leq 1,\text{ a.e. in
}\mathcal{G}\right\\}.\end{split}$

Here, ${u}_{\Omega}$ may be an arbitrary, essentially bounded scalar-valued
function on $\mathcal{C}$. The constant bound $1$ for the boundary controls is
in fact dictated by the application (recall the eigenvalue bounds discussed in
Section 2.1). Note that if $U_{\Gamma}$ is constant in time, then the
$H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$ space is replaced by
$\mathbf{H}_{\Gamma}$, and the $\alpha_{\mathcal{G}}$ term in Eq. 13 becomes
an $\mathbf{H}_{\Gamma}$ norm. In most applications, boundary controls are
_constant_ in time. Allowing the boundary control to vary in time is mainly
for the sake of generality.

Due to the similarities of Eq. 5 with the Allen-Cahn equation, there are a
wide array of relevant contributions in the literature, where optimal control
of Allen-Cahn and related equations, e.g., Cahn-Hilliard, have been studied.
We highlight here several early studies [33, 30, 31], which focused on the
optimal control of Cahn-Hilliard (phase field problems of Caginalp-type); more
recent work [26, 27], in which the author studied the optimal control of
scalar- and vector-valued Allen-Cahn equations with a nonsmooth bulk energy
term (obstacle potential), and [17]. In some sense, [17] is the most relevant.
However, there are several major differences. Our boundary condition has no
diffusive term, because it is not clear how that would manifest in an LC
system, which thus affects our solution’s regularity. Moreover, we are dealing
with a parabolic _system_ with _tensor-valued_ solutions and controls; the PDE
in [17] is scalar-valued. This affects several arguments needed to derive
first-order optimality conditions and greatly increases the difficulty for
numerical methods.

##  3 Well-posedness of the forward problem

To prove well-posedness of Eq. 5, we start with the usual arguments (cf.
[25]). The minimal regularity of the data is given by

(15) $U_{\Omega}\in L^{2}(0,t_{\mathrm{f}};\mathbf{H}),\quad U_{\Gamma}\in
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})\quad Q_{0}\in\mathbf{V},$

and the space of weak solutions we consider is

(16) $\mathcal{W}:=L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})\cap
L^{\infty}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})\cap
H^{1}(0,t_{\mathrm{f}};\mathbf{H})\cap L^{2}(0,t_{\mathrm{f}};\mathbf{V}).$

######  Remark 3.1.

Since $t_{\mathrm{f}}$ is finite and $L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})$
is continuously embedded into $L^{2}(0,t_{\mathrm{f}};\mathbf{V})$, the fourth
space in the definition of $\mathcal{W}$ is redundant. However, we keep it as
written to emphasize that the $L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})$ and
$L^{2}(0,t_{\mathrm{f}};\mathbf{V})$ norms are utilized at different points in
the analysis.

Our notion of weak solution is as follows. We say $Q\in\mathcal{W}$ is a weak
solution of Eq. 5, if $Q(0)\equiv Q(\cdot,t)|_{t=0}=Q_{0}$, and for a.e.
$s\in(0,t_{\mathrm{f}})$, we have

(17) $\begin{split}\left(Q_{t}(s),P(s)\right)_{\mathbf{H}}+&\left(\nabla
Q(s),\nabla
P(s)\right)_{\mathbf{H}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\left(\psi^{\prime}(Q(s)),P(s)\right)_{\mathbf{H}}\\\
+\eta_{\Gamma}\left(Q(s),P(s)\right)_{\mathbf{H}_{\Gamma}}&=\eta_{\Gamma}\left(U_{\Gamma}(s),P(s)\right)_{\mathbf{H}_{\Gamma}}+\eta_{\Omega}\left(U_{\Omega}(s),P(s)\right)_{\mathbf{H}},\end{split}$

for all $P\in H^{1}(0,t_{\mathrm{f}};\mathbf{H})\cap
L^{2}(0,t_{\mathrm{f}};\mathbf{V})$ with $P(0)=0$, where we introduced the
inner products on $\mathbf{H}$ and $\mathbf{H}_{\Gamma}$, respectively.

The solutions of the forward problem Eq. 17 are tensor-valued in space. There
is little work on such problems in the control literature. Nevertheless, in
many instances, we can exploit the Hilbert space structure on $\mathbf{V}$ or
$\mathbf{H}$ and extend the derivations of typical energy estimates and
Lipschitz continuity results. As a consequence, the proofs of any results that
follow the corresponding scalar or vector-valued cases without major changes
have been drastically shortened and placed in the appendix.

###  3.1 Uniqueness of the state and a Lipschitz bound

Under the assumption that solutions $Q$ with the appropriate regularity exist,
we can prove Lipschitz continuity with respect to the input controls and
therefore, a fortiori, uniqueness of solutions. Existence of solutions
ultimately follows from a standard Galerkin approach.

######  Theorem 3.2 (Continuous dependence on the data).

Let $Q_{1},Q_{2}\in\mathcal{W}$ be two solutions of Eq. 17 corresponding to
the input variables $U_{\Omega,i}$, $U_{\Gamma,i}$, $Q_{0,i}$, for $i=1,2$,
which satisfy Eq. 15. Then there exists a constant $c>0$, independent of the
input variables, such that

(18)
$\|Q_{1}-Q_{2}\|^{2}_{C([0,t_{\mathrm{f}}];\mathbf{H})}+\|Q_{1}-Q_{2}\|^{2}_{L^{2}(0,t_{\mathrm{f}};\mathbf{V})}\\\
\leq
c\left(\|Q_{0,1}-Q_{0,2}\|^{2}_{\mathbf{V}}+\eta_{\Omega}\|U_{\Omega,1}-U_{\Omega,2}\|^{2}_{L^{2}(\mathcal{C})}+\eta_{\Gamma}\|U_{\Gamma,1}-U_{\Gamma,2}\|^{2}_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}\right).$

######  Remark 3.3.

We use $c>0$ as a generic constant throughout the text. We also note that the
state space $\mathcal{W}$ can be compactly embedded into the space
$C([0,t_{\mathrm{f}}];\mathbf{H})$.

######  Proof 3.4.

The result follows by standard energy techniques, i.e. first test with the
difference of solutions. The lack of a monotone nonlinear operator is handled
using the convex splitting (9), which exploits the linear growth of
$\psi^{\prime}_{e}$.  Afterwards, we apply weighted Young’s inequalities and
the classic Gronwall lemma. We omit the details.

###  3.2 Existence and Energy Estimates

This section is concerned with the existence of weak solutions. We use a
Faedo-Galerkin approach, for which we require the following assumption. This
condition will be tacitly assumed throughout the remainder of the text.

For each $n\in\mathbb{N}$ (sufficiently large) there is an $n$-dimensional
subspace $\mathbf{V}_{n}$ of $\mathbf{V}$ such that if $\\{Y_{k}\\}_{k=1}^{n}$
is the basis of $\mathbf{V}_{n}$ and $\Pi_{n}:\mathbf{V}\to\mathbf{V}_{n}$ is
the linear projection onto $\mathbf{V}_{n}$, then $\Pi_{n}$ satisfies the
following convergence property: $\|\Pi_{n}P-P\|_{\mathbf{V}}\to 0$, as
$n\to\infty$, for all $P\in\mathbf{V}$.

Typically $\mathbf{V}_{n}$ is based on the eigenvectors of the Laplacian
associated to (in this case) homogeneous Robin boundary conditions. As another
example, when $\Omega$ has a piecewise $C^{2}$ boundary, $\mathbf{V}_{n}$ can
be a conforming finite element space with $n$ nodal degrees-of-freedom defined
over a conforming (curvilinear) mesh of $\Omega$. Next, we define
$Q^{n}_{0}:=\Pi_{n}Q_{0},$ which under Section 3.2 means $Q^{n}_{0}$ converges
strongly to $Q_{0}$ in $V$, as $n\to+\infty$. Set
$Q^{n}_{0,k}:=(Q_{0},Y_{k})_{\mathbf{H}}$ for each $k=1,\dots,n$.

####  3.2.1 Existence of a discrete solution in
$H^{1}(0,t_{\mathrm{f}};\mathbf{V}_{n})$

We start with the existence of unique solutions to the semi-discrete system.

######  Proposition 3.5.

There exists a unique solution $Q^{n}\in
H^{1}(0,t_{\mathrm{f}};\mathbf{V}_{n})$ such that

(19) $\displaystyle\left(Q^{n}_{t}(s),P\right)_{\mathbf{H}}+\left(\nabla
Q^{n}(s),\nabla
P\right)_{\mathbf{H}}+\eta_{\Gamma}\left(Q^{n}(s),P\right)_{\mathbf{H}_{\Gamma}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\psi^{\prime}(Q^{n}(s)):P$
$\displaystyle=\eta_{\Omega}\left(U_{\Omega}(s),P\right)_{\mathbf{H}}+\eta_{\Gamma}\left(U_{\Gamma}(s),P\right)_{\mathbf{H}_{\Gamma}},$

for all $P\in\mathbf{V}_{n}$ and for a.e. $s\in(0,t_{\mathrm{f}})$, with
$Q^{n}(0)=Q^{n}_{0}$.

######  Proof 3.6.

The proof involves a standard application of the Carathéodory existence
theorem, i.e. Eq. 19 reduces to a system of coupled ODEs, for which existence
and uniqueness is straightforward to show.

####  3.2.2 A priori estimates

Given the existence of finite dimensional solutions, we now consider energy
estimates. Let $M_{1}:\mathbb{R}^{3}\to\mathbb{R}$ be given by
$M_{1}(x_{1},x_{2},x_{3})=c(x_{1}+\eta_{\Omega}x_{2}+\eta_{\Gamma}x_{3})$. For
readability, we will often leave off the arguments of $M_{1}$ when it is clear
in context; $c>0$ is a generic constant that can be updated as needed. It will
never depend on $n$, the controls, or the input data.

######  Proposition 3.7.

Suppose that $Q_{0},U_{\Omega},U_{\Gamma}$ satisfy Eq. 15. Then for all
$n\in\mathbb{N}$, the solutions $Q^{n}$ from Proposition 3.5 satisfy the bound

(20) $\|Q^{n}\|^{2}_{L^{2}(0,t_{\mathrm{f}};\mathbf{V})}\leq
M_{1}(\|Q_{0}\|^{2}_{\mathbf{H}},\|U_{\Omega}\|^{2}_{L^{2}(\mathcal{C})},\|U_{\Gamma}\|^{2}_{L^{2}(\mathcal{G})}).$

######  Proof 3.8.

The result follows by similar energy techniques as in the proof of Theorem
3.2; we omit the details.

By exploiting the Hilbert space structure and the nature of the finite
dimensional inner products, we can again use standard derivation techniques to
derive further energy estimates and ultimately prove that the sequence of
finite dimensional solutions is bounded in $\mathcal{W}$. To indicate the
dependence of the bound on the controls, we define
$M_{2}:\mathbb{R}^{4}\to\mathbb{R}$ by
$M_{2}(x_{1},x_{2},x_{3},x_{4})=(1/2)\left(\delta_{1}x_{1}+\delta_{2}x_{2}+\delta_{3}x_{3}+\delta_{4}x_{4}\right)$.
The positive constants $\delta_{1},\dots,\delta_{4}$ are arbitrary and can be
adjusted as needed. Given input controls, we leave off the arguments and
abbreviate both $M_{1}$ and $M_{2}$ by setting

$\displaystyle M_{1}$ $\displaystyle\equiv
M_{1}(\|Q_{0}\|^{2}_{\mathbf{H}},\|U_{\Omega}\|^{2}_{L^{2}(\mathcal{C})},\|U_{\Gamma}\|^{2}_{L^{2}(\mathcal{G})}),$
$\displaystyle M_{2}$ $\displaystyle\equiv
M_{2}(\|U_{\Gamma}(t)\|^{2}_{\mathbf{H}_{\Gamma}},\|U_{\Gamma}(0)\|^{2}_{\mathbf{H}_{\Gamma}},\|(U_{\Gamma})_{t}\|^{2}_{L^{2}(\mathcal{G}_{t})},\|U_{\Omega}\|^{2}_{L^{2}(\mathcal{C}_{t})}).$

######  Proposition 3.9.

Suppose that $Q_{0},U_{\Omega},U_{\Gamma}$ satisfy Eq. 15. Then there exists
an $M_{0}\geq 0$ for all $n$ such that

(21) $\|\nabla
Q^{n}_{0}\|^{2}_{\mathbf{H}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\psi_{c}(Q^{n}_{0})-\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\psi_{e}(Q^{n}_{0})+\delta^{\prime}_{2}\|Q^{n}_{0}\|^{2}_{\mathbf{H}_{\Gamma}}\leq
M_{0},$

holds. Up to rescaling by a generic constant, it also holds for all $n$ a.e.
in $t$ that

(22) $\|Q^{n}_{t}\|^{2}_{L^{2}(\mathcal{C}_{t})}+\|\nabla
Q^{n}(t)\|^{2}_{\mathbf{H}}+\|Q^{n}(t)\|^{2}_{\mathbf{H}_{\Gamma}}\leq
M_{0}+M_{1}+M_{2}.$

Furthermore, up to rescaling by a generic constant, it holds for all $n$ a.e.
in $t$ that

(23) $\|Q^{n}_{t}\|^{2}_{L^{2}(\mathcal{C}_{t})}+\|\nabla
Q^{n}(t)\|^{2}_{\mathbf{H}}+\|Q^{n}(t)\|^{2}_{\mathbf{H}}\leq
M_{0}+M_{1}+M_{2}.$

Finally, as a consequence of Eq. 22, Eq. 23, and Eq. 20, the sequence of
solutions $\left\\{Q^{n}\right\\}$, with $Q^{n}$ from Proposition 3.5, is
bounded in $\mathcal{W}$.

######  Proof 3.10.

See Appendix A.

In order to obtain further properties of the control-to-state mapping, we need
a stronger Lipschitz continuity result which we first state for the semi-
discrete problem.

######  Proposition 3.11.

Suppose that $Q_{0,i},U_{\Omega,i},U_{\Gamma,i}$, for $i=1,2$, satisfy Eq. 15
and let $Q^{n}_{0,i}=\Pi Q_{0,i}$, where $\Pi$ is given in Section 3.2. Then
the corresponding solutions $Q^{n}_{i}$, for $i=1,2$, satisfy the bound

(24) $\|Q^{n}_{1}-Q^{n}_{2}\|_{\mathcal{W}}\leq\\\
c\Big{(}\|U_{\Omega,1}-U_{\Omega,2}\|_{L^{2}(0,t_{\mathrm{f}};\mathbf{H})}^{2}+\|U_{\Gamma,1}-U_{\Gamma,2}\|_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}^{2}+\|Q^{n}_{0,1}-Q^{n}_{0,2}\|^{2}_{\mathbf{V}}\Big{)}^{1/2},$

where $c>0$ is a generic constant that does not depend on the controls,
states, or $n$.

######  Proof 3.12.

The proof is straightforward and directly mirrors the proof of Theorem 3.2. In
particular, the Lipschitz continuity of the gradient of the bulk energy term
is essential. We omit the details.

####  3.2.3 Passage to the limit

In light of the uniform bounds and energy estimates on
$\left\\{Q^{n}\right\\}$ provided above, we can now prove the existence of a
solution.

######  Theorem 3.13.

For every triple $(Q_{0},U_{\Omega},U_{\Gamma})$ that satisfies Eq. 15, there
exists a unique solution $\bar{Q}\in\mathcal{W}$ of the weak form Eq. 17 and a
$\rho>0$ such that

(25) $\|\bar{Q}\|_{\mathcal{W}}\leq\rho.$

######  Proof 3.14.

See Appendix B.

######  Remark 3.15.

Note that despite the fact it depends on $U_{\Omega}$ and $U_{\Gamma}$, the
constant $\rho$ in Eq. 25 can be bounded from above by a uniform constant,
which is independent of $U_{\Omega}$ and $U_{\Gamma}$, provided the latter two
are taken over a bounded set in the space $L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$. This is a direct consequence of
the structure of the a priori estimates.

The previous results also allow us to pass to the limit along a subsequence to
obtain the following global Lipschitz bound from Eq. 24:

(26) $\|Q_{1}-Q_{2}\|_{\mathcal{W}}\leq\\\
c\Big{(}\|U_{\Omega,1}-U_{\Omega,2}\|_{L^{2}(0,t_{\mathrm{f}};\mathbf{H})}^{2}+\|U_{\Gamma,1}-U_{\Gamma,2}\|_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}^{2}+\|Q_{0,1}-Q_{0,2}\|^{2}_{\mathbf{V}}\Big{)}^{1/2},$

provided the controls are feasible.

The final result in this section involves the continuity of the states on the
full space-time cylinder. In contrast to the results above, our nonlinear
system in tensor-valued variables inhibits a direct application of the
standard techniques as can be found, e.g., in the relevant chapters in [54].
We require a few additional steps, which we provide here. The remainder can be
found in the appendix. Note that the following argument is unique to the one
elastic constant case. For more general elastic energy densities, we require
new techniques to derive such continuity results in future studies.

To begin, given the existence of a solution $\bar{Q}$ in $\mathcal{W}$, we
have $\bar{Q}\in L^{6}(\mathcal{C})$, which follows from the Sobolev embedding
theorem and the fact that $t_{\mathrm{f}}<+\infty$. Therefore, $\bar{Q}$ is
the unique solution of the system of linear parabolic equations given by

(27)
$\begin{split}\widehat{Q}_{t}-\Delta\widehat{Q}&=\eta_{\Omega}U_{\Omega}-\frac{1}{\eta_{\mathrm{dw}}^{2}}\psi^{\prime}(\bar{Q}),\text{
in
}\Omega,\quad\partial_{\bm{\nu}}\widehat{Q}+\eta_{\Gamma}\widehat{Q}=\eta_{\Gamma}U_{\Gamma},\text{
on }\Gamma,\end{split}$

with $\widehat{Q}(\cdot,0)=Q_{0}$ in $\Omega$. Next, we use the fact that
there exists a set of five symmetric, traceless, $3\times 3$ orthonormal
matrices $\\{E^{i}\\}_{i=1}^{5}\subset\mathbb{R}^{3\times 3}$ such that every
$Q\in\mathbf{V}$ admits the representation $Q=q_{i}E^{i}$ where $q_{i}\in
H^{1}(\Omega;\mathbb{R})$ for $i=1,\dots,5$. Let $\bar{q}_{i}$ denote the
scalar-valued functions for the solution $\bar{Q}$.

This decomposition was also exploited in [19]. It provides us with an
isometric isomorphism between $\mathbf{V}$ and $H^{1}(\Omega;\mathbb{R}^{5})$
and allows us to split the tensor-based problem into five separate scalar
parabolic equations with Robin boundary conditions. In addition, we note that
the second bound in Eq. 8 implies that $\psi^{\prime}(\bar{Q})$ is also in
$L^{6}(\mathcal{C})$.

Now that we can separate the system into independent scalar-valued equations,
and apply the well-known regularity theory to obtain continuity of $\bar{Q}$.
To be clear, we obtain continuity of each $\bar{q}_{i}$ via e.g., [54, Thm.
5.5], and consequently of $\bar{Q}$. The remainder of the proof is concerned
with removing the dependency on $\bar{Q}$ from the upper bound. Since this
does not require any special techniques for tensor-valued solutions, we have
placed it in Appendix C.

######  Theorem 3.16.

If, in addition to Eq. 15, we have $U_{\Omega}\in L^{r}(\mathcal{C})$,
$U_{\Gamma}\in L^{s}(\mathcal{G})$, $Q_{0}\in C(\overline{\Omega})$, with
$r\in(5/2,6]$ and $s>4$, then $\bar{Q}\in C(\overline{\mathcal{C}})$. Next,
let

(28) $\mathcal{U}_{\rm{reg}}:=[L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})]\cap[L^{r}(\mathcal{C})\times
L^{s}(\mathcal{G})],$

endowed with the natural norm

(29)
$\|(U_{\Omega},U_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}}=\max\\{\|(U_{\Omega},U_{\Gamma})\|_{L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})},\|(U_{\Omega},U_{\Gamma})\|_{L^{r}(\mathcal{C})\times
L^{s}(\mathcal{G})}\\}.$

Then there exists a constant $c>0$ independent of $\bar{Q}$,
$(U_{\Omega},U_{\Gamma})\in\mathcal{U}_{\rm{reg}}$, and a constant $M_{0}$
such that

(30) $\|\bar{Q}\|_{C(\overline{\mathcal{C}})}\leq
c\left(\sqrt{M_{0}}+\|(U_{\Omega},U_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}}+\|Q_{0}\|_{C(\overline{\Omega})}\right).$

######  Remark 3.17.

If $Q_{0}\equiv 0$, then $M_{0}=0$. Moreover, since $c$ is independent of
$(U_{\Omega},U_{\Gamma})$, we can vary on a ball in $\mathcal{U}_{\rm{reg}}$
and obtain a uniform bound on the solution operator $S(U_{\Omega},U_{\Gamma})$
in Eq. 31 below as a mapping from $\mathcal{U}_{\rm{reg}}$ into
$C(\overline{\mathcal{C}})$.

######  Proof 3.18.

See Appendix C.

######  Remark 3.19.

In optimal control, proving the state variable is continuous on
$\overline{\mathcal{C}}$ is often useful for the derivation of optimality
conditions for zero-order bound constraints, as it provides an essential
constraint qualification in the convex setting. However, a constraint of the
type “$Q\geq 0$” is not interesting for the current application. Nevertheless,
the continuity of $\bar{Q}$ provides a justification for constraints of the
type “$|\bar{Q}|=0$” on lower-dimensional manifolds embedded in $\Omega$,
which would correspond to the placement of defects. The analysis of this
challenging type of constraint will be part of future research.

######  Remark 3.20.

The energy estimates and related bounds derived in this section (and the
associated appendices) can be useful for future work on the a priori numerical
analysis of the optimal control problem as they would remain true if we
replace the controls $U_{\Omega}$ and $U_{\Gamma}$ by, e.g., finite element
approximations. At several points, we adjust the coefficients in $M_{0}$,
$M_{1}$, $M_{2}$, the generic constant $c$, and $\rho$ in (25) used throughout
the text. Nevertheless, these arguments should be largely unaffected by the
usage of discrete controls and ultimately stable bounds.

##  4 Existence of optimal controls

We denote the control-to-state operator for the forward problem Eq. 17 by

(31) $S:L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})\to\mathcal{W},$

i.e. $S(U_{\Omega},U_{\Gamma})\in\mathcal{W}$ solves Eq. 17 for any controls
$(U_{\Omega},U_{\Gamma})\in L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$.

######  Theorem 4.1.

Suppose that the set of control constraints $\mathcal{U}_{\rm{ad}}$ is a
nonempty, closed, and convex subset of $L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$. Then the optimal control problem
Eq. 13-Eq. 14 admits a solution.

######  Remark 4.2.

The assumption that $\mathcal{U}_{\rm{ad}}\cap[L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})]$ is closed can be guaranteed in a
variety of contexts, e.g. for pointwise a.e. bound constraints. Moreover, if
the boundary control is independent of time or only applied at a _finite_
number of points in time, as in many applications, then this assumption is
fulfilled.

######  Proof 4.3.

For readability, we set
$\widetilde{\mathcal{U}}_{\rm{ad}}:=\mathcal{U}_{\rm{ad}}\cap[L^{2}(\mathcal{C})\times
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})]$. By assumption, this set is
nonempty, closed, and convex and therefore weakly closed in
$L^{2}(\mathcal{C})\times H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$. In
addition, we restrict the control-to-state operator $S(U_{\Omega},U_{\Gamma})$
to $\widetilde{\mathcal{U}}_{\rm{ad}}$.

By hypothesis, $\widetilde{\mathcal{U}}_{\rm{ad}}\neq\emptyset$. Consequently,
there exists a minimizing sequence
$\left\\{(U_{\Omega,n},U_{\Gamma,n})\right\\}\subset\widetilde{\mathcal{U}}_{\rm{ad}}$
for Eq. 13-Eq. 14. Clearly, $\left\\{(U_{\Omega,n},U_{\Gamma,n})\right\\}$ is
uniformly bounded in $L^{2}(\mathcal{C})\times L^{2}(\mathcal{G})$. Moreover,
since $\left\\{(U_{\Omega,n},U_{\Gamma,n})\right\\}$ is a minimizing sequence,
there exists an $n_{0}\in\mathbb{N}$ such that for all $n\geq n_{0}$:

$\left\\{(U_{\Omega,n},U_{\Gamma,n})\right\\}\subset\\\
\left\\{(U_{\Omega},U_{\Gamma})\in\mathcal{U}_{\rm{ad}}\left|J(S(U_{\Omega},U_{\Gamma}),U_{\Omega},U_{\Gamma})\leq
J(S(U_{\Omega,n_{0}},U_{\Gamma,n_{0}}),U_{\Omega,n_{0}},U_{\Gamma,n_{0}})\right.\right\\}.$

By the definition of $J$, there is a constant $c>0$, such that
$\|U_{\Gamma,n}\|_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}\leq c$, for
all $n\geq 1$. It follows that there exists a subsequence
$\left\\{(U_{\Omega,n_{k}},U_{\Gamma,n_{k}})\right\\}$ and (weak) limit point
$(\bar{U_{\Omega}},\bar{U}_{\Gamma})\in\widetilde{\mathcal{U}}_{\rm{ad}}$.
Finally, it follows in light of Theorem 3.13, Remark 3.15, equation Eq. 25,
and the Aubin-Lions-Lemma that there exists a subsequence $\\{Q^{l}\\}$ with
$Q^{l}:=S(U_{\Omega,n_{k_{l}}},U_{\Gamma,n_{k_{l}}})$ that converges to
$\bar{Q}$ such that

  * •

$Q^{l}\to\bar{Q}$ weakly∗ in $L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})$,

  * •

$Q^{l}\to\bar{Q}$ weakly in $L^{2}(0,t_{\mathrm{f}};\mathbf{V})$,

  * •

$Q^{l}\to\bar{Q}$ weakly∗ in
$L^{\infty}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$,

  * •

$Q^{l}\to\bar{Q}$ weakly in $H^{1}(0,t_{\mathrm{f}};\mathbf{H})$,

  * •

$Q^{l}\to\bar{Q}$ strongly in
$C([0,t_{\mathrm{f}}];L^{6-\varepsilon}(\Omega,\mathbf{S}_{0}))$ for
$\varepsilon\in(0,5]$.

Using analogous arguments to those in Appendix B, it follows that
$\bar{Q}=S(\bar{U_{\Omega}},\bar{U}_{\Gamma})$. Finally, the weak lower-
semicontinuity of $J$ along with the properties of $Q^{l}$ and
$(U_{\Omega,n_{k_{l}}},U_{\Gamma,n_{k_{l}}})$ guarantee that
$(\bar{U_{\Omega}},\bar{U}_{\Gamma})$ is an optimal solution of Eq. 13-Eq. 14.

##  5 First-Order Optimality Conditions and the Adjoint Equation

We first derive a differentiability result for $S$.

######  Theorem 5.1.

Under the hypotheses of Theorem 3.16, the control-to-state mapping $S$ from
$\mathcal{U}_{\rm{reg}}$ into $\mathcal{W}$ is Fréchet differentiable.
Moreover, given $(U_{\Omega},U_{\Gamma})$, $(H_{\Omega},H_{\Gamma})$ in
$\mathcal{U}_{\rm{reg}}$, the derivative of $S$ at $(U_{\Omega},U_{\Gamma})$
in direction $(H_{\Omega},H_{\Gamma})$ is given by the unique solution
$\Xi=S^{\prime}_{U_{\Omega},U_{\Gamma}}(H_{\Omega},H_{\Gamma})$ of Eq. 32.

######  Proof 5.2.

Let $(U_{\Omega},U_{\Gamma}),(H_{\Omega},H_{\Gamma})\in
L^{2}(\mathcal{C})\times L^{2}(\mathcal{G})$, denote
$\bar{Q}:=S(U_{\Omega},U_{\Gamma})$,
$Q_{H}:=S(U_{\Omega}+H_{\Omega},U_{\Gamma}+H_{\Gamma})$, and let $\Xi$ be the
solution to

(32)
$\begin{split}\Xi_{t}-\Delta\Xi+\frac{1}{\eta_{\mathrm{dw}}^{2}}\psi^{\prime\prime}(\bar{Q}):\Xi&=\eta_{\Omega}H_{\Omega},\text{
in
}\Omega,\quad\partial_{\bm{\nu}}\Xi+\eta_{\Gamma}\Xi=\eta_{\Gamma}H_{\Gamma},\text{
on }\Gamma,\end{split}$

with $\Xi(\cdot,0)=0$. It follows from Theorem 3.16 that
$\psi^{\prime\prime}(\bar{Q})$ is continuous on $\overline{\mathcal{C}}$.
Next, we set $A_{H}:=Q_{H}-\bar{Q}-\Xi$ and show that it behaves like
$o(\|(H_{\Omega},H_{\Gamma})\|)$ in the appropriate norms. Almost everywhere
on $\mathcal{C}$, we have

$\begin{split}\psi^{\prime}(Q_{H})-\psi^{\prime}(\bar{Q})-\psi^{\prime\prime}(\bar{Q}):\Xi=\psi^{\prime\prime}(\bar{Q}):A_{H}-X_{H},\end{split}$

where
$X_{H}:=-\int_{0}^{1}\psi^{\prime\prime}(\bar{Q}+\tau(Q_{H}-\bar{Q}))-\psi^{\prime\prime}(\bar{Q})d\tau:[Q_{H}-\bar{Q}]$.
With the extended assumptions of Theorem 3.13, $A_{H}$ satisfies

(33) $\begin{split}(A_{H})_{t}-\Delta
A_{H}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\psi^{\prime\prime}(\bar{Q}):A_{H}&=\frac{1}{\eta_{\mathrm{dw}}^{2}}X_{H},\text{
in }\Omega,~{}\partial_{\bm{\nu}}A_{H}+\eta_{\Gamma}A_{H}=0,\text{ on
}\Gamma,\end{split}$

with $A_{H}(\cdot,0)=0$. The system Eq. 33 is a simplified version of the
nonlinear forward problem. Therefore, using a slight modification of the same
arguments, we can prove that $A_{H}\in\mathcal{W}$. In particular, we readily
obtain the following bound (for a generic constant $c>0$ independent of
$H=(H_{\Omega},H_{\Gamma})$): $\|A_{H}\|_{\mathcal{W}}\leq
c\|(X_{H},0)\|_{L^{2}(\mathcal{C})\times L^{2}(\mathcal{G})}.$

As a consequence of the Lipschitz continuity of $\psi^{\prime\prime}$, we have
(a.e. on $\mathcal{C}$)

$|X_{H}|\leq\int_{0}^{1}|\psi^{\prime\prime}(\bar{Q}+\tau(Q_{H}-\bar{Q}))-\psi^{\prime\prime}(\bar{Q})|d\tau|Q_{H}-\bar{Q}|,\leq
c|Q_{H}-\bar{Q}|^{2},$

where $c$ represents the Lipschitz modulus for $\psi^{\prime\prime}$. Note
that $\psi^{\prime\prime\prime}(Q)$ is only nonzero on the set where
$\mathrm{tr}(Q^{2})=|Q|^{2}$ is between two fixed constants $b_{1}$ and
$b_{2}$, i.e., where $Q$ is bounded in space-time. Then for a.e.
$t\in(0,t_{\mathrm{f}})$ we obtain the bound: $\begin{aligned}
\|X_{H}(t)\|^{2}_{\mathbf{H}}\leq
c\|Q_{H}(t)-\bar{Q}(t)\|^{4}_{L^{4}(\Omega)}.\end{aligned}$ It follows from
Eq. 26, the definition of $\mathcal{W}$, and the Sobolev embedding theorem
that $\|Q_{H}(t)-\bar{Q}(t)\|_{L^{p}(\Omega)}\leq
c\|(H_{\Omega},H_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}},$ for a.e.
$t\in(0,t_{\mathrm{f}})$ and any $p\in[1,6]$. Consequently, we have

$\displaystyle\|X_{H}(t)\|^{2}_{\mathbf{H}}$ $\displaystyle\leq
c\|(H_{\Omega},H_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}}^{2}\|Q_{H}(t)-\bar{Q}(t)\|^{2}_{L^{4}(\Omega)}.$

Integrating in time and taking the square root, we obtain

$\displaystyle\|X_{H}\|_{L^{2}(\mathcal{G})}$ $\displaystyle\leq
c\|(H_{\Omega},H_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}}\|Q_{H}-\bar{Q}\|_{L^{2}(0,t_{\mathrm{f}};L^{4}(\Omega))},$

which immediately yields $\|(X_{H},0)\|_{L^{2}(\mathcal{C})\times
L^{2}(\mathcal{G})}=o(\|(H_{\Omega},H_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}})$.

######  Corollary 5.3.

Under the assumptions of Theorem 5.1, the reduced objective function
$\mathcal{J}:\mathcal{U}_{\rm{reg}}\to\overline{\mathbb{R}}$ defined by

(34) $\displaystyle\mathcal{J}(U_{\Omega},U_{\Gamma}):=$
$\displaystyle\frac{\beta_{\mathcal{C}}}{2}\|S(U_{\Omega},U_{\Gamma})-Z_{\mathcal{C}}\|^{2}_{L^{2}(\mathcal{C})}+\frac{\beta_{\mathcal{G}}}{2}\|S(U_{\Omega},U_{\Gamma})|_{\Gamma}-Z_{\mathcal{G}}\|^{2}_{L^{2}(\mathcal{G})}$
$\displaystyle+\frac{\beta_{t_{\mathrm{f}}}}{2}$
$\displaystyle\|S(U_{\Omega},U_{\Gamma})(\cdot,T)-Z_{t_{\mathrm{f}}}\|^{2}_{\mathbf{H}}+\frac{\alpha_{\mathcal{C}}}{2}\|U_{\Omega}\|^{2}_{L^{2}(\mathcal{C})}+\frac{\alpha_{\mathcal{G}}}{2}\|U_{\Gamma}\|^{2}_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})},$

is Fréchet differentiable. Furthermore, given
$\bar{Q}=S(\bar{U}_{\Omega},\bar{U}_{\Gamma})$, a direction
$(H_{\Omega},H_{\Gamma})\in\mathcal{U}_{\rm{reg}}$, and
$\Xi=S^{\prime}_{\bar{U}_{\Omega},\bar{U}_{\Gamma}}(H_{\Omega},H_{\Gamma})$,
the associated solution of Eq. 32, the directional derivative of $\mathcal{J}$
at $(\bar{U}_{\Omega},\bar{U}_{\Gamma})$ in direction
$(H_{\Omega},H_{\Gamma})$ is given by:

(35) $\displaystyle\mathcal{J}_{\bar{U}_{\Omega},\bar{U}_{\Gamma}}^{\prime}$
$\displaystyle(H_{\Omega},H_{\Gamma})=\beta_{\mathcal{C}}\left(\bar{Q}-Z_{\mathcal{C}},\Xi\right)_{L^{2}(\mathcal{C})}+\beta_{\mathcal{G}}\left(\bar{Q}_{\Gamma}-Z_{\mathcal{G}},\Xi\right)_{L^{2}(\mathcal{G})}$
$\displaystyle+\beta_{t_{\mathrm{f}}}$
$\displaystyle\left(\bar{Q}(t_{\mathrm{f}})-Z_{t_{\mathrm{f}}},\Xi(t_{\mathrm{f}})\right)_{\mathbf{H}}+\alpha_{\mathcal{C}}\left(\bar{U}_{\Omega},H_{\Omega}\right)_{L^{2}(\mathcal{C})}+\alpha_{\mathcal{G}}\left(\bar{U}_{\Gamma},H_{\Gamma}\right)_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}.$

Moreover, if $R$ is the unique weak solution of the linear parabolic (adjoint)
equation:

(36a) $\displaystyle-R_{t}-\Delta
R+\frac{1}{\eta_{\mathrm{dw}}^{2}}\psi^{\prime\prime}(\bar{Q}):R$
$\displaystyle=\beta_{\mathcal{C}}(\bar{Q}-Z_{\mathcal{C}}),\qquad~{}~{}\text{
in }\Omega\times(0,t_{\mathrm{f}})$ (36b)
$\displaystyle\partial_{\bm{\nu}}R+\eta_{\Gamma}R$
$\displaystyle=\beta_{\mathcal{G}}(\bar{Q}_{\Gamma}-Z_{\mathcal{G}}),\qquad\text{
on }\Gamma\times(0,t_{\mathrm{f}})$ (36c) $\displaystyle
R(\cdot,t_{\mathrm{f}})$
$\displaystyle=\beta_{t_{\mathrm{f}}}(\bar{Q}(t_{\mathrm{f}})-Z_{t_{\mathrm{f}}}),\quad\text{
in }\Omega,$

then we have

(37)
$\begin{split}\mathcal{J}^{\prime}_{\bar{U}_{\Omega},\bar{U}_{\Gamma}}(H_{\Omega},H_{\Gamma})&=\eta_{\Omega}\left(R,H_{\Omega}\right)_{L^{2}(\mathcal{C})}+\eta_{\Gamma}\left(R,H_{\Gamma}\right)_{L^{2}(\mathcal{G})}\\\
&\qquad\quad+\alpha_{\mathcal{C}}\left(\bar{U}_{\Omega},H_{\Omega}\right)_{L^{2}(\mathcal{C})}+\alpha_{\mathcal{G}}\left(\bar{U}_{\Gamma},H_{\Gamma}\right)_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}.\end{split}$

######  Proof 5.4.

The differentiability of the reduced objective functional is a consequence of
Theorem 5.1, the smoothness of the original tracking-type functional, and the
chain rule. This yields Eq. 35. For the equivalent characterization Eq. 37, we
use Eq. 35 and the adjoint equations Eq. 36 by following the standard
computations for the adjoint calculus, see e.g., [54].

Theorem 5.1 and Corollary 5.3 provide us with first-order optimality
conditions of primal and dual type and a means of efficiently calculating
derivatives of the reduced objective functional, which are needed for
numerical methods.

######  Theorem 5.5.

In addition to Eq. 15, suppose $Q_{0}\in\mathbf{V}\cap C(\overline{\Omega})$
and $\mathcal{U}_{\rm{ad}}\cap\mathcal{U}_{\rm{reg}}\neq\emptyset.$ If the
optimal solution $(\bar{U}_{\Omega},\bar{U}_{\Gamma})$ of Eq. 13-Eq. 14 is in
$\mathcal{U}_{\rm{ad}}\cap\mathcal{U}_{\rm{reg}}$, then the following
variational inequality holds:

(38)
$\eta_{\Omega}\left(R,U_{\Omega}-\bar{U}_{\Omega}\right)_{L^{2}(\mathcal{G})}+\eta_{\Gamma}\left(R,U_{\Gamma}-\bar{U}_{\Gamma}\right)_{L^{2}(\mathcal{G})}\\\
+\alpha_{\mathcal{C}}\left(\bar{U}_{\Omega},U_{\Omega}-\bar{U}_{\Omega}\right)_{L^{2}(\mathcal{C})}+\alpha_{\mathcal{G}}\left(\bar{U}_{\Gamma},U_{\Gamma}-\bar{U}_{\Gamma}\right)_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}\geq
0,$

for all
$(U_{\Omega},U_{\Gamma})\in\mathcal{U}_{\rm{ad}}\cap\mathcal{U}_{\rm{reg}}$,
where $R$ solves Eq. 36 with $\bar{Q}=S(\bar{U}_{\Omega},\bar{U}_{\Gamma})$.

######  Proof 5.6.

This is an immediate consequence of Theorem 5.1 and Corollary 5.3. To see
this, note that
$\mathcal{J}(\bar{U}_{\Omega},\bar{U}_{\Gamma})\leq\mathcal{J}(U_{\Omega},U_{\Gamma})\quad\forall(U_{\Omega},U_{\Gamma})\in\mathcal{U}_{\rm{ad}}$.
By assumption, $\mathcal{U}_{\rm{ad}}\cap\mathcal{U}_{\rm{reg}}$ is a nonempty
convex set. Therefore, the previous relation gives us the difference quotients

$0\leq\lambda^{-1}(\mathcal{J}(\bar{U}_{\Omega}+\lambda
H_{\Omega},\bar{U}_{\Gamma}+\lambda
H_{\Gamma})-\mathcal{J}(\bar{U}_{\Omega},\bar{U}_{\Gamma})),$

where $\lambda\in(0,1)$ and
$(H_{\Omega},H_{\Gamma})=(U_{\Omega},U_{\Gamma})-(\bar{U}_{\Omega},\bar{U}_{\Gamma})$
with
$(U_{\Omega},U_{\Gamma})\in\mathcal{U}_{\rm{ad}}\cap\mathcal{U}_{\rm{reg}}$.
The rest follows from Theorem 5.1 and Corollary 5.3; in particular Eq. 37.

##  6 Finite Element Approximation

We discretize Eq. 17 in the following way. First, we assume that $\Omega$ is
polyhedral so that it can be represented exactly by a conforming triangulation
$\mathcal{T}_{h}=\\{T_{i}\\}$ of shape regular simplices (e.g. tetrahedra),
where $h=\max_{T\in\mathcal{T}_{h}}\mathrm{diam}(T)$. In other words,
$\Omega\equiv\cup_{T\in\mathcal{T}_{h}}T$. Curved domains can also be
considered; the polyhedral assumption is only for simplicity.

Next, we define the space of continuous piecewise polynomial functions on
$\Omega$: $\mathbb{M}_{h}^{k}(\Omega):=\left\\{v\in C^{0}(\Omega)\mid
v|_{T}\in\mathcal{P}_{k}(T),~{}\forall T\in\mathcal{T}_{h}\right\\}$, for
$k\geq 1$, and we reserve $\mathbb{M}_{h}^{0}(\Omega)$ for piecewise constant
functions. Let $\\{E^{i}\\}_{i=1}^{5}$ be a basis of $\mathbf{S}_{0}$. We then
define the following continuous, piecewise linear approximation of
$\mathbf{V}$:

(39) $\begin{split}\mathbf{V}_{h}:=\left\\{P\in
C^{0}(\Omega;\mathbf{S}_{0})\mid
P=\sum_{i=1}^{5}p_{i,h}E^{i},~{}p_{i,h}\in\mathbb{M}_{h}^{1}(\Omega),~{}1\leq
i\leq 5\right\\}\subset\mathbf{V},\end{split}$

and denote by $\mathcal{I}_{h}$ the standard Lagrange interpolant on
$\mathbf{V}_{h}$. Therefore, we approximate $Q\in
L^{2}(0,t_{\mathrm{f}};\mathbf{V})$ by $Q_{h}\in
H^{1}(0,t_{\mathrm{f}};\mathbf{V}_{h})$, i.e. piecewise linear in time; c.f.
[19]. We also introduce the following piecewise constant approximations of
$\mathbf{H}$ and $\mathbf{H}_{\Gamma}$, respectively, for approximating the
controls $U_{\Omega}$, $U_{\Gamma}$:

(40) $\begin{split}\mathbf{H}_{h}&:=\left\\{P\in
L^{2}(\Omega;\mathbf{S}_{0})\mid
P=\sum_{i=1}^{5}p_{i,h}E^{i},~{}p_{i,h}\in\mathbb{M}_{h}^{0}(\Omega),~{}1\leq
i\leq 5\right\\}\subset\mathbf{H},\\\ \mathbf{H}_{\Gamma,h}&:=\left\\{P\in
L^{2}(\Gamma;\mathbf{S}_{0})\mid
P=\sum_{i=1}^{5}p_{i,h}E^{i},~{}p_{i,h}\in\mathbb{M}_{h}^{0}(\Gamma),~{}1\leq
i\leq 5\right\\}\subset\mathbf{H}_{\Gamma},\end{split}$

where $\mathbb{M}_{h}^{0}(\Gamma):=\left\\{v\in L^{2}(\Gamma)\mid
v|_{F}\in\mathcal{P}_{0}(F),~{}\forall F\in\mathcal{F}_{h}\right\\}$, where
$\mathcal{F}_{h}=\\{F\\}$ is the set of faces that make up $\partial\Omega$.
Thus, we approximate $U_{\Omega}\in L^{2}(0,t_{\mathrm{f}};\mathbf{H})$,
$U_{\Gamma}\in H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$ by
$U_{\Omega,h}\in H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{h})$, $U_{\Gamma,h}\in
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma,h})$, respectively. The control
bounds are enforced at the nodal degrees-of-freedom of $\mathbf{H}_{h}$ and
$\mathbf{H}_{\Gamma,h}$, i.e. at the centroid of the mesh elements.

Furthermore, we discretize the time interval $[0,t_{\mathrm{f}}]$ into a union
of $K$ sub-intervals of uniform length $\delta t$, i.e. time-steps. With this,
we write $Q_{h}^{k}(x):=Q_{h}(x,k\delta t)$, and approximate $Q_{t}(x,k\delta
t)$ by the finite difference quotient: $\delta
t^{-1}\left(Q_{h}^{k+1}(x)-Q_{h}^{k}(x)\right)$. In addition, the time-
dependence of the controls is written
$U_{\Omega,h}^{k}(x):=U_{\Omega,h}(x,k\delta t)$,
$U_{\Gamma,h}^{k}(x):=U_{\Gamma,h}(x,k\delta t)$.

The fully discrete version of Eq. 17 is as follows. Given the initial
condition $Q^{0}(\cdot):=\mathcal{I}_{h}Q(\cdot,t=0)$, and controls
$\\{U_{\Omega,h}^{k}\\}_{k=0}^{K}\subset\mathbf{H}_{h}$,
$\\{U_{\Gamma,h}^{k}\\}_{k=0}^{K}\in\mathbf{H}_{\Gamma,h}$, we iteratively
solve the following implicit equation for $k=0,...,K-1$: find
$Q_{h}^{k+1}\in\mathbf{V}_{h}$ such that

(41) $\begin{split}\delta
t^{-1}\left(Q_{h}^{k+1}-Q_{h}^{k},P_{h}\right)_{\mathbf{H}}+&\left(\nabla
Q_{h}^{k+1},\nabla
P_{h}\right)_{\mathbf{H}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\left(\psi^{\prime}(Q_{h}^{k+1}),P_{h}\right)_{\mathbf{H}}\\\
+\eta_{\Gamma}\left(Q_{h}^{k+1},P_{h}\right)_{\mathbf{H}_{\Gamma}}=\eta_{\Gamma}&\left(U_{\Gamma,h}^{k+1},P_{h}\right)_{\mathbf{H}_{\Gamma}}+\eta_{\Omega}\left(U_{\Omega,h}^{k+1},P_{h}\right)_{\mathbf{H}},~{}~{}\forall~{}P_{h}\in\mathbf{V}_{h}.\end{split}$

For $\delta t$ sufficiently small, depending on $\eta_{\mathrm{dw}}$, Eq. 41
is monotone at each time-step and can be effectively solved with Newton’s
method. Similar to Eq. 5, Eq. 41 is a tensor-valued version of a discrete
Allen-Cahn equation with Robin boundary conditions. Convergence of $Q_{h}$ to
the exact solution $Q$ of Eq. 5 follows from the standard theory for semi-
linear parabolic problems; c.f. [19, 52, 62, 61]. The adjoint problem is
solved in an analogous way, using a similar discretization; since the adjoint
PDE Eq. 36 is linear (variable coefficient), Newton’s method is not required.

##  7 Numerical Results

We approximate minimizers of Eq. 13, by discretizing the forward problem in
Eq. 17 with the finite element method described in Section 6. Moreover, the
time integrals present in Eq. 13 are discretized with the trapezoidal rule.
This leads to a discrete form of the adjoint problem in Eq. 36 along with the
corresponding discrete form of the derivative functional Eq. 37. Thus, we use
a projected gradient optimization method, with a back-tracking line-search,
see e.g., [6, 23], to compute (discrete) optimal solutions of Eq. 13. During
the line-search, we compute the projection onto the convex set in Eq. 14 by
straightforward normalization of the current guess for $U_{\Omega}$ and
$U_{\Gamma}$. The entire algorithm was implemented in NGSolve [49].

We present examples when the dimension $d$ is $2$ or $3$. Our experiments
involve tensor quantities that are uniaxial (recall Eq. 3 when $d=3$). For any
$d>1$, a uniaxial $Q$ has the form
$Q_{ij}=s_{*}\left(n_{i}n_{j}-\delta_{ij}/d\right)$, for $1\leq i,j\leq d$,
where $[n_{i}]_{i=1}^{d}\equiv n\in\mathbb{R}^{d}$ is a unit vector (often
called the _director_) and $s_{*}$ depends on the coefficients in $\psi$.

The concept of _defect_ is ubiquitous in liquid crystals and plays a critical
role in our numerical experiments. Assuming that $Q$ has a uniaxial form, a
defect corresponds to a discontinuity in the director $n$. Let $d=2$ and
suppose $n:\mathbb{R}^{2}\to\mathbb{R}^{2}$ is a vector field defined in the
plane, continuous everywhere except at isolated points. The _index_ of $n$,
about a point of discontinuity $p_{0}$, is simply the number of full rotations
of $n(a(s))$ along a closed path $a(s)$ that surrounds $p_{0}$ (see [22, pg.
280]). For vector fields, the index is always an integer. If
$\hat{n}:\mathbb{R}^{2}\to\mathbb{RP}^{1}$, i.e. $\hat{n}(x)\equiv-\hat{n}(x)$
(also known as a line field), then the index may be a half-integer. One can
represent a line field $\hat{n}$ with a vector field $n$ (see [3]), and vice-
versa, in the sense that $\hat{n}\otimes\hat{n}=n\otimes
n\equiv[n_{i}n_{j}]_{i,j=1}^{d}$. Thus, since all $Q$-tensors are uniaxial in
dimension $d=2$, and the algebraic form of a uniaxial $Q$ involves $n\otimes
n$, the _degree_ of the defect of $Q$ at $p_{0}$ is simply the index of the
director $n$ (or equivalently $\hat{n}$) about $p_{0}$.

In dimension $d=3$, the degree of a defect makes sense relative to a plane in
$\mathbb{R}^{3}$. For instance, if the set of defects (points of
discontinuity) forms a $C^{1}$ curve, $\alpha$, in $\mathbb{R}^{3}$, then the
degree of a point on that curve is computed relative to the normal plane of
the curve. In other words, let $n$ be the eigenvector of $Q$ with largest
eigenvalue and define the degree of the defect to be the index of $n$ with
respect to a closed curve (in the normal plane) around a point of
discontinuity in $\alpha$.

However, the LdG model will not create point (line) discontinuities in
dimension $d=2$ ($d=3$) because $Q(t,\cdot)\in H^{1}(\Omega;\mathbf{S}_{0})$.
Therefore, any potential discontinuities get smoothed out causing $Q$ to
vanish there (i.e. the liquid crystal “melts”). Thus, in the LdG model, the
location of defects are usually identified with regions where $Q_{ij}=0$ for
$1\leq i,j\leq d$. For more information on defects, see [13, 12, 50, 56, 29,
39, 45, 36, 10].

###  7.1 Control of a $+1/2$ degree point defect in two dimensions

The domain is the unit square $\Omega=(0,1)^{2}$ and the parameters of the
forward problem are as follows. The coefficients of the double well in Eq. 2
are

(42) $\begin{split}a_{0}=1,\quad a_{2}=16.32653061225,\quad a_{3}=0,\quad
a_{4}=66.63890045814,\end{split}$

and $\psi(Q)$ has a global minimum at
$Q_{*}=s_{*}\left[n_{i}n_{j}-\delta_{ij}/2\right]_{i,j=1}^{2}$, where
$n\in\mathbb{R}^{2}$ is any unit vector, and $s_{*}=0.7$. The other
coefficients are given by $\eta_{\mathrm{dw}}=0.2$, $\eta_{\Omega}=0$,
$\eta_{\Gamma}=100$.

The initial condition was defined as follows. First, let $n=n(x_{1},x_{2})$ be
given by

(43)
$\begin{split}n&=\left(\cos\frac{\theta[0.5,0.5]}{2},\sin\frac{\theta[0.5,0.5]}{2}\right),\quad\theta[a,b](x_{1},x_{2}):=\mbox{atan2}\left(\frac{x_{2}-b}{x_{1}-a}\right),\end{split}$

where atan2 is the four-quadrant inverse tangent function and brackets $[a,b]$
indicate parameters. In other words, $n\otimes n$ corresponds to a $+1/2$
degree defect centered at $(0.5,0.5)$. Next, we set
$r[a,b](x_{1},x_{2})=|(x_{1}-a,x_{2}-b)|$ and

(44)
$\begin{split}Q^{0}:=\frac{r^{2}[0.5,0.5]}{r^{2}[0.5,0.5]+\delta^{2}}s_{*}\left[n_{i}n_{j}-\delta_{ij}/2\right]_{i,j=1}^{2},\end{split}$

where $\delta=\eta_{\mathrm{dw}}/4$; this ensures that $Q^{0}\in
H^{1}(\Omega;\mathbf{S}_{0})\cap C^{0}(\Omega;\mathbf{S}_{0})$. The final time
is $t_{\mathrm{f}}=0.4$ and the time-step is $\delta t=0.004$.

The control parameters in Eq. 13 are

(45)
$\begin{split}\beta_{\mathcal{C}}=1.0,\quad\beta_{\mathcal{G}}=0.0,\quad\beta_{t_{\mathrm{f}}}=1.0,\quad\alpha_{\mathcal{C}}=0.0,\quad\alpha_{\mathcal{G}}=0.01,\end{split}$

and the targets are given by setting
$z=\left(\cos\frac{\theta[0.25,0.35]}{2},\sin\frac{\theta[0.25,0.35]}{2}\right)$
and

(46)
$\begin{split}Z_{\mathcal{C}}=Z_{t_{\mathrm{f}}}=\frac{r^{2}[0.25,0.35]}{r^{2}[0.25,0.35]+\delta^{2}}s_{*}\left[z_{i}z_{j}-\delta_{ij}/2\right]_{i,j=1}^{2},\quad
Z_{\mathcal{G}}=0.\end{split}$

In other words, the control objective is to drive $Q$ toward a state that has
a $+1/2$ degree defect located at the coordinates $(0.25,0.35)$.

In this example, we set $U_{\Omega}\equiv 0$, so we only optimize the boundary
control $U_{\Gamma}$ which we enforce to be time-independent. The initial
guess for optimizing the control is given by setting
$u=\left(\cos\theta[0.5,0.5],\sin\theta[0.5,0.5]\right)$ and

(47)
$U_{\Gamma}=\frac{r^{2}[0.5,0.5]}{r^{2}[0.5,0.5]+\delta^{2}}s_{*}\left[u_{i}u_{j}-\delta_{ij}/2\right]_{i,j=1}^{2}.$

Note: we enforce the convex constraint in Eq. 14 with a projected gradient
method.

Figure 1 shows the performance of our gradient descent method. The
$U_{\Gamma}$ residual is computed as follows. Let $P_{\Gamma}^{k}$ satisfy
$\left(P_{\Gamma}^{k},H_{\Gamma}\right)_{L^{2}(\Gamma)}=-\mathcal{J}_{0,U_{\Gamma}^{k}}^{\prime}(0,H_{\Gamma})$,
for all $H_{\Gamma}\in\mathbf{H}_{\Gamma,h}$ (note: we treat the control as
time-independent here), where $k$ is the optimization iteration. In other
words, $P_{\Gamma}^{k}$ is the $L^{2}(\Gamma)$ projection of the negative
gradient. Next, let $\Pi_{\Gamma}$ be the projection onto the boundary control
part of the convex set in Eq. 14. Then the $U_{\Gamma}$ residual, at the
$k$-th iteration, is defined as
$\|U_{\Gamma}^{k}-\Pi_{\Gamma}\left(U_{\Gamma}^{k}+P_{\Gamma}^{k}\right)\|_{L^{2}(\Gamma)}$.
The computed boundary controls $U_{\Gamma}$ at later iterations do not exhibit
any active set, i.e. the inequality constraint is not active. However, we
found that removing the constraint yielded an optimal $U_{\Gamma}$ that was
_not physical_ , i.e. the eigenvalues of $U_{\Gamma}$ were outside the
physical range (recall the discussion around Eq. 14). Thus, it is necessary to
enforce the inequality constraint during the line-search.

Figure 1: Optimization history (Section 7.1). The $U_{\Gamma}$ residual is
described in the text.

Figure 2 shows the target $Z_{\mathcal{C}}$ and optimized boundary control
$U_{\Gamma}$.

Figure 2: Target $Z_{\mathcal{C}}$ (left) and optimized boundary control
$U_{\Gamma}$ (right) (Section 7.1). We visualize $Z_{\mathcal{C}}$ by plotting
line segments that correspond to the eigenvector of $Z_{\mathcal{C}}$ with
maximum eigenvalue; $U_{\Gamma}$ is visualized similarly. Note how the
boundary control mimics the boundary conditions of the target.

Figure 3 shows the initial and final state of $Q$ that clearly demonstrates
the efficacy of the control.

Figure 3: Initial state for $Q(t=0)$ (left) and final state for
$Q(t=t_{\mathrm{f}})$ (right) (Section 7.1). Line segments correspond to the
eigenvector of $Q$ with maximum eigenvalue; the color scale is based on the
maximum eigenvalue. The position of the point defect (at the final time) is
almost exactly the same as in the target.

###  7.2 Prevent $+1/2$ and $-1/2$ degree point defects from annihilating in
two dimensions

Most of the parameters are the same as in Section 7.1 with the following
modifications. The initial condition is given by first defining:

(48)
$\begin{split}n&:=\left(\cos\left(\frac{\theta[0.4,0.505]}{2}+\frac{\pi}{2}\right),\sin\left(\frac{\theta[0.4,0.505]}{2}+\frac{\pi}{2}\right)\right),\\\
m&:=\left(\cos-\frac{\theta[0.6,0.495]}{2},\sin-\frac{\theta[0.6,0.495]}{2}\right),\end{split}$

i.e. $n\otimes n$ corresponds to a $+1/2$ degree defect centered at
$(0.4,0.505)$ and $m\otimes m$ corresponds to a $-1/2$ degree defect centered
at $(0.6,0.495)$. Then, we set

(49)
$\begin{split}Q_{n}^{0}&:=\frac{r^{2}[0.4,0.505]}{r^{2}[0.4,0.505]+\delta^{2}}s_{*}\left[n_{i}n_{j}-\delta_{ij}/2\right]_{i,j=1}^{2},\\\
Q_{m}^{0}&:=\frac{r^{2}[0.6,0.495]}{r^{2}[0.6,0.495]+\delta^{2}}s_{*}\left[m_{i}m_{j}-\delta_{ij}/2\right]_{i,j=1}^{2},\end{split}$

where $\delta=\eta_{\mathrm{dw}}/4$; this ensures that $Q_{n}^{0},Q_{m}^{0}\in
H^{1}(\Omega;\mathbf{S}_{0})\cap C^{0}(\Omega;\mathbf{S}_{0})$. Then, the
initial condition is given by the following interpolation:

(50) $\begin{split}Q^{0}&:=(1-x)Q_{n}^{0}+xQ_{m}^{0}.\end{split}$

The control parameters in Eq. 13 are the same as in Eq. 45, and the targets
$Z_{\mathcal{C}}$, $Z_{t_{\mathrm{f}}}$ have the same form as Eq. 50, except
the $+1/2$ defect is placed at $(0.2,0.6)$ and the $-1/2$ defect is placed at
$(0.8,0.4)$. Note that $Z_{\mathcal{G}}=0$ plays no role. In other words, the
control objective is to drive $Q$ toward a stable configuration of a $+1/2$
and $-1/2$ defect. In this example, we set $U_{\Omega}\equiv 0$, so we only
optimize the boundary control $U_{\Gamma}$ which we enforce to be time-
independent. The initial guess for optimizing the control is the constant
tensor
$U_{\Gamma}=s_{*}\left[\hat{n}_{i}\hat{n}_{j}-\delta_{ij}/2\right]_{i,j=1}^{2}$,
where $\hat{n}=(1,0)$. In this case, the $Q$ state evolves toward a constant
state identical to the initial boundary control, i.e. the two initial defects
annihilate.

In this example, we modify the inequality constraint in Eq. 14 to be $|P|\leq
0.6$ on $\Gamma$. Figure 4 shows the performance of our gradient descent
method. The $U_{\Gamma}$ residual is computed as in Section 7.1. The computed
boundary control $U_{\Gamma}$ does exhibit an active set. Indeed, it was
necessary to lower the bound to $0.6$ in order to ensure that the computed
control satisfied the eigenvalue bounds described in Section 2.1, which in two
dimensions is $-1/2\leq\lambda_{i}(U_{\Gamma})\leq 1/2$, for $i=1,2$. This
further emphasizes that the inequality constraint is needed to prevent
computing minimizers of the objective functional that are not physical (see
the discussion in Section 7.1).

Figure 4: Optimization history (Section 7.2). The $U_{\Gamma}$ residual is
described in the text.

Figure 5 shows the target $Z_{\mathcal{C}}$ and optimized boundary control
$U_{\Gamma}$; note that the maximum value of $\lambda_{1}(U_{\Gamma})$ is
$0.42$.

Figure 5: Target $Z_{\mathcal{C}}$ (left) and optimized boundary control
$U_{\Gamma}$ (right) (Section 7.2). We visualize $Z_{\mathcal{C}}$ by plotting
line segments that correspond to the eigenvector of $Z_{\mathcal{C}}$ with
maximum eigenvalue; $U_{\Gamma}$ is visualized similarly. Note how the
boundary control mimics the boundary conditions of the target.

Figure 6 shows the initial and final state of $Q$ that clearly demonstrates
the efficacy of the control.

Figure 6: Initial state for $Q(t=0)$ (left) and final state for
$Q(t=t_{\mathrm{f}})$ (right) (Section 7.2). Line segments correspond to the
eigenvector of $Q$ with maximum eigenvalue; the color scale is based on the
maximum eigenvalue. The position of the point defect (at the final time) is
almost exactly the same as in the target.

###  7.3 Control of a $+1/2$ line defect in three dimensions

The domain is the unit cube $\Omega=(0,1)^{3}$ and the parameters of the
forward problem are as follows. The coefficients of the double well in Eq. 2
are

(51) $\begin{split}a_{0}=1,\quad a_{2}=7.5021037403,\quad
a_{3}=60.975813166,\quad a_{4}=66.519068908,\end{split}$

and $\psi(Q)$ has a global minimum at
$Q_{*}=s_{*}\left[n_{i}n_{j}-\delta_{ij}/3\right]_{i,j=1}^{3}$, where
$n\in\mathbb{R}^{3}$ is any unit vector, and $s_{*}=0.700005531$. The other
coefficients are given by $\eta_{\mathrm{dw}}=0.2$, $\eta_{\Omega}=0$,
$\eta_{\Gamma}=100$.

The initial condition was defined as follows. First, let $n=n(x_{1},x_{2})$ be
given by

(52)
$\begin{split}n&:=\left(\cos\frac{\theta[0.5,0.5]}{2},\sin\frac{\theta[0.5,0.5]}{2},0\right),\quad\theta[a,b](x_{1},x_{2}):=\mbox{atan2}\left(\frac{x_{2}-b}{x_{1}-a}\right),\end{split}$

similar to Eq. 43. In other words, $n\otimes n$ corresponds to a $+1/2$ degree
defect, in any plane parallel to the $x_{3}=0$ plane, centered at $(0.5,0.5)$.
Then, we have

(53)
$\begin{split}Q^{0}:=\frac{r^{2}[0.5,0.5]}{r^{2}[0.5,0.5]+\delta^{2}}s_{*}\left[n_{i}n_{j}-\delta_{ij}/3\right]_{i,j=1}^{3},\end{split}$

where $\delta=\eta_{\mathrm{dw}}/4$; this ensures that $Q^{0}\in
H^{1}(\Omega;\mathbf{S}_{0})\cap C^{0}(\Omega;\mathbf{S}_{0})$. The final time
is $t_{\mathrm{f}}=0.3$ and the time-step is $\delta t=0.006$.

The control parameters in Eq. 13 are the same as in Eq. 45. The targets are
defined through a parameterized curve in $\mathbb{R}^{3}$, denoted
$\left(\tilde{x}_{1}(\xi),\tilde{x}_{2}(\xi),\tilde{x}_{3}(\xi)\right)$, given
by

(54) $\begin{split}f(\xi)&:=l_{2}\xi^{2}+l_{3}\xi^{3},\quad
l_{2}=3c_{0},~{}l_{3}=-2c_{0},~{}c_{0}=0.6,\\\
\left(\tilde{x}_{1}(\xi),\tilde{x}_{2}(\xi),\tilde{x}_{3}(\xi)\right)&:=\left(f(\xi)+0.2,f(\xi)+0.2,\xi\right),\text{
for }0\leq\xi\leq 1.\end{split}$

Next, we define
$\tilde{r}^{2}(x_{1},x_{2},x_{3}):=\left(x_{1}-\tilde{x}_{1}(x_{3})\right)^{2}+\left(x_{2}-\tilde{x}_{2}(x_{3})\right)^{2}$,

(55)
$\begin{split}\tilde{\theta}(x_{1},x_{2},x_{3})&:=\mbox{atan2}\left(\frac{x_{2}-\tilde{x}_{2}(x_{3})}{x_{1}-\tilde{x}_{1}(x_{3})}\right),\quad
z=\left(\cos\frac{\tilde{\theta}}{2},\sin\frac{\tilde{\theta}}{2},0\right),\end{split}$

and the targets are given by

(56)
$\begin{split}Z_{\mathcal{C}}=Z_{t_{\mathrm{f}}}=\frac{\tilde{r}^{2}}{\tilde{r}^{2}+\delta^{2}}s_{*}\left[z_{i}z_{j}-\delta_{ij}/3\right]_{i,j=1}^{3},\quad
Z_{\mathcal{G}}=0.\end{split}$

In other words, the control objective is to drive $Q$ toward a state that has
a $+1/2$ degree defect, with respect to the $x_{3}=l_{0}$ plane, located at
$(\tilde{x}_{1}(l_{0}),\tilde{x}_{2}(l_{0}),l_{0})$.

In this example, we set $U_{\Omega}\equiv 0$, so we only optimize the boundary
control $U_{\Gamma}$ which we enforce to be time-independent. The initial
guess for optimizing the control is given by setting $U_{\Gamma}=Q^{0}$.

Figure 7 shows the performance of our gradient descent method. The
$U_{\Gamma}$ residual is computed as in Section 7.1. The computed boundary
controls $U_{\Gamma}$ at later iterations do not exhibit any active set, i.e.
the inequality constraint is not active.

Figure 7: Optimization history (Section 7.3). The $U_{\Gamma}$ residual is
described in the text.

Figure 8 shows the target $Z_{\mathcal{C}}$ and optimized boundary control
$U_{\Gamma}$. We note, however, that the most negative eigenvalue,
$\lambda_{3}(Q)$ (not plotted), is approximately $-0.33$ at the core of the
defect in $U_{\Gamma}$ on the $x_{3}=1$ side of the cube (see middle plot of
Fig. 8). Again, it is necessary to enforce the inequality constraint during
the line-search in order to prevent computing minimizers of the objective
functional that are not physical (see the discussion in Section 7.1).

Figure 8: Target $Z_{\mathcal{C}}$ (left) and optimized boundary control
$U_{\Gamma}$ (middle, right) (Section 7.3). We visualize $Z_{\mathcal{C}}$ by
plotting line segments that correspond to the eigenvector of $Z_{\mathcal{C}}$
with maximum eigenvalue; $U_{\Gamma}$ is visualized similarly. Middle (Right)
view shows the front (back) three faces.

Figure 9 shows the initial and final state of $Q$ that clearly demonstrates
the efficacy of the control.

Figure 9: Initial state for $Q(t=0)$ (left) and final state for
$Q(t=t_{\mathrm{f}})$ (middle) (Section 7.3). Line segments correspond to the
eigenvector of $Q$ with maximum eigenvalue; the color scale is based on the
maximum eigenvalue. The position of the line defect (at the final time) is
very close to that in the target. Right plot: the color scale corresponds to
the biaxiality measure (see Remark 7.1). Away from the defect (not shown), $Q$
is essentially uniaxial with $\beta_{\mathrm{bi}}(Q)\approx 0$.

######  Remark 7.1.

In dimension $d=2$, all $Q$-tensors have a uniaxial form. For $d=3$, $Q$ is
uniaxial if and only if $Q$ has two repeated eigenvalues [56]. Moreover, even
if the initial condition $Q_{0}$ is uniaxial, the solution $Q(t,\cdot)$ of Eq.
5 will not be uniaxial in general, i.e. it will become biaxial with three
distinct eigenvalues. Typically, the solution is only biaxial near a defect;
away from a defect, it is (essentially) uniaxial because of the global minimum
properties of the bulk potential in Eq. 2 (see [44]).

The right plot of Fig. 9 shows the normalized biaxiality measure [42],
$\beta_{\mathrm{bi}}(Q)$, of the solution $Q$ on the iso-surface surrounding
the line defect. Note that $0\leq\beta_{\mathrm{bi}}(Q)\leq 1$, with
$\beta_{\mathrm{bi}}(Q)=0$ corresponding to a uniaxial state and
$\beta_{\mathrm{bi}}(Q)=1$ representing “maximal” biaxiality.

##  8 Conclusions

The main contribution of this paper is to show that optimal control of LC
devices, in the framework of the LdG model, is possible. Indeed, our numerical
study demonstrates this effectively by directly controlling the placement of
defects, which is of considerable interest in the LC scientific community. We
only consider boundary controls in our numerical study since that is most
relevant in applications. Further extensions of our framework, as related to
actual LC systems, would involve controls that are either finite-dimensional
(in space) or with a special restriction on the admissible controls, e.g.
homeotropic versus planar anchoring for boundary controls.

From an analytical standpoint, by restricting our study to the one-parameter
model (i.e. the only non-zero elastic constant is $\ell_{1}=1$), we were able
to exploit a large number of derivation techniques for the optimal control of
scalar Allen-Cahn equations. The rigorous proofs for the bounds and energy
estimates in the tensor-valued setting have therefore been relegated to
appendices. Nevertheless, there remain a number of analytical challenges if we
wish to go beyond the one-parameter model. For example, our current proof of
continuity in space-time may only work in the current setting and new
techniques or regularity results also appear necessary. This is because for
more general elastic constants, the Laplacian in Eq. 5 is replaced by a more
general elliptic operator that fully couples all components of the $Q$-tensor.

Finally, our numerical study made use of a basic optimization algorithm. A
more advanced scheme, e.g., one based on second-order information would
require an additional sensitivity result to derive an analytical formula for
second-order directional derivatives (Hessian-vector products) for use in
Newton-type methods. At least for the bulk energy term considered here, such a
result should be obtainable by modifying the proof of differentiability in
Section 5.

##  Appendix A Proof of Proposition 3.9

We start by using the test function $P=Q^{n}_{t}(s)$ for all
$s\in[0,t_{\mathrm{f}}]$ in (19). This leads to

(57) $\|Q^{n}_{t}(s)\|^{2}_{\mathbf{H}}+\frac{d}{dt}\|\nabla
Q^{n}(s)\|^{2}_{\mathbf{H}}+\frac{\eta_{\Gamma}}{2}\frac{d}{dt}\|Q^{n}(s)\|^{2}_{\mathbf{H}_{\Gamma}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\frac{d}{dt}\int_{\Omega}\psi(Q^{n}(s))\\\
=\eta_{\Omega}\left(U_{\Omega}(s),Q^{n}_{t}(s)\right)_{\mathbf{H}}+\eta_{\Gamma}\left(U_{\Gamma}(s),Q^{n}_{t}(s)\right)_{\mathbf{H}_{\Gamma}}.$

We continue (57) by using $U_{\Gamma}\in
H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$, integrating from $0$ to
$t\in(0,t_{\mathrm{f}}]$, and rearranging terms to obtain new constants
$\delta^{\prime}_{1},\dots,\delta^{\prime}_{4}>0$:

(58) $\delta^{\prime}_{4}\|Q^{n}_{t}\|^{2}_{L^{2}(\mathcal{C}_{t})}+\|\nabla
Q^{n}(t)\|^{2}_{\mathbf{H}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\psi(Q^{n}(t))+\delta^{\prime}_{1}\|Q^{n}(t)\|^{2}_{\mathbf{H}_{\Gamma}}\leq\\\
\|\nabla
Q^{n}_{0}\|^{2}_{\mathbf{H}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\psi(Q^{n}_{0})+\delta^{\prime}_{2}\|Q^{n}_{0}\|^{2}_{\mathbf{H}_{\Gamma}}+\delta^{\prime}_{3}\|Q^{n}\|^{2}_{L^{2}(\mathcal{C}_{t})}+M_{2}.$

We can bound the penultimate term in (58) by applying (20) and

$\|Q^{n}\|_{L^{2}(\mathcal{C}_{t})}^{2}\leq\|Q^{n}\|_{L^{2}(\mathcal{C})}^{2}\leq
c_{\rm emb}\|Q^{n}\|_{L^{2}(0,t_{\mathrm{f}};\mathbf{V})}^{2}\leq c_{\rm
emb}M_{1},$

where $c_{\rm emb}$ is an embedding constant. $M_{1}$ absorbs
$\delta^{\prime}_{3}$ and $c_{\rm emb}$ below.

Based on the order of the nonlinearity $\psi$, the continuity of the trace
operator, and the convex splitting $\psi=\psi_{c}-\psi_{e}$, there is a
constant $M_{0}\geq 0$ such that

$\|\nabla
Q^{n}_{0}\|^{2}_{\mathbf{H}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\psi_{c}(Q^{n}_{0})-\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{\Omega}\psi_{e}(Q^{n}_{0})+\delta^{\prime}_{2}\|Q^{n}_{0}\|^{2}_{\mathbf{H}_{\Gamma}}\leq
M_{0}.$

Since $\psi_{c}$ is continuous on $\mathbf{V}$ due to the Sobolev embedding
theorem and $\psi_{e}(Q)$ is quadratic in $Q$, we can pass to the limit in $n$
and thus obtain (21).

Next, since $\psi$ is bounded from below, we can adjust all the constants and
coefficients if necessary to obtain the bound

(59) $\delta^{\prime}_{4}\|Q^{n}_{t}\|^{2}_{L^{2}(\mathcal{C}_{t})}+\|\nabla
Q^{n}(t)\|^{2}_{\mathbf{H}}+\delta^{\prime}_{1}\|Q^{n}(t)\|^{2}_{\mathbf{H}_{\Gamma}}\leq
M_{0}+M_{1}+M_{2}.$

This yields (22). Now, by letting $\varepsilon>0$ be a small positive
constant, we can bound (59) from below, which yields

(60)
$\delta^{\prime}_{4}\|Q^{n}_{t}\|^{2}_{L^{2}(\mathcal{C}_{t})}+\min\\{1-\varepsilon\delta^{\prime}_{1},\varepsilon\frac{\delta^{\prime}_{1}}{\kappa_{0}}\\}\left(\|\nabla
Q^{n}(t)\|^{2}_{\mathbf{H}}+\|Q^{n}(t)\|^{2}_{\mathbf{H}}\right)\leq
M_{0}+M_{1}+M_{2}.$

Here, $\kappa_{0}$ comes from using a Poincaré type inequality. We can now
adjust the coefficients and constants to deduce the bound:

(61) $\|Q^{n}_{t}\|^{2}_{L^{2}(\mathcal{C}_{t})}+\|\nabla
Q^{n}(t)\|^{2}_{\mathbf{H}}+\|Q^{n}(t)\|^{2}_{\mathbf{H}}\leq
c\left(M_{0}+M_{1}+M_{2}\right).$

This yields (23). It follows from (59), (61), and (20) that
$\left\\{Q^{n}\right\\}$ is uniformly bounded in $\mathcal{W}$ (16).

##  Appendix B Proof of Theorem 3.13

The Aubin-Lions-Simon Lemma, see e.g., Theorem II.5.16, pages 102-103 in [11]
provides several helpful statements. We provide brief justifications
afterwards, as these are well-known embeddings.

  1. 1.

There exists a subsequence $\\{Q^{k}\\}$ with $Q^{k}:=Q^{n_{k}}$ that
converges strongly to the function $\bar{Q}$ in
$C([0,t_{\mathrm{f}}];L^{6-\varepsilon}(\Omega,\mathbf{S}_{0}))$ for
$\varepsilon\in(0,5]$.

  2. 2.

There exists a subsequence $\\{Q^{l}\\}$ with $Q^{l}:=Q^{k_{l}}$ that
converges weakly in

$\mathcal{W}_{1}:=\left\\{Q\in
L^{2}(0,t_{\mathrm{f}};\mathbf{V})\left|\;Q_{t}\in
L^{2}(0,t_{\mathrm{f}};\mathbf{H})\right.\right\\}~{}\text{ to }~{}\bar{Q}.$

  3. 3.

There exists a subsequence $\\{Q^{m}\\}$ with $Q^{m}:=Q^{l_{m}}$ that
converges weakly to $\bar{Q}$ in
$L^{2}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$.

The first subsequence exists by the Aubin-Lions Lemma, which implies that

$\mathcal{W}_{2}:=\left\\{Q\in
L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})\left|\;Q_{t}\in
L^{2}(0,t_{\mathrm{f}};\mathbf{H})\right.\right\\},$

is compactly embedded into the space
$C([0,t_{\mathrm{f}}];L^{6-\varepsilon}(\Omega,\mathbf{S}_{0}))$. Here, we
make use of the Sobolev embedding theorem to embed $\mathbf{V}$ into
$L^{6-\varepsilon}(\Omega,\mathbf{S}_{0})$. The second subsequence exists due
to the reflexivity of $\mathcal{W}_{1}$; likewise for the final subsequence.
Finally, we can also argue that $\bar{Q}\in\mathcal{W}$, by appealing to the
bounds in Appendix A, which are stable under passage to the limit in $n$. To
be more specific, we can find an independent constant $\rho>0$ such that

(62)
$\displaystyle\mathrm{ess\,sup}_{t\in[0,t_{\mathrm{f}}]}(\|\nabla\bar{Q}(t)\|^{2}_{\mathbf{H}}+\|\bar{Q}(t)\|^{2}_{\mathbf{H}})^{1/2}$
$\displaystyle\leq\rho,$
$\displaystyle\mathrm{ess\,sup}_{t\in[0,t_{\mathrm{f}}]}\|\bar{Q}(t)\|_{\mathbf{H}_{\Gamma}}$
$\displaystyle\leq\rho,$
$\displaystyle(\|\bar{Q}\|^{2}_{L^{2}(\mathcal{C})}+\|\bar{Q}_{t}\|^{2}_{L^{2}(\mathcal{C})})^{1/2}$
$\displaystyle\leq\rho$
$\displaystyle\|\bar{Q}\|_{L^{2}(0,t_{\mathrm{f}};\mathbf{V})}$
$\displaystyle\leq\rho.$

We arrive at the bound in the space $\mathcal{W}$, i.e. (25).

It remains to show that $\bar{Q}$ is a weak solution of (5). Uniqueness is a
consequence of Theorem 3.2. For arbitrarily fixed data
$(U_{\Omega},U_{\Gamma},Q_{0})$ that satisfies (15) and a test function
$P\in\mathcal{W}_{1}$ such that $P(0)=0$ a.e., we recall (19) and integrate in
$t$:

(63)
$\begin{split}\int_{0}^{t_{\mathrm{f}}}\left(Q^{n}_{t},P\right)_{\mathbf{H}}+&\int_{0}^{t_{\mathrm{f}}}\left(\nabla
Q^{n},\nabla
P\right)_{\mathbf{H}}+\frac{1}{\eta_{\mathrm{dw}}^{2}}\int_{0}^{t_{\mathrm{f}}}\int_{\Omega}\psi^{\prime}(Q^{n}):P\\\
+\eta_{\Gamma}&\int_{0}^{t_{\mathrm{f}}}\left(Q^{n},P\right)_{\mathbf{H}_{\Gamma}}=\eta_{\Omega}\int_{0}^{t_{\mathrm{f}}}\left(U_{\Omega},P\right)_{\mathbf{H}}+\eta_{\Gamma}\int_{0}^{t_{\mathrm{f}}}\left(U_{\Gamma},P\right)_{\mathbf{H}_{\Gamma}}.\end{split}$

The convergence of the linear terms in Eq. 63 follows by straightforward
arguments, e.g., weak convergence and use of compact embeddings. For the
nonlinear term, it suffices to note that $\psi^{\prime}$ is globally
Lipschitz, which provides, e.g., strong convergence in
$C([0,t_{f}];L^{7/4}(\Omega))$ of $\psi^{\prime}(Q^{k})$ to
$\psi^{\prime}(\bar{Q})$. Given $P\in
L^{2}([0,t_{\mathrm{f}}];\mathbf{V})\subset
L^{2}([0,t_{\mathrm{f}}];L^{7/3}(\Omega))$, we can then pass to the limit.

##  Appendix C Proof of Theorem 3.16

After using the bootstrapping and decomposition technique, we derive from [54,
Thm. 5.5] the existence of some constant $c>0$, independent of
$(U_{\Omega},U_{\Gamma})$, for which we have

(64) $\|\bar{Q}\|_{C(\overline{\mathcal{C}})}\leq
c\left(\|\eta_{\Omega}U_{\Omega}-\frac{1}{\eta_{\mathrm{dw}}^{2}}\psi(\bar{Q})\|_{L^{r}(\mathcal{C})}+\|\eta_{\Gamma}U_{\Gamma}\|_{L^{s}(\mathcal{G})}+\|Q_{0}\|_{C(\overline{\Omega})}\right).$

We remove the dependence on $\bar{Q}$ from the right-hand side, by noting that
(8) implies

(65) $\|\psi^{\prime}(\bar{Q})\|_{L^{r}(\mathcal{C})}\leq
c_{1}\|\bar{Q}\|_{L^{r}(\mathcal{C})}\leq c_{1}c_{\rm
emb}\|\bar{Q}\|_{L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})}\leq\rho,$

where the second inequality follows from the continuous embedding of
$L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})$ into $L^{r}(\mathcal{C})$ (provided
$r\in(5/2,6]$), $c_{\rm emb}$ is the associated embedding constant, and $\rho$
is from (62). Next, we derive an explicit bound for $\rho$. Starting from
(61), we note that before passing to the limit in $n$ we have

$\sqrt{\|\nabla
Q^{n}(t)\|^{2}_{\mathbf{H}}+\|Q^{n}(t)\|^{2}_{\mathbf{H}}}\leq\sqrt{M_{0}+M_{1}+M_{2}}\leq\sqrt{M_{0}}+\\\
M_{1}(\|Q^{n}_{0}\|_{\mathbf{H}},\|U_{\Omega}\|_{L^{2}(\mathcal{C})},\|U_{\Gamma}\|_{L^{2}(\mathcal{G})})+M_{2}(\|U_{\Gamma}(t)\|_{\mathbf{H}_{\Gamma}},\|U_{\Gamma}(0)\|_{\mathbf{H}_{\Gamma}},\|(U_{\Gamma})_{t}\|_{L^{2}(\mathcal{G}_{t})},\|u\|_{L^{2}(\mathcal{C}_{t})}).$

Here, we use the subadditivity of $\sqrt{\cdot}$ along with the fact that
$M_{1}$ and $M_{2}$ are simple multilinear maps of their arguments with
positive coefficients. We may then pass to the limit in $n$ along an
appropriate subsequence and obtain the same inequality independent of $n$. The
$M_{0}$-term is independent of $(U_{\Omega},U_{\Gamma})$. The $M_{1}$-term can
be bounded in the first argument by $\|Q_{0}\|_{C(\overline{\Omega})}$ and the
$(U_{\Omega},U_{\Gamma})$-terms by stronger norms. For the $M_{2}$-term we
have several possibilities. Since
$(U_{\Omega},U_{\Gamma})\in\mathcal{U}_{\rm{reg}}$, with norm given by Eq. 29,
and $H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})$ is continuously embedded
into $C([0,t_{\mathrm{f}}];H_{\Gamma})$, we can bound the first two arguments
in $M_{2}$ first by $\|U_{\Gamma}\|_{C([0,t_{\mathrm{f}}];H_{\Gamma})}$ and
then further from above by
$\|U_{\Gamma}\|_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}$. The latter two
arguments can be bounded from above by the norms
$\|(U_{\Gamma})_{t}\|_{L^{2}(\mathcal{G})}$ and
$\|U_{\Omega}\|_{L^{2}(\mathcal{C})}$, respectively. Clearly, the third
argument can be bounded from above by
$\|U_{\Omega}\|_{H^{1}(0,t_{\mathrm{f}};\mathbf{H}_{\Gamma})}$. Since $r>2$,
the fourth argument can be bounded from above by
$\|U_{\Omega}\|_{L^{r}(\mathcal{C})}$. By combining all of these observations,
we deduce the existence of a constant $c>0$, independent of
$\bar{Q},(U_{\Omega},U_{\Gamma})$ and $Q_{0}$ such that for all
$t\in[0,t_{\mathrm{f}}]$ we have

$\sqrt{\|\nabla\bar{Q}(t)\|^{2}_{\mathbf{H}}+\|\bar{Q}(t)\|^{2}_{\mathbf{H}}}\leq\sqrt{M_{0}}+c(\|(U_{\Omega},U_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}}+\|Q_{0}\|_{C(\overline{\Omega})}),$

which implies that
$\|\bar{Q}\|_{L^{\infty}(0,t_{\mathrm{f}};\mathbf{V})}\leq\sqrt{M_{0}}+c(\|(U_{\Omega},U_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}}+\|Q_{0}\|_{C(\overline{\Omega})}).$
Combining this bound with (64) and (65), there exists a constant $c>0$,
independent of $\bar{Q},(U_{\Omega},U_{\Gamma})$ and $Q_{0}$, such that
$\|\bar{Q}\|_{C(\overline{\mathcal{C}})}\leq
c\left(\sqrt{M_{0}}+\|(U_{\Omega},U_{\Gamma})\|_{\mathcal{U}_{\rm{reg}}}+\|Q_{0}\|_{C(\overline{\Omega})}\right)$.
The assertion then follows.

## References

  * [1] Ágnes Buka and N. Éber, eds., Flexoelectricity in Liquid Crystals: Theory, Experiments and Applications, World Scientific, 2012.
  * [2] I. Bajc, F. Hecht, and S. Žumer, A mesh adaptivity scheme on the Landau–de Gennes functional minimization case in 3D, and its driving efficiency, Journal of Computational Physics, 321 (2016), pp. 981 – 996, https://doi.org/https://doi.org/10.1016/j.jcp.2016.02.072, http://www.sciencedirect.com/science/article/pii/S0021999116001443.
  * [3] J. M. Ball and A. Zarnescu, Orientable and non-orientable director fields for liquid crystals, Proceedings in Applied Mathematics and Mechanics (PAMM), 7 (2007), pp. 1050701–1050704, https://doi.org/10.1002/pamm.200700489.
  * [4] G. Barbero and G. Durand, On the validity of the rapini-papoular surface anchoring energy form in nematic liquid crystals, J. Phys. France, 47 (1986), pp. 2129–2134, https://doi.org/10.1051/jphys:0198600470120212900, https://doi.org/10.1051/jphys:0198600470120212900.
  * [5] S. Bartels and A. Raisch, Simulation of Q-tensor fields with constant orientational order parameter in the theory of uniaxial nematic liquid crystals, in Singular Phenomena and Scaling in Mathematical Models, M. Griebel, ed., Springer International Publishing, 2014, pp. 383–412, https://doi.org/10.1007/978-3-319-00786-1_17.
  * [6] D. P. Bertsekas, On the goldstein-levitin-polyak gradient projection method, IEEE Transactions on automatic control, 21 (1976), pp. 174–184.
  * [7] J. Biggins, M. Warner, and K. Bhattacharya, Elasticity of polydomain liquid crystal elastomers, Journal of the Mechanics and Physics of Solids, 60 (2012), pp. 573 – 590, https://doi.org/10.1016/j.jmps.2012.01.008, http://www.sciencedirect.com/science/article/pii/S0022509612000166.
  * [8] L. Blinov, Electro-optical and magneto-optical properties of liquid crystals, Wiley, 1983.
  * [9] J.-P. Borthagaray, R. H. Nochetto, and S. W. Walker, A structure-preserving FEM for the uniaxially constrained $\mathbf{Q}$-tensor model of nematic liquid crystals, Numerische Mathematik, 145 (2020), pp. 837 – 881, https://doi.org/10.1007/s00211-020-01133-z, https://doi.org/10.1007/s00211-020-01133-z.
  * [10] J. P. Borthagaray and S. W. Walker, Chapter 5 - the Q-tensor model with uniaxial constraint, in Geometric Partial Differential Equations - Part II, A. Bonito and R. H. Nochetto, eds., vol. 22 of Handbook of Numerical Analysis, Elsevier, 2021, pp. 313 – 382, https://doi.org/https://doi.org/10.1016/bs.hna.2020.09.001, http://www.sciencedirect.com/science/article/pii/S1570865920300132.
  * [11] F. Boyer and P. Fabrie, Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models, Springer New York, 2013, https://doi.org/10.1007/978-1-4614-5975-0.
  * [12] H. Brezis, J.-M. Coron, and E. H. Lieb, Harmonic maps with defects, Communications in Mathematical Physics, 107 (1986), pp. 649–705, https://doi.org/10.1007/BF01205490.
  * [13] W. F. Brinkman and P. E. Cladis, Defects in liquid crystals, Physics Today, 35 (1982), pp. 48–56.
  * [14] F. Brochard, L. Léger, and R. B. Meyer, Freedericksz transition of a homeotropic nematic liquid crystal in rotating magnetic fields, J. Phys. Colloques, 36 (1975), pp. C1–209–C1–213, https://doi.org/10.1051/jphyscol:1975139.
  * [15] M. Camacho-Lopez, H. Finkelmann, P. Palffy-Muhoray, and M. Shelley, Fast liquid-crystal elastomer swims into the dark, Nature Materials, 3 (2004), pp. 307–310.
  * [16] H. Coles and S. Morris, Liquid-crystal lasers, Nature Photonics, 4 (2010), pp. 676–685.
  * [17] P. Colli and J. Sprekels, Optimal control of an Allen–Cahn equation with singular potentials and dynamic boundary condition, SIAM Journal on Control and Optimization, 53 (2015), pp. 213–234, https://doi.org/10.1137/120902422, https://doi.org/10.1137/120902422, https://arxiv.org/abs/https://doi.org/10.1137/120902422.
  * [18] P. Dasgupta, M. K. Das, and B. Das, Fast switching negative dielectric anisotropic multicomponent mixtures for vertically aligned liquid crystal displays, Materials Research Express, 2 (2015), p. 045015, http://stacks.iop.org/2053-1591/2/i=4/a=045015.
  * [19] T. Davis and E. Gartland, Finite element analysis of the Landau-de Gennes minimization problem for liquid crystals, SIAM Journal on Numerical Analysis, 35 (1998), pp. 336–362, https://doi.org/10.1137/S0036142996297448, https://doi.org/10.1137/S0036142996297448, https://arxiv.org/abs/https://doi.org/10.1137/S0036142996297448.
  * [20] P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, vol. 83 of International Series of Monographs on Physics, Oxford Science Publication, Oxford, UK, 2nd ed., 1995.
  * [21] W. H. de Jeu, ed., Liquid Crystal Elastomers: Materials and Applications, Advances in Polymer Science, Springer, 2012.
  * [22] M. P. do Carmo, Differential Geometry of Curves and Surfaces, Prentice Hall, Upper Saddle River, New Jersey, 1976.
  * [23] J. C. Dunn, Global and asymptotic convergence rate estimates for a class of projected gradient processes, SIAM Journal on Control and Optimization, 19 (1981), pp. 368–400, https://doi.org/10.1137/0319022, https://doi.org/10.1137/0319022, https://arxiv.org/abs/https://doi.org/10.1137/0319022.
  * [24] J. Eugene C. Gartland, Scalings and limits of landau-de gennes models for liquid crystals: a comment on some recent analytical papers, Mathematical Modelling and Analysis, 23 (2018), pp. 414 – 432, https://doi.org/https://doi.org/10.3846/mma.2018.025.
  * [25] L. C. Evans, Partial Differential Equations, American Mathematical Society, Providence, Rhode Island, 1998.
  * [26] M. Farshbaf-Shaker, A penalty approach to optimal control of Allen–Cahn variational inequalities: MPEC-view, Numer. Func. Anal. Opt., (2012), p. 1321–1349.
  * [27] M. Farshbaf-Shaker, A relaxation approach to vector-valued Allen–Cahn MPEC problems, Appl. Math. Optim., (2015), pp. 325–351.
  * [28] J. W. Goodby, Handbook of Visual Display Technology (Editors: Chen, Janglin, Cranton, Wayne, Fihn, Mark), Springer, 2012, ch. Introduction to Defect Textures in Liquid Crystals, pp. 1290–1314.
  * [29] Y. Gu and N. L. Abbott, Observation of saturn-ring defects around solid microspheres in nematic liquid crystals, Phys. Rev. Lett., 85 (2000), pp. 4719–4722, https://doi.org/10.1103/PhysRevLett.85.4719.
  * [30] M. Heinkenschloss, The numerical solution of a control problem governed by a phase filed model, Optimization Methods and Software, 7 (1997), pp. 211–263, https://doi.org/10.1080/10556789708805656, https://doi.org/10.1080/10556789708805656, https://arxiv.org/abs/https://doi.org/10.1080/10556789708805656.
  * [31] M. Heinkenschloss and F. Tröltzsch, Analysis of the lagrange-sqp-newton method for the control of a phase field equation, Control Cybernet., 28 (1999), pp. 178–211.
  * [32] J. Heo, J.-W. Huh, and T.-H. Yoon, Fast-switching initially-transparent liquid crystal light shutter with crossed patterned electrodes, AIP Advances, 5 (2015), 047118, pp. –, https://doi.org/10.1063/1.4918277.
  * [33] K.-H. Hoffman and L. Jiang, Optimal control of a phase field model for solidification, Numerical Functional Analysis and Optimization, 13 (1992), pp. 11–27, https://doi.org/10.1080/01630569208816458, https://doi.org/10.1080/01630569208816458, https://arxiv.org/abs/https://doi.org/10.1080/01630569208816458.
  * [34] J. Hoogboom, J. A. Elemans, A. E. Rowan, T. H. Rasing, and R. J. Nolte, The development of self-assembled liquid crystal display alignment layers, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 365 (2007), pp. 1553–1576, https://doi.org/10.1098/rsta.2007.2031.
  * [35] M. Humar and I. Muševič, 3D microlasers from self-assembled cholesteric liquid-crystal microdroplets, Opt. Express, 18 (2010), pp. 26995–27003, https://doi.org/10.1364/OE.18.026995, http://www.opticsexpress.org/abstract.cfm?URI=oe-18-26-26995.
  * [36] S. Kralj and A. Majumdar, Order reconstruction patterns in nematic liquid crystal wells, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 470 (2014), https://doi.org/10.1098/rspa.2014.0276, http://rspa.royalsocietypublishing.org/content/470/2169/20140276, https://arxiv.org/abs/http://rspa.royalsocietypublishing.org/content/470/2169/20140276.full.pdf.
  * [37] J. P. Lagerwall and G. Scalia, A new era for liquid crystal research: Applications of liquid crystals in soft matter nano-, bio- and microtechnology, Current Applied Physics, 12 (2012), pp. 1387 – 1412, https://doi.org/https://doi.org/10.1016/j.cap.2012.03.019, http://www.sciencedirect.com/science/article/pii/S1567173912001113.
  * [38] G.-D. Lee, J. Anderson, and P. J. Bos, Fast Q-tensor method for modeling liquid crystal director configurations with defects, Applied Physics Letters, 81 (2002), pp. 3951–3953, https://doi.org/10.1063/1.1523157, https://doi.org/10.1063/1.1523157, https://arxiv.org/abs/https://doi.org/10.1063/1.1523157.
  * [39] F.-H. Lin and C. Liu, Static and dynamic theories of liquid crystals, Journal of Partial Differential Equations, 14 (2001), pp. 289–330.
  * [40] T. Lopez-Leon and A. Fernandez-Nieves, Drops and shells of liquid crystal, Colloid and Polymer Science, 289 (2011), pp. 345–359, https://doi.org/10.1007/s00396-010-2367-7.
  * [41] A. Majumdar, Equilibrium order parameters of nematic liquid crystals in the landau-de gennes theory, European Journal of Applied Mathematics, 21 (2010), pp. 181–203, https://doi.org/10.1017/S0956792509990210.
  * [42] A. Majumdar and A. Zarnescu, Landau-de gennes theory of nematic liquid crystals: the oseen-frank limit and beyond, Archive for rational mechanics and analysis, 196 (2010), pp. 227–280.
  * [43] H. Mori, J. Eugene C. Gartland, J. R. Kelly, and P. J. Bos, Multidimensional director modeling using the Q-tensor representation in a liquid crystal cell and its application to the $\pi$-cell with patterned electrodes, Japanese Journal of Applied Physics, 38 (1999), p. 135, http://stacks.iop.org/1347-4065/38/i=1R/a=135.
  * [44] N. J. Mottram and C. J. P. Newton, Introduction to Q-tensor theory, ArXiv e-prints, (2014), https://arxiv.org/abs/1409.3542.
  * [45] I. Muševič, M. Škarabot, U. Tkalec, M. Ravnik, and S. Žumer, Two-dimensional nematic colloidal crystals self-assembled by topological defects, Science, 313 (2006), pp. 954–958, https://doi.org/10.1126/science.1129660, http://www.sciencemag.org/content/313/5789/954.abstract, https://arxiv.org/abs/http://www.sciencemag.org/content/313/5789/954.full.pdf.
  * [46] I. Muševič and S. Žumer, Liquid crystals: Maximizing memory, Nature Materials, 10 (2011), pp. 266–268.
  * [47] M. Ravnik and S. Žumer, Landau-degennes modelling of nematic liquid crystal colloids, Liquid Crystals, 36 (2009), pp. 1201–1214, https://doi.org/10.1080/02678290903056095, https://doi.org/10.1080/02678290903056095, https://arxiv.org/abs/https://doi.org/10.1080/02678290903056095.
  * [48] A. Rešetič, J. Milavec, B. Zupančič, V. Domenici, and B. Zalar, Polymer-dispersed liquid crystal elastomers, Nature Communications, 7 (2016), p. 13140, https://doi.org/10.1038/ncomms13140.
  * [49] J. Schöberl, C++11 implementation of finite elements in NGSolve, Tech. Report ASC-2014-30, Institute for Analysis and Scientific Computing, September 2014, http://www.asc.tuwien.ac.at/~schoeberl/wiki/publications/ngs-cpp11.pdf.
  * [50] N. Schopohl and T. Sluckin, Defect core structure in nematic liquid crystals, Physical review letters, 59 (1987), p. 2582.
  * [51] A. A. Shah, H. Kang, K. L. Kohlstedt, K. H. Ahn, S. C. Glotzer, C. W. Monroe, and M. J. Solomon, Self-assembly: Liquid crystal order in colloidal suspensions of spheroidal particles by direct current electric field assembly (small 10/2012), Small, 8 (2012), pp. 1551–1562, https://doi.org/10.1002/smll.201290056.
  * [52] J. Shen and X. Yang, Numerical approximations of Allen-Cahn and Cahn-Hilliard equations, Discrete Contin. Dyn. Syst., 28 (2010), pp. 1669 – 1691.
  * [53] J. Sun, H. Wang, L. Wang, H. Cao, H. Xie, X. Luo, J. Xiao, H. Ding, Z. Yang, and H. Yang, Preparation and thermo-optical characteristics of a smart polymer-stabilized liquid crystal thin film based on smectic A-chiral nematic phase transition, Smart Materials and Structures, 23 (2014), p. 125038, http://stacks.iop.org/0964-1726/23/i=12/a=125038.
  * [54] F. Tröltzsch, Optimal Control of Partial Differential Equations, Graduate Studies in Mathematics, American Mathematical Society, April 2010.
  * [55] S. Čopar, U. Tkalec, I. Muševič, and S. Žumer, Knot theory realizations in nematic colloids, Proceedings of the National Academy of Sciences, 112 (2015), pp. 1675–1680, https://doi.org/10.1073/pnas.1417178112, http://www.pnas.org/content/112/6/1675.abstract, https://arxiv.org/abs/http://www.pnas.org/content/112/6/1675.full.pdf.
  * [56] E. G. Virga, Variational Theories for Liquid Crystals, vol. 8, Chapman and Hall, London, 1st ed., 1994.
  * [57] M. Wang, L. He, S. Zorba, and Y. Yin, Magnetically actuated liquid crystals, Nano Letters, 14 (2014), pp. 3966–3971, https://doi.org/10.1021/nl501302s, http://dx.doi.org/10.1021/nl501302s, https://arxiv.org/abs/http://dx.doi.org/10.1021/nl501302s.  PMID: 24914876.
  * [58] T. H. Ware, M. E. McConney, J. J. Wie, V. P. Tondiglia, and T. J. White, Voxelated liquid crystal elastomers, Science, 347 (2015), pp. 982–984, https://doi.org/10.1126/science.1261019, http://www.sciencemag.org/content/347/6225/982.abstract, https://arxiv.org/abs/http://www.sciencemag.org/content/347/6225/982.full.pdf.
  * [59] J. K. Whitmer, X. Wang, F. Mondiot, D. S. Miller, N. L. Abbott, and J. J. de Pablo, Nematic-field-driven positioning of particles in liquid crystal droplets, Phys. Rev. Lett., 111 (2013), p. 227801, https://doi.org/10.1103/PhysRevLett.111.227801.
  * [60] R. Zhang, A. Mozaffari, and J. J. de Pablo, Logic operations with active topological defects, Science Advances, 8 (2022), p. eabg9060, https://doi.org/10.1126/sciadv.abg9060, https://www.science.org/doi/abs/10.1126/sciadv.abg9060, https://arxiv.org/abs/https://www.science.org/doi/pdf/10.1126/sciadv.abg9060.
  * [61] J. Zhao and Q. Wang, Semi-discrete energy-stable schemes for a tensor-based hydrodynamic model of nematic liquid crystal flows, Journal of Scientific Computing, 68 (2016), pp. 1241–1266, https://doi.org/10.1007/s10915-016-0177-x, https://doi.org/10.1007/s10915-016-0177-x.
  * [62] J. Zhao, X. Yang, J. Shen, and Q. Wang, A decoupled energy stable scheme for a hydrodynamic phase-field model of mixtures of nematic liquid crystals and viscous fluids, Journal of Computational Physics, 305 (2016), pp. 539 – 556, https://doi.org/https://doi.org/10.1016/j.jcp.2015.09.044, http://www.sciencedirect.com/science/article/pii/S0021999115006439.
  * [63] W. Zhu, M. Shelley, and P. Palffy-Muhoray, Modeling and simulation of liquid-crystal elastomers, Phys. Rev. E, 83 (2011), p. 051703, https://doi.org/10.1103/PhysRevE.83.051703.