$\displaystyle=B_{2}({\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},\Phi_{h})-B_{2}({\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},\Theta_{h},\Phi_{h})=B_{2}({\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},\Phi_{h})$
$\displaystyle(\text{resp.
}B_{3}(\textrm{I}_{h}\Psi_{\mathbf{g}},{\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi},\Phi_{h})-2B_{3}({\rm{I}}_{h}\Psi_{\mathbf{g}},{\rm{I}}_{h}\widetilde{\Psi},\Theta_{h},\Phi_{h})+B_{3}(\textrm{I}_{h}\Psi_{\mathbf{g}},\Theta_{h},\Theta_{h},\Phi_{h})$
$\displaystyle=B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},{\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},\Phi_{h})-B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},\Theta_{h},\Phi_{h})=B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},\Phi_{h}))$

and achieved by utilizing the linearity and symmetry of
$B_{2}(\cdot,\cdot,\cdot)$ (resp. $B_{3}(\cdot,\cdot,\cdot,\cdot)$) in first
two variables (resp. second and third variables). Lemma 3.9$(i)$ (resp.
$(ii)$) for $\boldsymbol{\eta}:=\widetilde{\Psi}$ (resp.
$\boldsymbol{\xi}:=\boldsymbol{\eta}:=\widetilde{\Psi}$) shows

$\displaystyle T_{1}:=A({\rm
I}_{h}\widetilde{\Psi}-\widetilde{\Psi},\Phi_{h})+B_{1}({\rm
I}_{h}\widetilde{\Psi}-\widetilde{\Psi},\Phi_{h})\lesssim h+\ell^{-1}h^{2}.$
$\displaystyle(\text{resp.
}T_{2}:=B_{2}({\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi},\Phi_{h})-B_{2}(\widetilde{\Psi},\widetilde{\Psi},\Phi_{h})\lesssim\ell^{-1}h^{2}.)$

For $\mathbf{e}:={\rm I}_{h}\widetilde{\Psi}-\Theta_{h},$ Lemma 3.7$(ii)$ with
${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}=1$
shows

$T_{3}:=B_{2}(\textrm{I}_{h}\widetilde{\Psi}-{\Theta}_{h},\textrm{I}_{h}\widetilde{\Psi}-{\Theta}_{h},\Phi_{h})\lesssim\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\widetilde{\Psi}-{\Theta}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}=\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|{\mathbf{e}}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}.$

The estimate for $T_{4}$ (resp. $T_{5}$) follows from the second (resp. third)
inequality of Lemma 3.9$(iii)$ for
$\boldsymbol{\xi}:=\boldsymbol{\eta}:=\widetilde{\Psi}$ (resp.
$\boldsymbol{\eta}:={\rm{I}}_{h}\widetilde{\Psi}$), and Lemma 3.8.

$\displaystyle
T_{4}:=B_{3}({\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi},{\rm
I}_{h}\widetilde{\Psi},\Phi_{h})-B_{3}(\widetilde{\Psi},\widetilde{\Psi},\widetilde{\Psi},\Phi_{h})\lesssim\ell^{-1}h^{2}.$
$\displaystyle(\text{resp.
}T_{5}:=2B_{3}({\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi},{\rm
I}_{h}\widetilde{\Psi},\Phi_{h})-3B_{3}({\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi},\Theta_{h},\Phi_{h})+B_{3}(\Theta_{h},\Theta_{h},\Theta_{h},\Phi_{h})\lesssim\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}({|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+1).)$

Lemma 3.9$(ii)$ (resp. $(iii)$) with $\boldsymbol{\xi}:={\Psi}_{\mathbf{g}}$,
$\boldsymbol{\eta}:=\widetilde{\Psi}$ and (3.7) shows

$\displaystyle\frac{1}{2}T_{6}:=B_{2}({\rm{I}}_{h}{\Psi}_{\mathbf{g}},\textrm{I}_{h}\widetilde{\Psi},\Phi_{h})-B_{2}({\Psi}_{\mathbf{g}},\widetilde{\Psi},\Phi_{h})\lesssim\ell^{-1}h^{2}.$
$\displaystyle(\text{resp.
}\frac{1}{3}T_{7}:=B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},{\rm{I}}_{h}\widetilde{\Psi},{\rm{I}}_{h}\widetilde{\Psi},\Phi_{h})-B_{3}({\Psi}_{\mathbf{g}},\widetilde{\Psi},\widetilde{\Psi},\Phi_{h})\lesssim\ell^{-1}h^{2}.)$

Apply Lemma 3.9$(iii)$ with $\boldsymbol{\xi}:=\widetilde{\Psi}$,
$\boldsymbol{\eta}:={\Psi}_{\mathbf{g}}$ and (3.7) to obtain

$\displaystyle\frac{1}{3}T_{8}:$
$\displaystyle=B_{3}({\rm{I}}_{h}{\Psi}_{\mathbf{g}},{\rm{I}}_{h}{\Psi}_{\mathbf{g}},\textrm{I}_{h}\widetilde{\Psi},\Phi_{h})-B_{3}({\Psi}_{\mathbf{g}},{\Psi}_{\mathbf{g}},\widetilde{\Psi},\Phi_{h})\lesssim\ell^{-1}h^{2}.$

Lemmas 3.7(iii), 3.8 and (3.7) leads to

$\displaystyle\frac{1}{3}T_{9}:=B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},{\rm{I}}_{h}\widetilde{\Psi}-\Theta_{h},\Phi_{h})\lesssim\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}.$

The definition of $N(\cdot;\cdot)$ in (3.2), a re-arrangement of terms, Lemma
3.9$(i)$-$(iii)$ with
$\boldsymbol{\xi}:=\boldsymbol{\eta}:={\Psi}_{\mathbf{g}}$, and (3.7) yields

$\displaystyle T_{10}$
$\displaystyle:=N(\textrm{I}_{h}{\Psi}_{\mathbf{g}};\Phi_{h})-N({\Psi}_{\mathbf{g}};\Phi_{h})=\big{(}A(\textrm{I}_{h}{\Psi}_{\mathbf{g}}-{\Psi}_{\mathbf{g}},\Phi_{h})+B_{1}(\textrm{I}_{h}{\Psi}_{\mathbf{g}}-{\Psi}_{\mathbf{g}},\Phi_{h})\big{)}+\big{(}B_{2}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},\textrm{I}_{h}{\Psi}_{\mathbf{g}},\Phi_{h})$
$\displaystyle\quad-
B_{2}({\Psi}_{\mathbf{g}},{\Psi}_{\mathbf{g}},\Phi_{h})\big{)}+\big{(}B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},\textrm{I}_{h}{\Psi}_{\mathbf{g}},\textrm{I}_{h}{\Psi}_{\mathbf{g}},\Phi_{h})-B_{3}({\Psi}_{\mathbf{g}},{\Psi}_{\mathbf{g}},{\Psi}_{\mathbf{g}},\Phi_{h})\big{)}\lesssim
h+\ell^{-1}h^{2}.$

Substitute the estimates for $T_{1},\ldots,T_{10}$ in (3.4) and utilize
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
R(h)$ to obtain

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|{{\rm
I}_{h}\widetilde{\Psi}}-\mu_{h}(\Theta_{h})|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
C_{3}(h+\ell^{-1}h^{2}+\ell^{-1}R(h)^{2}(R(h)+1)),$

where the constant $C_{3}$ is independent of $h$ and $\ell.$ Assume
$h\leq\ell^{1+\varsigma}$ with $\varsigma>0$ so that $\ell^{-1}h\leq
h^{\frac{\varsigma}{1+\varsigma}}.$ Choose $R(h)=2C_{3}h.$ For
$h<h_{4}:=\min({h_{3},h_{2}})$ with
$h_{3}^{\frac{\varsigma}{1+\varsigma}}<\frac{1}{2(1+4C^{3}_{2})^{2}}<\frac{1}{2}$,

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\widetilde{\Psi}-\mu_{h}(\Theta_{h})|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
C_{3}h\left(1+h^{\frac{\varsigma}{1+\varsigma}}(1+4C_{3}^{2})+8C_{3}^{3}hh^{\frac{\varsigma}{1+\varsigma}}\right)\leq{C_{3}h\bigg{(}1+\frac{1}{2}+\frac{1}{2}\frac{8C_{3}^{3}h}{(1+4C_{3}^{2})^{2}}\bigg{)}}.$

Since $h<h_{3}<\frac{1}{2^{\frac{1+\varsigma}{\varsigma}}}<1$ and
$\frac{8C_{3}^{3}}{(1+4C_{3}^{2})^{2}}<1,$
${|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\widetilde{\Psi}-\mu_{h}(\Theta_{h})|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
2C_{3}h=R(h).$

Step 3 ($\mu_{h}$ is a contraction). Let
$\Theta_{1},\Theta_{2}\in\mathbb{B}_{R(h)}({\rm I}_{h}\widetilde{\Psi})$ and
set $\mathbf{e}_{1}:=\textrm{I}_{h}\widetilde{\Psi}-\Theta_{1},$
$\mathbf{e}_{2}:=\textrm{I}_{h}\widetilde{\Psi}-\Theta_{2},$ and
$\mathbf{e}_{3}:=\Theta_{1}-\Theta_{2}.$ The linearity of $\langle
D\widetilde{N}({{\rm I}_{h}\widetilde{\Psi}})\cdot,\cdot\rangle,$ (3.4), and a
re-arrangement of terms yields that for all $\Phi_{h}\in\mathbf{V}_{h},$

$\displaystyle\langle D\widetilde{N}({{\rm
I}_{h}\widetilde{\Psi}})(\mu_{h}(\Theta_{1})-\mu_{h}(\Theta_{2})),\Phi_{h}\rangle=\langle
D\widetilde{N}({{\rm
I}_{h}\widetilde{\Psi}})\mu_{h}(\Theta_{1}),\Phi_{h}\rangle-\langle
D\widetilde{N}({{\rm
I}_{h}\widetilde{\Psi}})\mu_{h}(\Theta_{2}),\Phi_{h}\rangle$
$\displaystyle=\big{(}B_{2}(\mathbf{e}_{1},\mathbf{e}_{3},\Phi_{h})+B_{2}(\mathbf{e}_{2},\mathbf{e}_{3},\Phi_{h})\big{)}+3\big{(}B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},\mathbf{e}_{1},\mathbf{e}_{3},\Phi_{h})+B_{3}(\textrm{I}_{h}{\Psi}_{\mathbf{g}},\mathbf{e}_{2},\mathbf{e}_{3},\Phi_{h})\big{)}$
$\displaystyle\quad+\big{(}3B_{3}(\textrm{I}_{h}\widetilde{\Psi},\textrm{I}_{h}\widetilde{\Psi},\Theta_{1},\Phi_{h})-B_{3}(\Theta_{1},\Theta_{1},\Theta_{1},\Phi_{h})-3B_{3}(\textrm{I}_{h}\widetilde{\Psi},\textrm{I}_{h}\widetilde{\Psi},\Theta_{2},\Phi_{h})+B_{3}(\Theta_{2},\Theta_{2},\Theta_{2},\Phi_{h})\big{)}$
$\displaystyle=:T_{1}^{\prime}+T_{2}^{\prime}+T_{3}^{\prime}.$ (3.20)

Note that
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{1}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
2C_{3}h$,
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{2}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
2C_{3}h$. This combined with Lemma 3.7$(ii)$-$(iii)$ and (3.7) implies

$\displaystyle
T_{1}^{\prime}\lesssim\ell^{-1}h{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{3}|\kern-1.07639pt|\kern-1.07639pt|}_{1}{|\kern-1.07639pt|\kern-1.07639pt|\Phi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\text{
and
}T_{2}^{\prime}\lesssim\ell^{-1}h{|\kern-1.07639pt|\kern-1.07639pt|{\mathbf{e}_{3}}|\kern-1.07639pt|\kern-1.07639pt|}_{1}{|\kern-1.07639pt|\kern-1.07639pt|\Phi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}.$

The estimation of the term $T_{3}^{\prime}$ utilizes fourth inequality of
Lemma 3.9$(iii)$ with $\boldsymbol{\eta}:=\textrm{I}_{h}\widetilde{\Psi},$
Lemma 3.8, and
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{1}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
2C_{3}h$,
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{2}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
2C_{3}h$.

$\displaystyle
T_{3}^{\prime}\lesssim\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{3}|\kern-1.07639pt|\kern-1.07639pt|}_{1}({|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{1}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{2}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}+({|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{1}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{2}|\kern-1.07639pt|\kern-1.07639pt|}_{1}){|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\widetilde{\Psi}|\kern-1.07639pt|\kern-1.07639pt|}_{1}){|\kern-1.07639pt|\kern-1.07639pt|\Phi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim\ell^{-1}h(1+h){|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{3}|\kern-1.07639pt|\kern-1.07639pt|}_{1}{|\kern-1.07639pt|\kern-1.07639pt|\Phi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}.$

Substitute the above displayed three estimates for
$T_{1}^{\prime},\cdots,T_{3}^{\prime}$ in (3.4). This plus the discrete inf-
sup condition in Lemma 3.11$(ii)$ implies that there exists a
$\Phi_{h}\in\mathbf{V}_{h}$ with
${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}=1$
such that

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\mu_{h}(\Theta_{1})-\mu_{h}(\Theta_{2})|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
4\beta^{-1}\langle DN({{\rm
I}_{h}\widetilde{\Psi}})(\mu_{h}(\Theta_{1})-\mu_{h}(\Theta_{2})),\Phi_{h}\rangle\lesssim\ell^{-1}h(1+h){|\kern-1.07639pt|\kern-1.07639pt|\Theta_{1}-\Theta_{2}|\kern-1.07639pt|\kern-1.07639pt|}_{1}.$

The assumption $h\leq\ell^{1+\varsigma}$ allows

${|\kern-1.07639pt|\kern-1.07639pt|\mu_{h}(\Theta_{1})-\mu_{h}(\Theta_{2})|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim
h^{\frac{\varsigma}{1+\varsigma}}{|\kern-1.07639pt|\kern-1.07639pt|\Theta_{1}-\Theta_{2}|\kern-1.07639pt|\kern-1.07639pt|}_{1},$

where the hidden constant in $"\lesssim"$ depends on
${|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\Psi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}$,
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{3}{2}}$,
$\beta$, $C_{I}$, $C_{3}$.

Step 4 (Existence and uniqueness). The nonlinear map $\mu_{h}$ is well-
defined, continuous and maps a closed convex subset
$\mathbb{B}_{R}(\textrm{I}_{h}\widetilde{\Psi})$ of a Hilbert space
$\mathbf{V}_{h}$ to itself. Therefore, Brouwer’s fixed point theorem [23] and
Step 3 imply the existence and uniqueness of the fixed point, say
$\widetilde{\Psi}_{h}$ in the ball
$\mathbb{B}_{R}(\textrm{I}_{h}\widetilde{\Psi})$. A triangle inequality,
${|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\widetilde{\Psi}-\widetilde{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq
2C_{3}h$ and Lemma 3.8 show
${|\kern-1.07639pt|\kern-1.07639pt|{\widetilde{\Psi}}-{\widetilde{\Psi}}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq(C_{I}+2C_{3})h$.
This concludes the proof. ∎

###### Proof of Theorem 3.3.

Recall that
$\widetilde{\Psi}_{h}={\Psi}_{h}-\textrm{I}_{h}{\Psi}_{\mathbf{g}}$, where
${\Psi}_{h}$ satisfies the discrete non-linear system (3.3). Theorem 3.13
shows the existence and local uniqueness of the discrete solution
${\Psi}_{h}$. Moreover,
$\widetilde{\Psi}_{h}\in\mathbb{B}_{R}(\textrm{I}_{h}\widetilde{\Psi})$ yields
that ${\Psi}_{h}\in\mathbb{B}_{R}(\textrm{I}_{h}{\Psi}).$ This, Lemma 3.8 and
triangle inequality lead to

${|\kern-1.07639pt|\kern-1.07639pt|{\Psi}-{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim{|\kern-1.07639pt|\kern-1.07639pt|{\Psi}-\textrm{I}_{h}{\Psi}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}{\Psi}-{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim
h.\qed$

The best approximation result presented in Theorem 3.4 is established next.
The technique used in [33, Theorem 3.3] yields a best approximation result in
$\widehat{\mathbf{X}}_{h}:=\\{\Theta_{h}\in\mathbf{X}_{h}|\,\Theta_{h}=\mathbf{g}_{h}\text{
on }\partial\Omega\\}\subset\mathbf{X}_{h}.$ In this article, we use elliptic
projection of $\Psi$ onto $\mathbf{X}_{h}$ to establish the best approximation
result in $\mathbf{X}_{h}$.

###### Proof of Theorem 3.4.

Step 1 (Best approximation on
$\widehat{\mathbf{X}}_{h}:=\\{\Theta_{h}\in\mathbf{X}_{h}|\,\Theta_{h}=\mathbf{g}_{h}\text{
on }\partial\Omega\\}$). Set $\mathbf{e}:=\Psi-\Psi_{h}.$ The Taylor series
expansion of ${N}(\cdot;\cdot)$ around $\Psi$ imply

$\displaystyle{N}(\Psi_{h};\Phi_{h})$
$\displaystyle={N}(\Psi;\Phi_{h})-\langle{DN}(\Psi)\mathbf{e},\Phi_{h}\rangle+\frac{1}{2}\langle{D}^{2}N(\Psi)(\mathbf{e})\mathbf{e},\Phi_{h}\rangle-\frac{1}{6}\langle{D}^{3}N(\Psi)(\mathbf{e})(\mathbf{e})\mathbf{e},\Phi_{h}\rangle.$

For $\widehat{\Psi}_{h}\in\widehat{\mathbf{X}}_{h}$, use
$\mathbf{e}=(\Psi-\widehat{\Psi}_{h})+(\widehat{\Psi}_{h}-\Psi_{h})$, the
linearity of $\langle DN(\Psi)\cdot,\cdot\rangle$, and (3.2), (3.3) to obtain

$\displaystyle\langle
DN(\Psi)(\widehat{\Psi}_{h}-\Psi_{h}),\Phi_{h}\rangle=\langle
DN(\Psi)(\widehat{\Psi}_{h}-\Psi),\Phi_{h}\rangle+\frac{1}{2}\langle{D}^{2}N(\Psi)(\mathbf{e})\mathbf{e},\Phi_{h}\rangle-\frac{1}{6}\langle{D}^{3}N(\Psi)(\mathbf{e})(\mathbf{e})\mathbf{e},\Phi_{h}\rangle.$
(3.21)

The identities $\langle
DN(\Psi)\mathbf{e},\Phi_{h}\rangle=A(\mathbf{e},\Phi_{h})+B_{1}(\mathbf{e},\Phi_{h})+2B_{2}(\Psi,\mathbf{e},\Phi_{h})+3B_{3}(\Psi,\Psi,\mathbf{e},\Phi_{h}),$
$\langle{D}^{2}N(\Psi)(\mathbf{e})\mathbf{e},\Phi_{h}\rangle=2B_{2}(\mathbf{e},\mathbf{e},\Phi_{h})+6B_{3}(\mathbf{e},\mathbf{e},\Psi,\Phi_{h})$,
$\langle{D}^{3}N(\Psi)(\mathbf{e})(\mathbf{e})\mathbf{e},\Phi_{h}\rangle=6B_{3}(\mathbf{e},\mathbf{e},\mathbf{e},\Phi_{h}),$
and Lemma 3.7 leads to

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|DN(\Psi)|\kern-1.07639pt|\kern-1.07639pt|}_{\mathbf{L}^{2}}\lesssim(1+\ell^{-1}),\,{|\kern-1.07639pt|\kern-1.07639pt|{D}^{2}N(\Psi)|\kern-1.07639pt|\kern-1.07639pt|}_{\mathbf{L}^{2}}\lesssim\ell^{-1}\text{
and
}{|\kern-1.07639pt|\kern-1.07639pt|{D}^{3}N(\Psi)|\kern-1.07639pt|\kern-1.07639pt|}_{\mathbf{L}^{2}}\lesssim\ell^{-1}.$
(3.22)

A triangle inequality yields
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\Psi}_{h}-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}$.
Since $(\widehat{\Psi}_{h}-{\Psi}_{h})|_{\partial\Omega}=0,$ the discrete inf-
sup condition from Remark 3.12 with
$\Theta_{h}:=\widehat{\Psi}_{h}-{\Psi}_{h}\in\mathbf{V}_{h}\subset\mathbf{V}$,
(3.21) and (3.22) yields

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}$
$\displaystyle\lesssim{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+\langle
DN(\Psi)(\widehat{\Psi}_{h}-{\Psi}_{h}),\Phi_{h}\rangle\lesssim(1+\ell^{-1}){|\kern-1.07639pt|\kern-1.07639pt|\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}(1+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}).$
(3.23)

For a sufficiently small choice of the discretization parameter
$h=O(\ell^{1+\varsigma})$ with $\varsigma>0$, use
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim
h$ from Theorem 3.3, $\ell^{-1}h\leq h^{\frac{\varsigma}{1+\varsigma}},$ and
$h<1$ (this holds for $h<h_{4}$) in (3.23) to obtain

$\displaystyle
C_{4}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq(1+\ell^{-1}){|\kern-1.07639pt|\kern-1.07639pt|\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+h^{\frac{\varsigma}{1+\varsigma}}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1},$

where the constant $C_{4}$ depends on $\left\lvert
c\right\rvert,{|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$,
$C_{I},C_{3}$ and $\beta$. A sufficiently small choice of
$h<h_{6}:=\min(h_{5},h_{4})$ with
$h_{5}^{\frac{\varsigma}{1+\varsigma}}=\frac{C_{4}}{2}$ leads to

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim(1+\ell^{-1}){|\kern-1.07639pt|\kern-1.07639pt|\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}.$
(3.24)

Step 2 (Best approximation on $\mathbf{X}_{h}$). Let $R_{h}\Psi$ be the
elliptic projection ($\mathbf{H}^{1}(\Omega)$ projection) of $\Psi$ onto
$\mathbf{X}_{h}$ defined by

$\displaystyle A(\Psi-R_{h}\Psi,\Phi_{h})+(\Psi-R_{h}\Psi,\Phi_{h})=0\text{
for all }\Phi_{h}\in\mathbf{X}_{h}.$

Then it holds

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\Psi-
R_{h}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{1}=\inf_{{\Psi}^{*}_{h}\in\mathbf{X}_{h}}{|\kern-1.07639pt|\kern-1.07639pt|\Psi-{\Psi}^{*}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}.$
(3.25)

Let $\mathbf{g}_{h}^{*}=R_{h}\Psi|_{\partial\Omega}.$ Choose
$\boldsymbol{\eta}_{h}\in\mathbf{X}_{h}$ such that
$\boldsymbol{\eta}_{h}|_{\partial\Omega}=\mathbf{g}_{h}-\mathbf{g}_{h}^{*}$
and
${|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\eta}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}_{h}-\mathbf{g}_{h}^{*}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}$.
For
$\widehat{\Psi}_{h}=\boldsymbol{\eta}_{h}+R_{h}\Psi\in\widehat{\mathbf{X}}_{h}$,
this plus a triangle inequality lead to

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq{|\kern-1.07639pt|\kern-1.07639pt|\Psi-
R_{h}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|R_{h}\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq{|\kern-1.07639pt|\kern-1.07639pt|\Psi-
R_{h}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}_{h}-\mathbf{g}_{h}^{*}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}.$

A triangle inequality
${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}_{h}-\mathbf{g}_{h}^{*}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}\leq{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}^{*}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}$,
trace inequality

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}^{*}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}\leq{|\kern-1.07639pt|\kern-1.07639pt|{\Psi}-R_{h}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{1},$
(3.26)

and the definition of elliptic projection in (3.25) show

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\widehat{\Psi}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\leq\inf_{{\Psi}^{*}_{h}\in\mathbf{X}_{h}}{|\kern-1.07639pt|\kern-1.07639pt|\Psi-{\Psi}^{*}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}.$

This combined with the best approximation result (3.24) obtained in Step 1
leads to the desired estimate. ∎

######  Lemma 3.14 (Estimate for boundary term).

[19] Let $\widehat{\Psi}\in\mathbf{X}$ with
$\widehat{\Psi}|_{\partial\Omega}=\widehat{\mathbf{g}}$ be given. Let
$R_{h}\widehat{\Psi}\in\mathbf{X}_{h}$ be the
$\mathbf{H}^{1}(\Omega)$-projection of $\widehat{\Psi}$ onto $\mathbf{X}_{h}$
and $\widehat{\mathbf{g}}_{h}:=(R_{h}\widehat{\Psi})|_{\partial\Omega}.$ Then

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathbf{g}}-\widehat{\mathbf{g}}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\partial\Omega}\leq
Ch^{\frac{1}{2}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\Psi}-R_{h}\widehat{\Psi}|\kern-1.07639pt|\kern-1.07639pt|}_{1}.$

######  Remark 3.15 (Best approximation on $\mathbf{X}_{h}$).

A triangle inequality and an inverse inequality yield

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}\leq{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-{\mathbf{g}}^{*}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}^{*}_{h}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}\lesssim{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}^{*}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}+h^{-\frac{1}{2}}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}_{h}-\mathbf{g}^{*}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\partial\Omega}.$

The trace inequality (3.26), a triangle inequality and Lemma 3.14 applied to
the above inequality yields

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}\lesssim{|\kern-1.07639pt|\kern-1.07639pt|\Psi-
R_{h}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{1}+h^{-\frac{1}{2}}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\partial\Omega}.$

This combined with the best approximation result obtained in Theorem 3.4 and
(3.25) leads to

${|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim(1+\ell^{-1})(\min_{\Psi^{*}_{h}\in\mathbf{X}_{h}}{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi^{*}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+h^{-\frac{1}{2}}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\partial\Omega}).\qed$

The proof of $\mathbf{L}^{2}$ norm error in Theorem 3.5 is established next.

###### Proof of Theorem 3.5.

Set $G=\Psi-\Psi_{h}\in\mathbf{L}^{2}(\Omega)$. Consider the well-posed dual
linear problem that seeks
$\boldsymbol{\chi}\in\mathbf{H}^{2}(\Omega)\cap\mathbf{V}$ such that

$\displaystyle
A(\Phi,\boldsymbol{\chi})+B_{1}(\Phi,\boldsymbol{\chi})+2B_{2}(\Psi,\Phi,\boldsymbol{\chi})+3B_{3}(\Psi,\Psi,\Phi,\boldsymbol{\chi})=(G,\Phi)\text{
for all }\Phi\in\mathbf{V}$ (3.27)

and satisfies

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\lesssim(1+\ell^{-1}){|\kern-1.07639pt|\kern-1.07639pt|G|\kern-1.07639pt|\kern-1.07639pt|}_{0},$
(3.28)

where the hidden constant in $"\lesssim"$ depends on $\left\lvert
c\right\rvert,{|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$,
$\beta$ and the constants in Sobolev embedding results. For
$G:=(w_{1},w_{2},w_{3},w_{4})$ and
$\boldsymbol{\chi}:=(\chi_{1},\chi_{2},\chi_{3},\chi_{4})$, the strong form of
the dual linear problem (3.27) is defined as

$\left.\begin{split}&-\Delta\chi_{1}+\ell^{-1}\big{(}-\chi_{1}-\frac{c}{2}(u_{3}\chi_{3}-u_{4}\chi_{4})+\frac{1}{3}((3u_{1}^{2}+u_{2}^{2})\chi_{1}+2u_{1}u_{2}\chi_{2})\big{)}=w_{1},\\\
&-\Delta\chi_{2}+\ell^{-1}\big{(}-\chi_{2}-\frac{c}{2}(u_{3}\chi_{4}+u_{4}\chi_{3})+\frac{1}{3}((u_{1}^{2}+3u_{2}^{2})\chi_{2}+2u_{1}u_{2}\chi_{1})\big{)}=w_{2},\\\
&-\Delta\chi_{3}+\ell^{-1}\big{(}-\chi_{3}-\frac{c}{2}(u_{3}\chi_{1}+u_{1}\chi_{3}+u_{2}\chi_{4}+u_{4}\chi_{2})+\frac{1}{3}((3u_{3}^{2}+u_{4}^{2})\chi_{3}+2u_{3}u_{4}\chi_{4})\big{)}=w_{3},\\\
&-\Delta\chi_{4}+\ell^{-1}\big{(}-\chi_{4}-\frac{c}{2}(u_{2}\chi_{3}+u_{3}\chi_{2}-u_{1}\chi_{4}-u_{4}\chi_{1})+\frac{1}{3}((u_{3}^{2}+3u_{4}^{2})\chi_{4}+2u_{3}u_{4}\chi_{3})\big{)}=w_{4}.\end{split}\right\\}$
(3.29)

Let $\Psi^{*}\in\mathbf{X}$ (resp. $\Psi_{h}^{*}\in\mathbf{X}_{h}$) be
extension of $\mathbf{g}$ (resp. $\mathbf{g}_{h}$) such that
$\Psi^{*}|_{\partial\Omega}=\mathbf{g}$ (resp.
$\Psi^{*}_{h}|_{\partial\Omega}=\mathbf{g}_{h}$). Let
$\Psi^{0}:=\Psi-\Psi^{*},$ and $\Psi^{0}_{h}:=\Psi_{h}-\Psi_{h}^{*}.$ Set
$\Phi=\Psi^{0}-\Psi^{0}_{h}\in\mathbf{V}$ in (3.27) to obtain

$\displaystyle(\Psi-\Psi_{h},\Psi^{0}-\Psi_{h}^{0})=\langle{DN}(\Psi)(\Psi^{0}-\Psi_{h}^{0}),\boldsymbol{\chi}\rangle.$
(3.30)

Test (3.29) with $\Psi_{h}^{*}-\Psi^{*}$ and use integration by parts to
obtain

$\displaystyle(\Psi-\Psi_{h},\Psi_{h}^{*}-\Psi^{*})=\langle{DN}(\Psi)(\Psi_{h}^{*}-\Psi^{*}),\boldsymbol{\chi}\rangle+\int_{\partial\Omega}(\mathbf{g}-\mathbf{g}_{h})\cdot\frac{\partial\boldsymbol{\chi}}{\partial\nu}\,{\rm
ds}.$ (3.31)

The trace inequality
${|\kern-1.07639pt|\kern-1.07639pt|\frac{\partial\boldsymbol{\chi}}{\partial\nu}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}\lesssim{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}$
leads to

$\int_{\partial\Omega}(\mathbf{g}-\mathbf{g}_{h})\cdot\frac{\partial\boldsymbol{\chi}}{\partial\nu}\,{\rm
ds}\leq{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{-\frac{1}{2},\partial\Omega}{|\kern-1.07639pt|\kern-1.07639pt|\frac{\partial\boldsymbol{\chi}}{\partial\nu}|\kern-1.07639pt|\kern-1.07639pt|}_{\frac{1}{2},\partial\Omega}\lesssim{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{-\frac{1}{2},\partial\Omega}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}.$

This and $\Psi^{0}-\Psi^{0}_{h}=(\Psi-\Psi_{h})+(\Psi_{h}^{*}-\Psi^{*})$
applied to (3.30) leads to

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{0}\leq\langle{DN}(\Psi)(\Psi-\Psi_{h}),\boldsymbol{\chi}\rangle+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{-\frac{1}{2},\partial\Omega}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}.$
(3.32)

A term
$\langle{DN}(\Psi)(\Psi-\Psi_{h}),\textrm{I}_{h}\boldsymbol{\chi}\rangle$ with
$\textrm{I}_{h}\boldsymbol{\chi}\in\mathbf{V}_{h}\subset\mathbf{V}$ is added
and subtracted to the first term on the right hand side of (3.32), (3.2)-(3.3)
are utilized and simple manipulations are performed to arrive at

$\displaystyle\langle{DN}(\Psi)(\Psi-\Psi_{h}),\boldsymbol{\chi}\rangle=\langle{DN}(\Psi)(\Psi-\Psi_{h}),\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi}\rangle+\langle{DN}(\Psi)(\Psi-\Psi_{h}),\textrm{I}_{h}\boldsymbol{\chi}\rangle+N(\Psi_{h},{\rm
I}_{h}\boldsymbol{\chi})-N(\Psi,{\rm I}_{h}\boldsymbol{\chi})$
$\displaystyle=\langle{DN}(\Psi)(\Psi-\Psi_{h}),\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi}\rangle+B_{2}({\Psi}_{h}-{\Psi},{\Psi}_{h}-{\Psi},\textrm{I}_{h}\boldsymbol{\chi})+(2B_{3}(\Psi,\Psi,\Psi,{\rm
I}_{h}\boldsymbol{\chi})-3B_{3}(\Psi,\Psi,\Psi_{h},{\rm
I}_{h}\boldsymbol{\chi})$
$\displaystyle\quad+B_{3}(\Psi_{h},\Psi_{h},\Psi_{h},{\rm
I}_{h}\boldsymbol{\chi})):=T_{1}+T_{2}+T_{3}.$

Here the term $T_{2}$ is a re-grouping of

$2B_{2}({\Psi},{\Psi}-{\Psi}_{h},\textrm{I}_{h}\boldsymbol{\chi})+B_{2}({\Psi}_{h},{\Psi}_{h},\textrm{I}_{h}\boldsymbol{\chi})-B_{2}({\Psi},{\Psi},\textrm{I}_{h}\boldsymbol{\chi})=B_{2}({\Psi}_{h}-{\Psi},{\Psi}_{h},\textrm{I}_{h}\boldsymbol{\chi})-B_{2}({\Psi},{\Psi}_{h}-{\Psi},\textrm{I}_{h}\boldsymbol{\chi})$

and obtained by applying the linearity and symmetry of
$B_{2}(\cdot,\cdot,\cdot)$ in first two variables. The definition of
$\langle{DN}(\Psi)\cdot,\cdot\rangle$, Lemmas 3.7$(i)$-$(iii)$, 3.8 and
Theorem 3.3 show

$\displaystyle T_{1}:$
$\displaystyle=\langle{DN}(\Psi)(\Psi-\Psi_{h}),\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi}\rangle$
$\displaystyle=A(\Psi-\Psi_{h},\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi})+B_{1}(\Psi-\Psi_{h},\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi})+2B_{2}(\Psi,\Psi-\Psi_{h},\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi})+3B_{3}(\Psi,\Psi,\Psi-\Psi_{h},\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi})$
$\displaystyle\lesssim(1+\ell^{-1}){|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{1}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim(1+\ell^{-1})h^{2}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}.$

Lemmas 3.7$(ii)$, 3.8 and Theorem 3.3 yield

$\displaystyle
T_{2}:=B_{2}({\Psi}_{h}-{\Psi},{\Psi}_{h}-{\Psi},\textrm{I}_{h}\boldsymbol{\chi})\lesssim\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|{\Psi}_{h}-{\Psi}|\kern-1.07639pt|\kern-1.07639pt|}_{1}^{2}{|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{1}\lesssim\ell^{-1}h^{2}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}.$

Set $\mathbf{e}:=\Psi-\Psi_{h}$. Utilize the third inequality of Lemma
3.9$(iii)$ for $\boldsymbol{\eta}:=\Psi$, $\Theta_{h}:=\Psi_{h}$,
$\Phi_{h}:={\rm I}_{h}\boldsymbol{\chi}$, and then apply Lemma 3.8, Theorem
3.3 to estimate $T_{3}$.

$\displaystyle T_{3}$ $\displaystyle:=2B_{3}(\Psi,\Psi,\Psi,{\rm
I}_{h}\boldsymbol{\chi})-3B_{3}(\Psi,\Psi,\Psi_{h},{\rm
I}_{h}\boldsymbol{\chi})+B_{3}(\Psi_{h},\Psi_{h},\Psi_{h},{\rm
I}_{h}\boldsymbol{\chi})$
$\displaystyle\lesssim\ell^{-1}{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{1}({|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{1}){|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\lesssim\ell^{-1}h^{2}(h+1){|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}.$

The estimates of $T_{1},$ $T_{2}$ and $T_{3}$ in (3.32) and (3.28) yield

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{0}\lesssim(1+\ell^{-1})\big{(}h^{2}(1+\ell^{-1})+{|\kern-1.07639pt|\kern-1.07639pt|\mathbf{g}-\mathbf{g}_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{-\frac{1}{2},\partial\Omega}\big{)},$

where the constants suppressed in $"\lesssim"$ depends on $\left\lvert
c\right\rvert,{|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$,
$\beta$, $C_{I},C_{3}$ and the constants in Sobolev embedding results. This
completes the proof. ∎

######  Remark 3.16.

Integration by parts in (3.31) leads to the boundary term
$\int_{\partial\Omega}(\mathbf{g}-\mathbf{g}_{h})\cdot\frac{\partial\boldsymbol{\chi}}{\partial\nu}\,{\rm
ds}$, which gives a sub-optimal convergence rate in the $\mathbf{L}^{2}$ norm.
An optimal convergence rate $\mathcal{O}(h^{2})$ in $\mathbf{L}^{2}$ norm is
obtained using Nitsche’s method discussed below.∎

###  3.5 Nitsche’s method

Let $\mathcal{E}$ (resp. $\mathcal{E}(\Omega)$ or
$\mathcal{E}(\partial\Omega)$) denote the set of all (resp. interior or
boundary) edges in $\mathcal{T}.$ For Nitsche’s method, the finite element
space $\mathbf{X}_{h}$ associated with the triangulation $\mathcal{T}$ of the
convex polygonal domain $\Omega\subset\mathbb{R}^{2}$ into triangles is
endowed with the mesh dependent norm defined by
${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{h}:={|\kern-1.07639pt|\varphi_{1}|\kern-1.07639pt|}_{h}+{|\kern-1.07639pt|\varphi_{2}|\kern-1.07639pt|}_{h}+{|\kern-1.07639pt|\varphi_{3}|\kern-1.07639pt|}_{h}+{|\kern-1.07639pt|\varphi_{4}|\kern-1.07639pt|}_{h}$
for all
$\Phi_{h}=(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4})\in\mathbf{X}_{h}$,
where for all $v\in X_{h},$

$\displaystyle{|\kern-1.07639pt|v|\kern-1.07639pt|}^{2}_{h}:=\int_{\Omega}\left\lvert\nabla
v\right\rvert^{2}\,{\rm
dx}+\sum_{E\in\mathcal{E}(\partial\Omega)}\frac{\sigma}{h_{E}}\int_{E}v^{2}\,{\rm
ds}.$

Here $\sigma>0$ is the penalty parameter and $h_{E}$ denote length of an edge
$E$. Let $\nu_{T}$ denotes the unit outward normal along $\partial T$ of
$T\in\mathcal{T}$. The jump $[\varphi]_{E}$ of piecewise $H^{1}$ function
$\varphi,$ i.e, $\varphi\in H^{1}(\mathcal{T}):=\\{v\in L^{2}(\Omega)|\,\,v\in
H^{1}(T)\text{ for all }T\in\mathcal{T}\\},$ across $E\in\mathcal{E}$ is
defined by

$[\varphi]_{E}(x):=\begin{cases}v|_{T_{+}}(x)-v|_{T_{-}}(x)&\text{ for }x\in
E=\partial T_{+}\cap\partial T_{-}\in\mathcal{E}(\Omega),\\\ v(x)&\text{ for
}x\in E\in\mathcal{E}(\partial\Omega),\end{cases}$

where for the interior edge $E=\partial T_{+}\cap\partial
T_{-}\in\mathcal{E}(\Omega)$ with unit normal $\nu_{E}$ of fixed orientation,
the adjacent triangles $T_{\pm}\in\mathcal{T}$ are in an order such that
$\nu_{E}=\nu_{T_{+}}|_{E}=-\nu_{T_{-}}|_{E}$. The discrete formulation for
Nitsche’s method that corresponds to (3.2) seeks
$\Psi_{h}\\!\in\\!\mathbf{X}_{h}$ such that for all
$\Phi_{h}\in\mathbf{X}_{h},$

$\displaystyle
N_{h}(\Psi_{h};\Phi_{h}):=A_{h}(\Psi_{h},\Phi_{h})+B_{1}(\Psi_{h},\Phi_{h})+B_{2}(\Psi_{h},\Psi_{h},\Phi_{h})+B_{3}(\Psi_{h},\Psi_{h},\Psi_{h},\Phi_{h})-L_{h}(\Phi_{h})=0.$
(3.33)

For $\Theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})$ and
$\Phi=(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4})\in\mathbf{X},$
$A_{h}(\Theta,\Phi):=\sum_{i=1}^{4}a_{h}(\theta_{i},\varphi_{i}),\text{ and
}L_{h}(\Phi)=\sum_{i=1}^{4}L^{i}_{h}(\varphi_{i}).$ Note that the forms
$B_{i},i=1,2,3$ are defined in Section 3.1. For $\theta,\varphi\in
H^{1}(\Omega)$, $\mathbf{g}=(g_{1},g_{2},g_{3},g_{4})$,

$\displaystyle
a_{h}(\theta,\varphi):=\int_{\Omega}\nabla\theta\cdot\nabla\varphi\,{\rm
dx}-\int_{\partial\Omega}\frac{\partial\theta}{\partial\nu}\varphi\,{\rm
ds}-\int_{\partial\Omega}\theta\frac{\partial\varphi}{\partial\nu}\,{\rm
ds}+\sum_{E\in\mathcal{E}(\partial\Omega)}\frac{\sigma}{h_{E}}\int_{E}\theta\varphi\,{\rm
ds},$ $\displaystyle\text{ and
}\,\,L^{i}_{h}(\varphi):=-\int_{\partial\Omega}g_{i}\frac{\partial\varphi}{\partial\nu}\,{\rm
ds}+\sum_{E\in\mathcal{E}(\partial\Omega)}\frac{\sigma}{h_{E}}\int_{E}g_{i}\varphi\,{\rm
ds}\text{ for }1\leq i\leq 4,$

where $\nu$ denotes the outward unit normal associated to $\partial\Omega.$
For all $\Theta_{h},\Phi_{h}\in\mathbf{X}_{h}$, define the discrete bilinear
form in this case as

$\displaystyle\langle{DN}_{h}(\Psi)\Theta_{h},\Phi_{h}\rangle:=A_{h}(\Theta_{h},\Phi_{h})+B_{1}(\Theta_{h},\Phi_{h})+2B_{2}(\Psi,\Theta_{h},\Phi_{h})+3B_{3}(\Psi,\Psi,\Theta_{h},\Phi_{h})$

and the perturbed bilinear form as

$\langle{DN}_{h}(\textrm{I}_{h}\Psi)\Theta_{h},\Phi_{h}\rangle:=A_{h}(\Theta_{h},\Phi_{h})+B_{1}(\Theta_{h},\Phi_{h})+2B_{2}(\textrm{I}_{h}\Psi,\Theta_{h},\Phi_{h})+3B_{3}(\textrm{I}_{h}\Psi,\textrm{I}_{h}\Psi,\Theta_{h},\Phi_{h}).$

######  Theorem 3.17 (Existence, uniqueness and error estimates).

Let $\Psi$ be a regular solution of the non-linear system (3.2) such that
(3.1) holds. For a given fixed $\ell>0,$ a sufficiently large $\sigma$ and a
sufficiently small discretization parameter chosen as
$h=O(\ell^{1+\varsigma})$ for any $\varsigma>0$, there exists a unique
solution $\Psi_{h}$ of the discrete non-linear problem (3.33) that
approximates $\Psi$ such that

$\displaystyle(i)\,{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{h}\lesssim
h,\,\text{ and
}\,\,(ii)\>{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{0}\lesssim
h^{2}(1+(1+\ell^{-1})^{2}).$

The proof of Theorem 3.17$(i)$ follows similar methodology of Theorem 3.3 with
the choice of the non-linear map
$\mu_{h}:\mathbf{X}_{h}\rightarrow\mathbf{X}_{h}$ defined by: for
$\Theta_{h},\Phi_{h}\in\mathbf{X}_{h}$,

$\displaystyle\langle{DN}_{h}({{\rm
I}_{h}\Psi})\mu_{h}(\Theta_{h}),\Phi_{h}\rangle$
$\displaystyle=2B_{2}(\textrm{I}_{h}\Psi,\Theta_{h},\Phi_{h})+3B_{3}(\textrm{I}_{h}\Psi,\textrm{I}_{h}\Psi,\Theta_{h},\Phi_{h})-B_{2}(\Theta_{h},\Theta_{h},\Phi_{h})$
$\displaystyle\quad-
B_{3}(\Theta_{h},\Theta_{h},\Theta_{h},\Phi_{h})+L_{h}(\Phi_{h}).$

###### Proof of Theorem 3.17$(ii)$.

Set $G:=\textrm{I}_{h}\Psi-\Psi_{h}\in\mathbf{L}^{2}(\Omega).$ Multiply (3.29)
by $\Phi_{h}=\textrm{I}_{h}\Psi-\Psi_{h}$, use integration by parts, and then
add and subtract an intermediate term as

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{h}\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{0}^{2}=\langle{DN}_{h}(\Psi)\textrm{I}_{h}\Psi-\Psi_{h},\boldsymbol{\chi}\rangle$
$\displaystyle=\langle{DN}_{h}(\Psi)\textrm{I}_{h}\Psi-\Psi,\boldsymbol{\chi}\rangle+\langle{DN}_{h}(\Psi)\Psi-\Psi_{h},\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi}\rangle+\langle{DN}_{h}(\Psi)\Psi-\Psi_{h},\textrm{I}_{h}\boldsymbol{\chi}\rangle=:T_{1}+T_{2}+T_{3}.$
(3.34)

Set
$B_{L_{1}}(\cdot,\cdot):=B_{1}(\cdot,\cdot)+2B_{2}(\Psi,\cdot,\cdot)+3B_{3}(\Psi,\Psi,\cdot,\cdot),$
then the terms

$\displaystyle
T_{1}:=A_{h}(\textrm{I}_{h}\Psi-\Psi,\boldsymbol{\chi})+B_{L_{1}}(\textrm{I}_{h}\Psi-\Psi,\boldsymbol{\chi}),\,\,T_{2}:=A_{h}(\Psi-\Psi_{h},\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi})+B_{L_{1}}(\Psi-\Psi_{h},\boldsymbol{\chi}-\textrm{I}_{h}\boldsymbol{\chi}).$

A use of definition of $A_{h}(\cdot,\cdot)$, $\boldsymbol{\chi}=0$ on
$\partial\Omega$, integration by parts and cancellation of terms yields

$A_{h}({\rm I}_{h}\Psi-\Psi,\boldsymbol{\chi})=\int_{\Omega}\nabla({\rm
I}_{h}\Psi-\Psi)\cdot\nabla\boldsymbol{\chi}\,{\rm
dx}-\int_{\partial\Omega}\frac{\partial\boldsymbol{\chi}}{\partial\nu}\cdot({\rm
I}_{h}\Psi-\Psi)\,{\rm ds}=-\int_{\Omega}({\rm
I}_{h}\Psi-\Psi)\cdot\Delta\boldsymbol{\chi}\,{\rm dx}.$

This plus Hölder’s inequality and Lemma 3.8 show

$A_{h}({\rm
I}_{h}\Psi-\Psi,\boldsymbol{\chi})\leq{|\kern-1.07639pt|\kern-1.07639pt|{\rm
I}_{h}\Psi-\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{0}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\lesssim
h^{2}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}.$

The remaining $B_{i},i=1,2,3$ terms in $T_{1}$ are estimated using Lemma
3.7$(i)$-$(iii)$, and Lemma 3.8. The fact that $\boldsymbol{\chi}-{\rm
I}_{h}\boldsymbol{\chi}=0$ on $\partial\Omega,$ Lemma 3.8 and Theorem
3.17$(i)$ leads to the estimate of $A_{h}(\cdot,\cdot)$ term in $T_{2}$ as
$A_{h}(\Psi-\Psi_{h},\boldsymbol{\chi}-{\rm I}_{h}\boldsymbol{\chi})\lesssim
h{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}{|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{h}|\kern-1.07639pt|\kern-1.07639pt|}_{h}\lesssim
h^{2}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}.$
The term $T_{3}$ with $\eqref{discrete nonlinear problem ferro Nitsche's
method}$ and $N_{h}(\Psi,\textrm{I}_{h}\boldsymbol{\chi})=0$ reads

$\displaystyle
T_{3}:=\langle{DN}_{h}(\Psi)\Psi-\Psi_{h},\textrm{I}_{h}\boldsymbol{\chi}\rangle=\langle{DN}_{h}(\Psi)\Psi-\Psi_{h},\textrm{I}_{h}\boldsymbol{\chi}\rangle+{N}_{h}(\Psi_{h},\textrm{I}_{h}\boldsymbol{\chi})-{N}_{h}(\Psi,\textrm{I}_{h}\boldsymbol{\chi})$
$\displaystyle\quad=B_{2}({\Psi}_{h}-{\Psi},{\Psi}_{h}-{\Psi},\textrm{I}_{h}\boldsymbol{\chi})+(2B_{3}(\Psi,\Psi,\Psi,{\rm
I}_{h}\boldsymbol{\chi})-3B_{3}(\Psi,\Psi,\Psi_{h},{\rm
I}_{h}\boldsymbol{\chi})+B_{3}(\Psi_{h},\Psi_{h},\Psi_{h},{\rm
I}_{h}\boldsymbol{\chi})).$

The term $T_{3}$ and the rest of The terms in $T_{2}$ are directly comparable
to the terms of Theorem 3.5. The rest of the proof utilizes similar
methodology of Theorem 3.5 and hence is skipped here. ∎

######  Remark 3.18.

The boundary term in $T_{1}$ of Theorem 3.17$(ii)$ appearing due to the
integration by parts gets cancelled with the boundary term in $A_{h}({\rm
I}_{h}\Psi-\Psi,\boldsymbol{\chi})$, whereas similar type of boundary term in
(3.32) of Theorem 3.5 leads to the sub-optimal convergence rate
$\mathcal{O}(h^{\frac{3}{2}})$ in $\mathbf{L}^{2}$ norm. ∎

##  4 Numerical experiments

This section reports on numerical experiments for the benchmark problem [3]
for dilute ferronematic suspensions, on a re-scaled two-dimensional square
domain $\Omega=(0,1)\times(0,1)$ with a uniform refinement strategy. Numerical
solutions approximate the regular solutions of (3.2) for a fixed value of the
parameters $\ell$ and $c.$ The discrete solution landscapes of (3.2), for
various parameter ($\ell,c$) values, the associated computational errors and
convergence rates are explored for conforming FEM. Let $e_{i}$ and $h_{i}$ be
the error and the mesh parameter at $i$-th level, respectively. The $i$-th
level experimental order of convergence is defined by
$\displaystyle\alpha_{i}:=\log(e_{i}/e_{i+1})/\log(h_{i}/h_{i+1})$ for
$i=1,\ldots,n-1$ and $n$ is the final iteration considered in numerical
experiments. Newton’s method is applied to approximate the solutions of (3.2).
For detailed construction of the initial conditions/profiles, we refer to [20,
26, 33].

Solution | $x=0$ | $x=1$ | $y=0$ | $y=1$
---|---|---|---|---
$Q_{11}$ | -1 | -1 | 1 | 1
$Q_{12}$ | 0 | 0 | 0 | 0
$M_{1}$ | 0 | 0 | -1 | 1
$M_{2}$ | 1 | -1 | 0 | 0
Table 1:  Tangential boundary conditions for solution components
$Q_{11},Q_{12},M_{1},M_{2}.$

The Dirichlet tangent boundary conditions [3, 26, 34, 35] are detailed in
Table 1. The natural mismatch in the tangent boundary conditions of the
director $\mathbf{n}$ and magnetization vector leads to the corner defects. We
construct a Lipschitz continuous boundary condition $\mathbf{g}$ using the
tangential boundary condition in Table 1 and trapezoidal shape functions, [26]
$\textit{T}_{d}:[0,1]\rightarrow{\mathbb{R}}$ defined as

$\mathbf{g}=\begin{cases}(\textit{T}_{d}(x),0,-\textit{T}_{d}(x),0)&\text{on}\,\,\,\,y=0,\\\
(\textit{T}_{d}(x),0,\textit{T}_{d}(x),0)&\text{on}\,\,\,\,y=1,\\\
(-\textit{T}_{d}(y),0,0,\textit{T}_{d}(y))&\text{on}\,\,\,\,x=0,\\\
(-\textit{T}_{d}(y),0,0,-\textit{T}_{d}(y))&\text{on}\,\,\,\,x=1,\end{cases}\,\,\text{
and }\,\,\textit{T}_{d}(t)=\begin{cases}t/d,&0\leq t\leq d,\\\ 1,&d\leq t\leq
1-d,\\\ (1-t)/d,&1-d\leq t\leq 1,\end{cases}$

where the parameter $d=3\sqrt{\ell}$, is the size of mismatch region. For
small choices of the parameter $\ell,$ the numerical results are divided into
three categories according to the positive, negative and zero value of the
coupling parameter $c$.

(a) $\mathbf{Q}_{D1}$ and $\mathbf{M}$ profile

(b) $\mathbf{Q}_{R4}$ and $\mathbf{M}$ profile

(c)

(d)

Figure 1: Nematic $\mathbf{Q}$ and magnetic $\mathbf{M}$ configurations: (a)
D1 diagonal nematic and uncoupled magnetic profile and (b) R4 rotated nematic
and uncoupled magnetic profile for the parameter values $\ell=0.001,c=0;$
Energy and $\mathbf{L}^{2}$ norm error versus discretization parameter $h$
plots for the discrete solutions (c) $\Psi_{h}=(\mathbf{Q}_{D1},\mathbf{M})$,
(d) $\Psi_{h}=(\mathbf{Q}_{R4},\mathbf{M})$ for two sets of parameter values
$\ell=0.001,c=0$ and $\ell=0.002,c=0$.

Case I: the coupling parameter $\mathbf{c=0}$

The two-dimensional planar bistable nematic device [26] (uncoupled system
i.e., $c=0$) exhibits two sets of equilibrium configurations- 1) diagonally
stable: the nematic directors roughly align along one of the square diagonals
and there are two classes of diagonal solutions: D1 and D2, one for each
square diagonal; 2) rotated states: here, the nematic director rotates by
$\pi$ radians between a pair of opposite parallel edges, and there are $4$
classes of rotated solutions labelled by R1, R2, R3 and R4 respectively,
related to each other by $\frac{\pi}{2}$ radians. The diagonal and rotated
solutions are distinguished by the locations of the splay vertices; a splay
vertex being a vertex such that the nematic director splays around the vertex
and a bend vertex being such that the nematic director bends around the vertex
in question. Each diagonal solution has a pair of diagonally opposite splay
vertices and each rotated solution has a pair of adjacent splay vertices,
connected by a square edge. In Figure 1(a) and Figure 1(b), the discrete
solutions, $\Psi_{h}=(\mathbf{Q}_{D1},\mathbf{M})$ and
$\Psi_{h}=(\mathbf{Q}_{R4},\mathbf{M})$, are plotted for $\ell=0.001,c=0$.
Here $\mathbf{Q}_{D1}$ (resp. $\mathbf{Q}_{R4}$) is the D1 diagonal (resp. R4
rotated) solution with defects at vertices, and the corresponding nematic
director $\mathbf{n}=(\cos\theta,\sin\theta)$ where
$\theta=\frac{1}{2}\textrm{atan}\frac{Q_{D1,12}}{Q_{D1,11}}$ and
$Q_{D1,11}=\frac{|\mathbf{Q}_{D1}|}{\sqrt{2}}\cos
2\theta,Q_{D1,12}=\frac{|\mathbf{Q}_{D1}|}{\sqrt{2}}\sin 2\theta$ are the two
independent components of $\mathbf{Q}_{D1}$. Analogous remarks apply to
$\mathbf{Q}_{R4}$. $\mathbf{M}$ labels the uncoupled magnetic profile with a
$+1$-degree vortex at the square center consistent with topologically non-
trivial boundary conditions. The magnetization vector $\mathbf{M}$ has a
direction whereas the nematic director field, $\mathbf{n}$, is plotted without
a direction since $\mathbf{n}$ and $-\mathbf{n}$ are physically equivalent.
Figure 1(c) (resp. Figure 1(d)) demonstrate the convergence history of the
discrete solutions, computed using piecewise polynomials of degree $1$,
associated with D1 ( resp. R4) nematic solutions, in energy and
$\mathbf{L}^{2}$ norms for the parameter values $\ell=0.001,c=0$ and
$\ell=0.002,c=0$. The order $1$ convergence in energy norm and order $2$
convergence in $\mathbf{L}^{2}$\- norm are obtained for both sets of parameter
values. The color bars for nematic and magnetic profiles plot the values of
$s=\sqrt{Q_{11}^{2}+Q_{12}^{2}}$ and
$\left\lvert\mathbf{M}\right\rvert=\sqrt{M_{1}^{2}+M_{2}^{2}},$ respectively.
The lines and arrows depict $\mathbf{n}$ and $\mathbf{M}$ respectively. Note
that all subsequent discrete solution profiles, $\Psi_{h}$, have the nematic
director field plot on the left and and magnetization vector plot on the
right.

(a) $\mathbf{Q}_{D1}$ and $\mathbf{M}_{D1}$ profile

(b) $\mathbf{Q}_{D2}$ and $\mathbf{M}_{D2}$ profile

(c) $\mathbf{Q}_{R1}$ and $\mathbf{M}_{R1}$ profile

(d) $\mathbf{Q}_{R2}$ and $\mathbf{M}_{R2}$ profile

(e) $\mathbf{Q}_{R3}$ and $\mathbf{M}_{R3}$ profile

(f) $\mathbf{Q}_{R4}$ and $\mathbf{M}_{R4}$ profile

Figure 2: Nematic $\mathbf{Q}$ and magnetic $\mathbf{M}$ configurations for
$\ell=0.001$ and $c=0.25.$ Left column: solution profiles
$\Psi_{h}=(\mathbf{Q}_{D1},\mathbf{M}_{D1})$ (top) and
$\Psi_{h}=(\mathbf{Q}_{D2},\mathbf{M}_{D2})$ (bottom) corresponding to
diagonal D1 and D2 nematic stable solutions, respectively; Middle column:
solution profiles $\Psi_{h}=(\mathbf{Q}_{R1},\mathbf{M}_{R1})$ (top) and
$\Psi_{h}=(\mathbf{Q}_{R2},\mathbf{M}_{R2})$ (bottom) corresponding to rotated
R1 and R2 nematic stable solutions, respectively; Right column: solution
profiles $\Psi_{h}=(\mathbf{Q}_{R3},\mathbf{M}_{R3})$ (top) and
$\Psi_{h}=(\mathbf{Q}_{R4},\mathbf{M}_{R4})$ (bottom) corresponding to rotated
R3 and R4 nematic stable solutions, respectively.

Case II: the coupling parameter $\mathbf{c>0}$

Figure 2 plots the numerically computed stable solution profiles for parameter
values $\ell=0.001,c=0.25$, represented by
$\Psi_{h}=(\mathbf{Q}_{D1},\mathbf{M}_{D1}),$
$\Psi_{h}=(\mathbf{Q}_{D2},\mathbf{M}_{D2}),$
$\Psi_{h}=(\mathbf{Q}_{R1},\mathbf{M}_{R1}),$
$\Psi_{h}=(\mathbf{Q}_{R2},\mathbf{M}_{R2}),$
$\Psi_{h}=(\mathbf{Q}_{R3},\mathbf{M}_{R3}),$
$\Psi_{h}=(\mathbf{Q}_{R4},\mathbf{M}_{R4}),$ corresponding to D1, D2 diagonal
and the R1, R2, R3, R4 rotated stable nematic equillibria, respectively. For
small $\ell$, and $c>0,$ the coupling energy favors the co-alignment of
$\mathbf{n}$ and $\mathbf{M}$, i.e., $\mathbf{n}\cdot\mathbf{M}=\pm 1.$ In
this case, the nematic profiles (both diagonal and rotated) do not exhibit any
interior vortices whereas the magnetic profiles develop an interior line of
reduced $|\mathbf{M}|$, analogous to a domain wall, smeared out along the
square diagonals/ near one of the square edges. The magnetic profiles,
$\mathbf{M}_{D1},\mathbf{M}_{D2}$ (resp. $\mathbf{M}_{R1},\mathbf{M}_{R2}$,
$\mathbf{M}_{R3},\mathbf{M}_{R4}$) exhibit $\pi$-walls [20] along the square
diagonals $y=x$ for D1 and along $y=-x+1$ for D2 (resp. along the square edges
$x=1$ for R1, $x=0$ for R2, $y=0$ for R3, $y=1$ for R4) nematic solutions. The
domain walls are created to ensure the compatibility between the angle
constraint in (2.2), the condition necessary to be minimizer for $c>0$, and
the tangent boundary conditions. Recall that for
$\mathbf{n}=(\cos\theta,\sin\theta)$, $Q_{11}=s\cos 2\theta$ and $Q_{12}=s\sin
2\theta.$ For a stable stationary point $\Psi=(s\cos 2\theta,s\sin
2\theta,\left\lvert\mathbf{M}\right\rvert\cos\varphi,\left\lvert\mathbf{M}\right\rvert\sin\varphi)$
with magnetization angle $\varphi$, (2.2) implies that

$\theta\approx\varphi+k\pi,\,k\in\mathbb{Z}$

almost everywhere in the domain interior , for sufficiently small values of
$\ell$.

(a)

(b)

Figure 3: Domain wall formation of magnetic profiles (a) $\mathbf{M}_{D1}$ and
(b) $\mathbf{M}_{D2}$.

(a)

(b)

(c)

Figure 4: Energy and $\mathbf{L}^{2}$ norm error versus discretization
parameter $h$ plots for the discrete solutions (a)
$\Psi_{h}=(\mathbf{Q}_{D1},\mathbf{M}_{D1})$ and
$\Psi_{h}=(\mathbf{Q}_{D2},\mathbf{M}_{D2})$, (b)
$\Psi_{h}=(\mathbf{Q}_{R1},\mathbf{M}_{R1})$ and
$\Psi_{h}=(\mathbf{Q}_{R2},\mathbf{M}_{R2})$, (c)
$\Psi_{h}=(\mathbf{Q}_{R3},\mathbf{M}_{R3})$ and
$\Psi_{h}=(\mathbf{Q}_{R4},\mathbf{M}_{R4}),$ for the two sets of parameter
values $\ell=0.001,c=0.25,$ and $\ell=0.002,c=0.25,$

The tangent boundary condition for D1 ( resp. D2) nematic profile is coded in
the boundary conditions: $\theta=0$ along $y=0,$ $\theta=0$ along $y=1,$
$\theta=\frac{\pi}{2}$ along $x=0,$ $\theta=\frac{\pi}{2}$ along $x=1$ (resp.
$\theta=\pi$ along $y=0,$ $\theta=\pi$ along $y=1,$ $\theta=\frac{\pi}{2}$
along $x=0,$ $\theta=\frac{\pi}{2}$ along $x=1$). Figure (3) shows that
$\theta\approx\varphi$ for $y\geq x$ (resp. $y\leq-x+1$) and
$\theta\approx\varphi+\pi$ for $y<x$ (resp. $y>-x+1$) in $\mathbf{M}_{D1}$
(resp. $\mathbf{M}_{D2}$) profile. The domain walls for
$\mathbf{M}_{R1},\mathbf{M}_{R2}$, $\mathbf{M}_{R3},\mathbf{M}_{R4}$ can be
interpreted similarly. Figure 4 illustrates the numerical errors and orders of
convergence, computed using piecewise polynomials of degree $1$, for the
discrete solutions (a) $\Psi_{h}=(\mathbf{Q}_{D1},\mathbf{M}_{D1})$ and
$\Psi_{h}=(\mathbf{Q}_{D2},\mathbf{M}_{D2})$, (b)
$\Psi_{h}=(\mathbf{Q}_{R1},\mathbf{M}_{R1})$ and
$\Psi_{h}=(\mathbf{Q}_{R2},\mathbf{M}_{R2})$, (c)
$\Psi_{h}=(\mathbf{Q}_{R3},\mathbf{M}_{R3})$ and
$\Psi_{h}=(\mathbf{Q}_{R4},\mathbf{M}_{R4})$, respectively, in energy and
$\mathbf{L}^{2}$ norms for the parameter values $\ell=0.001,c=0.25$ and
$\ell=0.002,c=0.25$. The convergence rates obtained in energy and
$\mathbf{L}^{2}$ norms are of $\mathcal{O}(h)$ and $\mathcal{O}(h^{2})$,
respectively.

(a) $\mathbf{Q}_{D1}^{1}$ and $\mathbf{M}^{1}_{D1}$ profile

(b) $\mathbf{Q}_{D1}^{2}$ and $\mathbf{M}^{2}_{D1}$ profile

(c) $\mathbf{Q}^{1}_{D2}$ and $\mathbf{M}^{1}_{D2}$ profile

(d) $\mathbf{Q}^{2}_{D2}$ and $\mathbf{M}^{2}_{D2}$ profile

(e) $\mathbf{Q}^{1}_{R1}$ and $\mathbf{M}^{1}_{R1}$ profile

(f) $\mathbf{Q}^{2}_{R1}$ and $\mathbf{M}^{2}_{R1}$ profile

Figure 5: Nematic $\mathbf{Q}$ and magnetic $\mathbf{M}$ configurations for
$\ell=0.001$ and $c=-0.25.$ Left column: two solution profiles
$\Psi_{h}=(\mathbf{Q}^{1}_{D1},\mathbf{M}^{1}_{D1})$ and
$\Psi_{h}=(\mathbf{Q}^{2}_{D1},\mathbf{M}^{2}_{D1})$ corresponding to diagonal
D1 nematic stable solution; Middle column: two solution profiles
$\Psi_{h}=(\mathbf{Q}^{1}_{D2},\mathbf{M}^{1}_{D2})$ and
$\Psi_{h}=(\mathbf{Q}^{2}_{D2},\mathbf{M}^{2}_{D2})$ corresponding to diagonal
D2 nematic stable solution; Right column: two solution profiles
$\Psi_{h}=(\mathbf{Q}^{1}_{R1},\mathbf{M}^{1}_{R1})$ and
$\Psi_{h}=(\mathbf{Q}^{2}_{R1},\mathbf{M}^{2}_{R2})$ corresponding to rotated
R1 nematic stable solution.

(a) $\mathbf{Q}^{1}_{R2}$ and $\mathbf{M}^{1}_{R2}$ profile

(b) $\mathbf{Q}^{2}_{R2}$ and $\mathbf{M}^{2}_{R2}$ profile

(c) $\mathbf{Q}^{1}_{R3}$ and $\mathbf{M}^{1}_{R3}$ profile

(d) $\mathbf{Q}^{2}_{R3}$ and $\mathbf{M}^{2}_{R3}$ profile

(e) $\mathbf{Q}^{1}_{R4}$ and $\mathbf{M}^{1}_{R4}$ profile

(f) $\mathbf{Q}^{2}_{R4}$ and $\mathbf{M}^{2}_{R4}$ profile

Figure 6: Nematic $\mathbf{Q}$ and magnetic $\mathbf{M}$ configurations for
$\ell=0.001$ and $c=-0.25.$ Left column: two solution profiles
$\Psi_{h}=(\mathbf{Q}^{1}_{R2},\mathbf{M}^{1}_{R2})$ and
$\Psi_{h}=(\mathbf{Q}^{2}_{R2},\mathbf{M}^{2}_{R2})$ corresponding to rotated
R2 nematic stable solution; Middle column: two solution profiles
$\Psi_{h}=(Q^{1}_{R3},\mathbf{M}^{1}_{R3})$ and
$\Psi_{h}=(\mathbf{Q}^{2}_{R3},\mathbf{M}^{2}_{R3})$ corresponding to rotated
R3 nematic stable solution; Right column: two solution profiles
$\Psi_{h}=(\mathbf{Q}^{1}_{R4},\mathbf{M}^{1}_{R4})$ and
$\Psi_{h}=(\mathbf{Q}^{2}_{R4},\mathbf{M}^{2}_{R4})$ corresponding to rotated
R4 nematic stable solution.

Case III: the coupling parameter $c<0$

Now, we discuss the discrete solution profiles for negative coupling, which
favours perpendicular alignment of $\mathbf{n}$ and $\mathbf{M}$, i.e,
$\mathbf{n}\cdot\mathbf{M}=0,$ for the parameter values $\ell=0.001$ and
$c=-0.25.$ For diagonal nematic solutions (resp. rotated solutions), the
symmetry between the diagonally opposite (resp. square edge) splay vertices is
broken, so that there are $4$ distinct diagonal solutions. By similar
reasoning, there are $8$ distinct rotated solutions, so that the number of
stable admissible equilibria is doubled.

(a)

(b)

Figure 7: Energy and $\mathbf{L}^{2}$ norm errors versus discretization
parameter $h$ plots for the discrete solutions (a)
$\Psi_{h}=(\mathbf{Q}_{D1}^{1},\mathbf{M}_{D1}^{1})$ and
$\Psi_{h}=(\mathbf{Q}_{D1}^{2},\mathbf{M}_{D1}^{2})$, (b)
$\Psi_{h}=(\mathbf{Q}_{D2}^{1},\mathbf{M}_{D2}^{1})$ and
$\Psi_{h}=(\mathbf{Q}_{D2}^{2},\mathbf{M}_{D2}^{2})$, for two sets of
parameter values $\ell=0.001,c=-0.25,$ and $\ell=0.002,c=-0.25.$

(a)

(b)

Figure 8: Energy and $\mathbf{L}^{2}$ norm errors versus discretization
parameter $h$ plots for the discrete solutions (a)
$\Psi_{h}=(\mathbf{Q}_{R1}^{1},\mathbf{M}_{R1}^{1})$ and
$\Psi_{h}=(\mathbf{Q}_{R1}^{2},\mathbf{M}_{R1}^{2})$, (b)
$\Psi_{h}=(\mathbf{Q}_{R2}^{1},\mathbf{M}_{R2}^{1})$ and
$\Psi_{h}=(\mathbf{Q}_{R2}^{2},\mathbf{M}_{R2}^{2})$, for two sets of
parameter values $\ell=0.001,c=-0.25,$ and $\ell=0.002,c=-0.25.$

(a)

(b)

Figure 9: Energy and $\mathbf{L}^{2}$ norm errors versus discretization
parameter $h$ plots for the discrete solutions, (a)
$\Psi_{h}=(\mathbf{Q}_{R3}^{1},\mathbf{M}_{R3}^{1})$ and
$\Psi_{h}=(\mathbf{Q}_{R3}^{2},\mathbf{M}_{R3}^{2})$, (b)
$\Psi_{h}=(\mathbf{Q}_{R4}^{1},\mathbf{M}_{R4}^{1})$ and
$\Psi_{h}=(\mathbf{Q}_{R4}^{2},\mathbf{M}_{R4}^{2}),$ for two sets of
parameter values $\ell=0.001,c=-0.25,$ and $\ell=0.002,c=-0.25.$

For instance, $\Psi_{h}=(\mathbf{Q}_{D1}^{1},\mathbf{M}_{D1}^{1})$ and
$\Psi_{h}=(\mathbf{Q}_{D1}^{2},\mathbf{M}_{D1}^{2})$ in Figure 5 are two
distinct stable, numerically computed solutions corresponding to standard D1
diagonal nematic profile (for $c=0$). Similarly, there are ten pairs of
distinct stable solution profiles, corresponding to the standard D2, R1, R2,
R3 and R4 profiles; see Figures 5 and 6. The numerical errors and orders of
convergence of the discrete solutions associated with the diagonal (D1, D2),
and rotated (R1, R2) and (R3, R4) nematic equilibria are plotted in Figures 7,
8 and 9, respectively, for two sets of parameter values $\ell=0.001,c=-0.25$
and $\ell=0.002,c=-0.25$. The convergence rates in energy and $\mathbf{L}^{2}$
norms are noted to be of order, $O(h)$ and $O(h^{2})$, respectively.

Parameter dependent plots

Figure 10 (resp. Figure 11) presents discretization parameter $h$ versus
energy and $\mathbf{L}^{2}$ norm error plots, for various values of $\ell$ and
positive coupling (resp. negative coupling) parameter, for the discrete
solution corresponding to D1 diagonal nematic equillibria. We observe that
both the energy and $\mathbf{L}^{2}$ norm errors are sensitive to the choice
of the small parameter $\ell$, for both instances of positive and negative
nemato-magnetic coupling. For instance, fix $h=10^{-2},$ the energy norm error
for $\ell=0.01$ at $h=10^{-2}$ is smaller than the error for $\ell=0.006$ and
similarly, the error increases as $\ell$ further decreases.

(a)

(b)

Figure 10: Convergence behavior plots of error in the (a) energy norm and (b)
$\mathbf{L}^{2}$ norm versus the discretization parameter $h$ for
$\Psi_{h}:=(\mathbf{Q}_{D1},\mathbf{M}_{D1})$ solution for various values of
$\ell$ and $c=0.25.$

(a)

(b)

Figure 11: Convergence behavior plots of error in the (a) energy norm and (b)
$\mathbf{L}^{2}$ norm versus the discretization parameter $h$ for
$\Psi_{h}:=(\mathbf{Q}_{D1}^{1},\mathbf{M}_{D1}^{1})$ solution for various
values of $\ell$ and $c=-0.25.$

######  Remark 4.1.

The Landau-de Gennes energy for nematic liquid crystal is defined [26] as

$\displaystyle
F_{\text{nem}}(\bar{{\Psi}}_{\epsilon})=\int_{\Omega}(\left\lvert\nabla\bar{{\Psi}}_{\epsilon}\right\rvert^{2}+\epsilon^{-2}(\left\lvert\bar{{\Psi}}_{\epsilon}\right\rvert^{2}-1)^{2})\,{\rm
dx},$

where $\bar{{\Psi}}_{\epsilon}:=(u_{1},u_{2})=\bar{g}$ on $\partial\Omega$ and
$\epsilon$ is a material-dependent parameter that depends on the elastic
constant, domain size and temperature. The Euler-Lagrange equations are a
system of second order non-linear elliptic partial differential equations that
seeks $\bar{{\Psi}}_{\epsilon}\in X\times X$ such that for all
$\bar{\varphi}:=(\varphi_{1},\varphi_{2})\in V\times V$

$\displaystyle\int_{\Omega}\nabla u_{1}\cdot\nabla\varphi_{1}\,{\rm
dx}+2\epsilon^{-2}\int_{\Omega}(u_{1}^{2}+u_{2}^{2}-1)u_{1}\varphi_{1}\,{\rm
dx}=0,$ $\displaystyle\int_{\Omega}\nabla u_{2}\cdot\nabla\varphi_{2}\,{\rm
dx}+2\epsilon^{-2}\int_{\Omega}(u_{1}^{2}+u_{2}^{2}-1)u_{2}\varphi_{2}\,{\rm
dx}=0,$

which is (3.2) for ${\Psi}:=(u_{1},u_{2},0,0)$ and
${\Phi}:=(\varphi_{1},\varphi_{2},0,0)$ with the parameter values
$\ell=\frac{\epsilon^{2}}{2}.$ Note that the non-linearity in reduced Landau-
de Gennes minimization problem is cubic. The quadratic non-linear term
$B_{2}(\cdot,\cdot,\cdot)$ in (3.2) is zero here. An a priori error analysis
with $h-\epsilon$ dependency has been discussed for this model for
discontinuous Galerkin method in [33]. The analysis for conforming finite
element method is a special case of the problem considered in this paper. ∎

##  5 Conclusions

We study the minimizers of a $\mathbf{Q}$-tensor -$\mathbf{M}$ model for
dilute ferronematic suspensions in 2-D framework. The energy functional has
two parameters- a scaled elastic parameter $\ell$ and nemato-magnetic coupling
parameter $c$. We analyze the asymptotic behavior of the global minimizers of
$\tilde{\mathcal{E}}$, as $\ell\rightarrow 0$ and establish that
${|\kern-1.07639pt|\kern-1.07639pt|\Psi^{\ell}|\kern-1.07639pt|\kern-1.07639pt|}_{2}$
is bounded, independent of $\ell.$ This result plays a key role in the
$h-\ell$ dependent finite element analysis for regular solutions of the
corresponding Euler-Lagrange PDEs such that (3.1) holds. Whether the analysis
holds for all regular solutions of the Euler-Lagrange PDEs is proposed as
future work. The numerical results focus on the solution landscapes for
$\ell=0.001,c=0,\pm 0.25$ in a square domain, the convergence rates in energy
and $\mathbf{L}^{2}$ norms, the convergence behavior of discrete solutions for
various values of $\ell$. The numerical results in this manuscript can be
extended to stable solutions, for other values of $c$ and $\ell$ as reported
in [3]. The convergence of minimizers in $\mathbf{L}^{\infty}$ norm i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\Psi^{\ell}-\Psi_{0}|\kern-1.07639pt|\kern-1.07639pt|}_{\infty}$
estimates, the analysis for three-dimensional geometries, the finite element
analysis for polygonal domains with re-entrant corners and Dirichlet boundary
data with lesser regularity, a posteriori error analysis to investigate the
effects of defects on numerical errors are interesting and challenging
extensions of this work. Moreover, the asymptotic analysis of minimizers with
topologically non-trivial boundary conditions and/or $\ell$-dependent
Dirichlet boundary data $\mathbf{g}_{\ell}$, including star-shaped domains,
are further interesting areas to be investigated.

## Acknowledgements

R. R. Maity gratefully acknowledges Professor Yiwei Wang for his illuminating
suggestions in numerical computations, as well as Professor Giacomo Canevari
for helpful discussions. R. R. Maity also acknowledges the support from
institute Ph.D. fellowship. N. N. gratefully acknowledges SERB POWER
Fellowship SPF/2020/000019. A.M. acknowledges support from the University of
Strathclyde Global Engagement Fund and an OCIAM Visiting Fellowship, Visiting
Professorship from the University of Bath. A.M. acknowledges support from a
Leverhulme International Academic Fellowship. A.M. also acknowledges support
from the DST-UKIERI for the project on "Theoretical and experimental studies
of suspensions of magnetic nanoparticles, their applications and
generalizations".

## References

  * [1] F. Bethuel, H. Brezis, and F. Hélein, _Asymptotics for the minimization of a Ginzburg-Landau functional_ , Calculus of Variations and Partial Differential Equations 1 (1993), no. 2, 123–148.
  * [2] , _Ginzburg-Landau vortices_ , Modern Birkhäuser Classics, Birkhäuser/Springer, Cham, 2017, Reprint of the 1994 edition.
  * [3] K. Bisht, Y. Wang, V. Banerjee, and A. Majumdar, _Tailored morphologies in two-dimensional ferronematic wells_ , Physical Review E 101 (2020), 022706\.
  * [4] J. P. Borthagaray, R. H. Nochetto, and S. W. Walker, _A structure-preserving FEM for the uniaxially constrained Q-tensor model of nematic liquid crystals_ , Numerische Mathematik 145 (2020), no. 4, 837–881.
  * [5] F. Brochard and P. G. de Gennes, _Theory of magnetic suspensions in liquid crystals_ , Journal De Physique 31 (1970), no. 7, 691–708.
  * [6] S. V. Burylov and Y. L. Raikher, _Macroscopic properties of ferronematics caused by orientational interactions on the particle surfaces. I. extended continuum model_ , Molecular Crystals and Liquid Crystals Science and Technology. Section A. 258 (1995), no. 1, 107–122.
  * [7] M. C. Calderer, A. DeSimone, D. Golovaty, and A. Panchenko, _An effective model for nematic liquid crystal composites with ferromagnetic inclusions_ , SIAM Journal on Applied Mathematics 74 (2014), no. 2, 237–262.
  * [8] P. G. Ciarlet, _The finite element method for elliptic problems_ , Classics in Applied Mathematics, vol. 40, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.
  * [9] C. Cîrtoaje, E. Petrescu, C. Stan, and D. Creangă, _Ferromagnetic nanoparticles suspensions in twisted nematic_ , Physica E: Low-dimensional Systems and Nanostructures 79 (2016), 38 – 43.
  * [10] J. Dalby, Farrell P. E., Majumdar A., and J. Xia, _One-dimensional ferronematics in a channel: order reconstruction, bifurcations and multistability_ , https://arxiv.org/abs/2102.06347 (2021).
  * [11] T. A. Davis and E. C. Gartland, Jr., _Finite element analysis of the Landau-de Gennes minimization problem for liquid crystals_ , SIAM Journal on Numerical Analysis 35 (1998), no. 1, 336–362.
  * [12] P. de Gennes and J. Prost, _The physics of liquid crystals_ , International Series of Monogr, Clarendon Press, 1993.
  * [13] D. A. Di Pietro and A. Ern, _Mathematical aspects of discontinuous Galerkin methods_ , Mathématiques & Applications (Berlin) [Mathematics & Applications], vol. 69, Springer, Heidelberg, 2012.
  * [14] A. Ern and J-L. Guermond, _Theory and practice of finite elements_ , Applied Mathematical Sciences, vol. 159, Springer-Verlag, New York, 2004.
  * [15] L. C. Evans, _Partial differential equations_ , second ed., Graduate Studies in Mathematics, vol. 19, American Mathematical Society, Providence, RI, 2010.
  * [16] C. Giacomo, _Biaxiality in the asymptotic analysis of a 2D Landau–de Gennes model for liquid crystals_ , ESAIM. Control, Optimisation and Calculus of Variations 21 (2015), no. 1, 101–137.
  * [17] D. Golovaty, J. Alberto Montero, and P. Sternberg, _Dimension reduction for the Landau-de Gennes model in planar nematic thin films_ , Journal of Nonlinear Science 25 (2015), 1431 – 1451.
  * [18] P. Grisvard, _Elliptic problems in nonsmooth domains_ , Classics in Applied Mathematics, vol. 69, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011.
  * [19] M. D. Gunzburger and S. L. Hou, _Treating inhomogeneous essential boundary conditions in finite element methods and the calculation of boundary stresses_ , SIAM Journal on Numerical Analysis 29 (1992), no. 2, 390–424.
  * [20] Y. Han, J. Harris, Majumdar A., and J. Walton, _Tailored morphologies in two-dimensional ferronematic wells_ , Physical Review E (2021).
  * [21] Y. Han, A. Majumdar, and L. Zhang, _A reduced study for nematic equilibria on two-dimensional polygons_ , SIAM Journal on Applied Mathematics 80 (2020), no. 4, 1678–1703.
  * [22] H. B. Keller, _Approximation methods for nonlinear problems with application to two-point boundary value problems_ , Mathematics of Computation 29 (1975), 464–474.
  * [23] S. Kesavan, _Topics in functional analysis and applications_ , John Wiley & Sons, Inc., New York, 1989.
  * [24] J. P.F. 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), no. 6, 1387–1412.
  * [25] J.-L. Lions and E. Magenes, _Non-homogeneous boundary value problems and applications. Vol. I_ , Springer-Verlag, New York-Heidelberg, 1972, Translated from the French by P. Kenneth, Die Grundlehren der mathematischen Wissenschaften, Band 181.
  * [26] C. Luo, A. Majumdar, and R. Erban, _Multistability in planar liquid crystal wells_ , Physical Review E 85 (2012), 061702.
  * [27] R. R. Maity, A. Majumdar, and N. Nataraj, _Error analysis of Nitsche’s and discontinuous Galerkin methods of a reduced Landau–de Gennes problem_ , Computational Methods in Applied Mathematics 21 (2021), no. 1, 179 – 209.
  * [28] 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), no. 1, 227–280.
  * [29] A. Mertelj and D. Lisjak, _Ferromagnetic nematic liquid crystals_ , Liquid Crystals Reviews 5 (2017), no. 1, 1–33.
  * [30] A. Mertelj, D. Lisjak, M. Drofenik, and M. Copič, _Ferromagnetism in suspensions of magnetic platelets in liquid crystal_ , Nature 504 (2013), no. 7479, 237—241.
  * [31] R. Moser, _Partial regularity for harmonic maps and related problems_ , World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2005.
  * [32] J. Nitsche, _Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind_ , Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg 36 (1971), no. 1, 9–15.
  * [33] R. R. Maity, A. Majumdar, and N. Nataraj, _Discontinuous Galerkin finite element methods for the Landau–de Gennes minimization problem of liquid crystals_ , IMA Journal of Numerical Analysis 41 (2021), no. 2, 1130–1163.
  * [34] M. Slavinec, E. Klemenčič, M. Ambrožič, and M. Krasna, _Impact of nanoparticles on nematic ordering in square wells_ , Advances in Condensed Matter Physics 2015 (2015), 1–11.
  * [35] C. Tsakonas, A. J. Davidson, C. V. Brown, and N. J. Mottram, _Multistable alignment states in nematic liquid crystal filled wells_ , Applied Physics Letters 90 (2007), Article 111913.
  * [36] W. Wang, L. Zhang, and P. Zhang, _Modeling and computation of liquid crystals_ , Acta Numerica (2022).
  * [37] Y. Wang, G. Canevari, and A. Majumdar, _Order reconstruction for nematics on squares with isotropic inclusions: a Landau–de Gennes study_ , SIAM Journal on Applied Mathematics 79 (2019), no. 4, 1314–1340.