$\Delta\varpi_{\text{D}}^{(1)}=-\boldsymbol{\mathbf{q}}_{\text{D}}\cdot\nabla\vartheta^{(0)}-{Q}_{\text{D}}\vartheta^{(0)}\qquad\text{and}\qquad\Delta\varpi_{\text{N}}^{(1)}=-\boldsymbol{\mathbf{q}}_{\text{N}}\cdot\nabla\vartheta^{(0)}-{Q}_{\text{N}}\vartheta^{(0)}.$

Hence, based on (63) and (E.11), we can get an estimate for
$(\varpi_{\text{D}}^{(1)},\varpi_{\text{N}}^{(1)},\vartheta^{(1)})$ given by

$\left\|{\varpi_{\text{D},\text{N}}^{(1)}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left\|{\vartheta^{(1)}}\right\|^{(2+\alpha)}_{{\Gamma;\,[0,t]}}\leqslant\underaccent{\tilde}{c}\left|{\vartheta^{(0)}}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}\leqslant\underaccent{\tilde}{c}_{1}K(t),$

for some constants $\underaccent{\tilde}{c}_{1}>0$.

Let us write $\sigma_{m}(t)=\sum_{j=1}^{m}a_{j}(t)$ where $a_{m}(t)$ is given
by

$a_{m}(t):=\left\|{\varpi_{\text{D},\text{N}}^{(m)}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left\|{\vartheta^{(m)}}\right\|^{(2+\alpha)}_{{\Gamma;\,[0,t]}}$

As shown previously, we have
$\sigma_{1}(t)=a_{1}(t)\leqslant\underaccent{\tilde}{c}_{1}K(t)$. So, from
(62) and (63) with $\varepsilon_{6}>0$ taken sufficiently small, we can get a
bound for the sum $\sigma_{m}(t)=\sum_{j=1}^{m}a_{j}(t)$ given as follows

$\sigma_{m}(t)\leqslant
c\left(K(t)+\sum_{j=1}^{m}\int_{0}^{t}{{\max_{\Gamma}}}\left|{\frac{\partial}{\partial
s}}\vartheta^{(j)}(\cdot,s)\right|{{{d}}s}\right),$

for some constant $c>0$ (independent of $m$), for all $m=2,3,\ldots$. Applying
Grönwall’s lemma, we deduce that the sums $\sigma_{m}(t)$, for $m=2,3,\ldots$,
are also uniformly bounded above by $K(t)$. Thus, the series sums
$\sum_{j=1}^{m}a_{j}(t)$ actually converges. As a result, the sequence
$\\{(\varpi_{\text{D}}^{(m)},\varpi_{\text{N}}^{(m)},\vartheta^{(m)})\\}_{m}$
converges in the corresponding norm. Passing to the limit as $m\to\infty$, we
finally get the solution to system (50).

The Schauder method used above only proves the existence of a solution to
system (50), which at the same time satisfies the estimate (51). Thus, it only
remains to address the issue of the uniqueness of the solution. To do this, we
assume that two solutions,
$(\varpi_{\text{D}}^{1},\varpi_{\text{N}}^{1},\vartheta^{1})$ and
$(\varpi_{\text{D}}^{2},\varpi_{\text{N}}^{2},\vartheta^{2})$, to system (50)
exist. Clearly, their difference is also a solution to (50) with
$F_{\text{D}}\equiv 0$, $F_{\text{N}}\equiv 0$, and $\psi\equiv 0$. Moreover,
the estimate (51) remains valid, from which we see that
$\varpi_{\text{D}}^{1}-\varpi_{\text{D}}^{2}\equiv 0$,
$\varpi_{\text{N}}^{1}-\varpi_{\text{N}}^{2}\equiv 0$, and
$\vartheta^{1}-\vartheta^{2}\equiv 0$. This completes the proof of Lemma 7.2.
∎

##  8 Proof of the Main Result

We are now in the position to prove Theorem 4.4 by adapting a technique used
in [14, Section 5]. Also, we will verify at the end of the section our claims
in Theorem 3.1 by using Theorem 4.4 and through the change of variables.
Regarding the latter objective, the interpolation inequalities stated in the
lemma below will be central to our proof.

######  Lemma 8.1.

For some $\alpha\in(0,1)$, let $\Omega\subset\mathbb{R}^{d}$ be an open,
connected, bounded set of class ${{C}}^{2+\alpha}$ and
$u\in{{C}}^{2+\alpha}(\overline{\Omega})$. Then, there exist a constant
$\varepsilon_{5}>0$ such that for any $\varepsilon\in(0,\varepsilon_{5})$ the
following inequalities hold

$\displaystyle{{\max_{\overline{\Omega}}}}\left|{\nabla u}\right|$
$\displaystyle\leqslant\varepsilon^{1+\alpha}\left|{u}\right|_{\overline{\Omega}}^{(2+\alpha)}+\frac{c_{7}}{\varepsilon}{{\max_{\overline{\Omega}}}}\left|{u}\right|,$
(64) $\displaystyle\left|{u}\right|_{\overline{\Omega}}^{(2)}$
$\displaystyle\leqslant\varepsilon^{\alpha}\left|{u}\right|_{\overline{\Omega}}^{(2+\alpha)}+\frac{c_{8}}{\varepsilon^{2}}{{\max_{\overline{\Omega}}}}\left|{u}\right|,$
(65)

where $c_{7}:=c_{7}(d,\Omega,\varepsilon)$ and
$c_{8}:=c_{8}(d,\Omega,\varepsilon)$ are positive constants that depend only
on the dimension $d$, the number $\varepsilon$, and $\Omega$.

See Appendix E.4 for the proof.

###### Proof of Theorem 4.4.

Let us now recall the nonlinear problem (49). We will prove its solvability by
giving an estimate for the nonlinear terms

$\displaystyle\pazocal{{N}}[{V}_{i},{\tilde{\rho}}]$
$\displaystyle=-(\pazocal{{L}}_{{\tilde{\rho}}}-\pazocal{{L}}_{0}){V}_{i}-(\pazocal{{L}}_{{\tilde{\rho}}}-\pazocal{{L}}_{0}-\delta\pazocal{{L}}_{0})U_{i}$
$\displaystyle\qquad-(\pazocal{{M}}_{{\tilde{\rho}}}-\pazocal{{M}}_{0}){V}_{i}-(\pazocal{{M}}_{{\tilde{\rho}}}-\pazocal{{M}}_{0}-\delta\pazocal{{M}}_{0})U_{i},$
$\displaystyle\qquad-(\pazocal{{K}}_{\tilde{\rho}}{{V}_{i}}+\pazocal{{K}}_{{\tilde{\rho}}}{U_{i}})(\pazocal{{L}}_{{\tilde{\rho}}}{{\tilde{\rho}}}-\pazocal{{L}}_{0}{{\tilde{\rho}}})$
$\displaystyle\text{in $\Omega$},$ $\displaystyle\qquad i=\text{D},\text{N},$
$\displaystyle\pazocal{{B}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}]$
$\displaystyle=-\left[(B_{{\tilde{\rho}}}-B_{0}-\delta
B_{0})\dfrac{\partial{{{U}_{\text{N}\text{D}}^{0}}}}{\partial{\nu_{\circ}}}+(B_{{\tilde{\rho}}}-B_{0})\dfrac{\partial{{{{V}_{\text{N}\text{D}}}}}}{\partial{\nu_{\circ}}}\right],$
$\displaystyle\text{on $\Gamma$}.$

We start by noting the following formulae:

$\displaystyle(\pazocal{{F}}_{{\tilde{\rho}}}-\pazocal{{F}}_{0}){{V}}$
$\displaystyle=\int_{0}^{1}\frac{d}{d\lambda}(\pazocal{{F}}_{\lambda{\tilde{\rho}}}{{V}}){{{d}}\lambda},$
$\displaystyle(\pazocal{{F}}_{{\tilde{\rho}}}-\pazocal{{F}}_{0}-\delta\pazocal{{F}}_{0}){{V}}$
$\displaystyle=\int_{0}^{1}\left\\{\frac{d}{d\lambda}(\pazocal{{F}}_{\lambda{\tilde{\rho}}}{{V}})-\left.\frac{d}{d\mu}(\pazocal{{F}}_{\mu{\tilde{\rho}}}{{V}})\right|_{\mu=0}\right\\}{{{d}}\lambda}$
$\displaystyle=\int_{0}^{1}(1-\mu)\frac{{{d}}^{2}}{{{d}}\mu^{2}}(\pazocal{{F}}_{\mu{\tilde{\rho}}}{{V}}){{{d}}\lambda},$

where the operator
$\pazocal{{F}}\in\\{\pazocal{{L}},\pazocal{{M}},\pazocal{{B}}\\}$. Using these
formulae, and by expanding the derivative
$\frac{d}{d\lambda}({{F}}(\lambda{\tilde{\rho}},\lambda\nabla{\tilde{\rho}}))$,
it can be verified that
$(\pazocal{{F}}_{{\tilde{\rho}}}-\pazocal{{F}}_{0}){{V}}$,
$\pazocal{{F}}\equiv\pazocal{{L}}$,
${{V}}\in\\{{V}_{\text{D}},{V}_{\text{N}}\\}$, can be expressed as a linear
combination of terms containing the products
${\tilde{\rho}}{{V}}_{y_{j}y_{i}}$ and
${\tilde{\rho}}_{y_{k}}{{V}}_{y_{j}y_{i}}$, for $i,j,k=1,2\ldots,d$. Here, for
notational convenience, ${\tilde{\rho}}_{y_{k}}$ is used to denote the partial
derivative of ${\tilde{\rho}}$ with respect to the variable $y_{k}$ while
${{V}}_{y_{j}y_{i}}$ stands for the second-order partial derivative of $W$
with respect to the variables $y_{i}$ and $y_{j}$. Meanwhile,
$(\pazocal{{F}}_{{\tilde{\rho}}}-\pazocal{{F}}_{0}-\delta\pazocal{{F}}_{0}){U}_{0}$,
$\pazocal{{F}}\in\\{\pazocal{{L}},\pazocal{{M}},\pazocal{{B}}\\}$,
${U}_{0}\in\\{U^{0}_{\text{D}},U^{0}_{\text{N}}\\}$, can be written as a
linear combination of terms containing the products ${\tilde{\rho}}^{2}$,
${\tilde{\rho}}{\tilde{\rho}}_{y_{j}}$, and
${\tilde{\rho}}_{y_{k}}{\tilde{\rho}}_{y_{j}}$ (cf. [14, p. 131]). Note also
that the terms involving $\pazocal{{F}}_{{\tilde{\rho}}}-\pazocal{{F}}_{0}$,
$\pazocal{{F}}\in\\{\pazocal{{L}},\pazocal{{M}},\pazocal{{B}}\\}$, share
similar structures. Indeed, the following expansions hold:

$\displaystyle\frac{d}{d\lambda}({{F}}(\lambda{\tilde{\rho}},\lambda\nabla{\tilde{\rho}}))$
$\displaystyle={{F}}^{\prime}_{{\tilde{\rho}}}(\lambda{\tilde{\rho}},\nabla\lambda{\tilde{\rho}}){\tilde{\rho}}+\sum_{j=1}^{d}{{F}}^{\prime}_{{\tilde{\rho}}_{y_{j}}}(\lambda{\tilde{\rho}},\nabla\lambda{\tilde{\rho}}){\tilde{\rho}}_{y_{j}},$
$\displaystyle\frac{{{d}}^{2}}{{{d}}\mu^{2}}({{F}}(\mu{\tilde{\rho}},\mu\nabla{\tilde{\rho}}))$
$\displaystyle={{F}}^{\prime\prime}_{{\tilde{\rho}}{\tilde{\rho}}}(\mu{\tilde{\rho}},\nabla\mu{\tilde{\rho}}){\tilde{\rho}}^{2}+2\sum_{j=1}^{d}{{F}}^{\prime\prime}_{{\tilde{\rho}}{\tilde{\rho}}_{y_{j}}}(\mu{\tilde{\rho}},\nabla\mu{\tilde{\rho}}){\tilde{\rho}}{\tilde{\rho}}_{y_{j}}$
$\displaystyle\qquad\qquad+\sum_{k,j=1}^{d}{{F}}^{\prime\prime}_{{\tilde{\rho}}_{y_{k}}{\tilde{\rho}}_{y_{j}}}(\mu{\tilde{\rho}},\nabla\mu{\tilde{\rho}}){\tilde{\rho}}_{y_{k}}{\tilde{\rho}}_{y_{j}}.$

Hence, all the nonlinear terms of
$\pazocal{{N}}[{V}_{\text{N}},{\tilde{\rho}}]$ and
$\pazocal{{N}}[{V}_{\text{D}},{\tilde{\rho}}]$ consist of a multiplier term
${\tilde{\rho}}$ or a first-order partial derivative ${\tilde{\rho}}_{y_{j}}$.
Moreover, they are linear with respect to the second-order derivatives of the
unknowns ${V}_{\text{N}}$ and ${V}_{\text{D}}$. Thus, we only need to estimate
the product of two functions, one of which contains ${\tilde{\rho}}$ or
${\tilde{\rho}}_{y_{j}}$.

Let us estimate, for instance, the product
${\tilde{\rho}}_{y_{k}}{{V}}_{y_{j}y_{i}}$ over the domain $\Omega$, for
$i,j,k=1,\ldots,d$, ${{V}}\in\\{{V}_{\text{D}},{V}_{\text{N}}\\}$. By using
the interpolation inequality (64) we can get the estimate

$\left|{{\tilde{\rho}}_{y_{k}}{{V}}_{y_{j}y_{i}}}\right|_{\overline{\Omega}}^{(\alpha)}\leqslant{{\max_{\overline{\Omega}}}}\left|{{\tilde{\rho}}_{y_{k}}}\right|\left|{{{V}}}\right|_{\overline{\Omega}}^{(2+\alpha)}\leqslant\left(\varepsilon^{1+\alpha}\left|{{\tilde{\rho}}}\right|_{\overline{\Omega}}^{(2+\alpha)}+\frac{c_{7}}{\varepsilon}{{\max_{\overline{\Omega}}}}\left|{{\tilde{\rho}}}\right|\right)\left|{{{V}}}\right|_{\overline{\Omega}}^{(2+\alpha)},$

for some constant $c_{7}:=c_{7}(d)>0$. Because ${\tilde{\rho}}$ vanishes on
the exterior boundary $\Sigma$, the maximum principle implies that
${{\max_{\overline{\Omega}}}}\left|{{\tilde{\rho}}}\right|\leqslant{{\max_{\Gamma}}}\left|{{\tilde{\rho}}}\right|$.
Then, from the previous estimate, and in view of (E.5), we get

$\displaystyle\left|{{\tilde{\rho}}_{y_{k}}{{V}}_{y_{j}y_{i}}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}$
$\displaystyle\leqslant\left(\varepsilon^{1+\alpha}\left|{{\tilde{\rho}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\frac{c_{7}}{\varepsilon}\int_{0}^{t}{{\max_{\Gamma}}}\left|{{\tilde{\rho}}_{s}(\cdot,s)}\right|{{{d}}s}\right)\left|{{{V}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}$
$\displaystyle\leqslant\left|{{{V}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\left(\varepsilon^{1+\alpha}\left|{{\tilde{\rho}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\frac{c_{7}}{\varepsilon}t\left|{{\tilde{\rho}}_{\tau}}\right|^{\infty}_{{\Gamma;\,[0,t]}}\right).$

Letting $\varepsilon=t^{\frac{1}{2+\alpha}}$, we get the estimate

$\left|{{\tilde{\rho}}_{y_{k}}{{V}}_{y_{j}y_{i}}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}\leqslant\tilde{c}_{0}t^{\eta}\left|{{{V}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\left(\left|{{\tilde{\rho}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{\tilde{\rho}}_{\tau}}\right|^{\infty}_{{\Gamma;\,[0,t]}}\right),$

for some constant $\tilde{c}_{0}>0$, where
$\eta=\dfrac{1+\alpha}{2+\alpha}\in(0,1)$, for each function
${{V}}\in\\{{V}_{\text{D}},{V}_{\text{N}}\\}$.

We next estimate the nonlinear terms in the boundary condition. Because these
terms do not contain the second-order derivatives of the unknown functions –
in fact, on the boundary $\Gamma$, we only have linear combinations of the
products of ${\tilde{\rho}}$ and the unknown functions ${V}_{\text{D}}$ and
${V}_{\text{N}}$, up to their first-order derivatives – we can estimate them
in the ${{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma))$-norm via the interpolation
inequality (cf. (E.1))

$\left|{{\tilde{\rho}}}\right|_{\Gamma}^{(2)}\leqslant\varepsilon^{\alpha}\left|{{\tilde{\rho}}}\right|_{\Gamma}^{(2+\alpha)}+\frac{c_{8}}{\varepsilon^{2}}{{\max_{\Gamma}}}\left|{{\tilde{\rho}}}\right|,$
(66)

for some constant $c_{8}:=c_{8}(d)>0$. The above estimates follows from (64),
the maximum principle, and the fact that ${\tilde{\rho}}$ vanishes on the
fixed boundary $\Sigma$. The end estimate in the
${{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma))$-norm is achieved through (66) in
combination with the identities in (E.3) and of the first estimate in (E.4).
Now, while keeping these informations in mind, we employ the aforementioned
estimates to obtain

$\displaystyle\left|{\pazocal{{N}}[{V}_{\text{D}},{\tilde{\rho}}]}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\pazocal{{N}}[{V}_{\text{N}},{\tilde{\rho}}]}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\pazocal{{B}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}]}\right|_{\Gamma}^{(1+\alpha)}$
(67)
$\displaystyle\qquad\leqslant\tilde{c}_{1}t^{\eta}\left(\left|{{V}_{\text{D}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{V}_{\text{N}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{\tilde{\rho}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right)\left(\left|{{\tilde{\rho}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{\tilde{\rho}}_{\tau}}\right|^{\infty}_{{\Gamma;\,[0,t]}}\right)$
$\displaystyle\qquad\leqslant\tilde{c}_{1}t^{\eta}\left({\mathbb{W}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}](t)\right)^{2},$

where

${\mathbb{W}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}](t):=\left\|{{V}_{{\text{D},\text{N}}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|{\tilde{\rho}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]},$

for some constant $\tilde{c}_{1}>0$ and $\eta>0$. Here, the notation
${\mathbb{W}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}](t)$ is introduced
for convenience of later use.

The rest of the proof applies the method of successive approximation. To this
end, we will need the following estimate (cf. [14, Eq. 5.3, p. 131]) in our
argumentation further below. Let us consider the functions

${V}_{\text{D}1},{V}_{\text{D}2},{V}_{\text{N}1},{V}_{\text{N}2},{\tilde{\rho}}_{1},{\tilde{\rho}}_{2}\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\overline{\Omega}))$

such that

${\tilde{\rho}}_{1t},{\tilde{\rho}}_{2t}\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma)).$

Then, by estimate (67), we have

$\displaystyle\left|{\pazocal{{N}}[{V}_{\text{D}2},{\tilde{\rho}}_{2}]-\pazocal{{N}}[{V}_{\text{D}1},{\tilde{\rho}}_{1}]}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\pazocal{{N}}[{V}_{\text{N}2},{\tilde{\rho}}_{2}]-\pazocal{{N}}[{V}_{\text{N}1},{\tilde{\rho}}_{1}]}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}$
(68)
$\displaystyle\quad+\left|{\pazocal{{B}}[{V}_{\text{D}2},{V}_{\text{N}2},{\tilde{\rho}}_{2}]-\pazocal{{B}}[{V}_{\text{D}1},{V}_{\text{N}1},{\tilde{\rho}}_{1}]}\right|_{\Gamma}^{(1+\alpha)}$
$\displaystyle\quad\quad\leqslant\tilde{c}_{2}t^{\eta}\left\\{\left(|{V}_{\text{D}2}-{V}_{\text{D}1}|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+|{V}_{\text{N}2}-{V}_{\text{N}1}|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{\tilde{\rho}}_{2}-{\tilde{\rho}}_{1}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right)\right.$
$\displaystyle\quad\qquad\qquad\qquad\qquad\times\sum_{i=1,2}\left(|{\tilde{\rho}}_{i}|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{\tilde{\rho}}_{i\tau}}\right|^{\infty}_{{\Gamma;\,[0,t]}}\right)$
$\displaystyle\quad\qquad\qquad\qquad+\left(\left|{{\tilde{\rho}}_{2}-{\tilde{\rho}}_{1}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{\tilde{\rho}}_{2\tau}-{\tilde{\rho}}_{1\tau}}\right|^{\infty}_{{\Gamma;\,[0,t]}}\right)$
$\displaystyle\quad\qquad\qquad\qquad\qquad\quad\times\left.\sum_{i=1,2}\left(\left|{{V}_{\text{D}i}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{V}_{\text{N}i}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{{\tilde{\rho}}_{i}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right)\right\\},$
$\displaystyle\quad\quad\leqslant\tilde{c}_{2}t^{\eta}\left({\mathbb{W}}[{W}_{\text{D}},{W}_{\text{N}},{\tilde{\varrho}}](t)\right)\left\\{{\mathbb{W}}[{V}_{\text{D}1},{V}_{\text{N}1},{\tilde{\rho}}_{1}](t)+{\mathbb{W}}[{V}_{\text{D}2},{V}_{\text{N}2},{\tilde{\rho}}_{2}](t)\right\\},$

where ${W}_{\text{D}}={V}_{\text{D}2}-{V}_{\text{D}1}$,
${W}_{\text{N}}={V}_{\text{N}2}-{V}_{\text{N}1}$, and
${\tilde{\varrho}}={\tilde{\rho}}_{2}-{\tilde{\rho}}_{1}$, for some constant
$\tilde{c}_{2}>0$.

Now, we consider the initial approximant
$({V}_{\text{D}}^{(0)},{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)})$ which we
supposed to satisfy the linear system (50) with $F_{\text{D}}(y,t)=0$ and
$F_{\text{N}}(y,t)=0$, for $y\in\Omega$, $t>0$, and
$\psi(y,t)=-B_{0}\dfrac{\partial{{{{U}_{\text{N}\text{D}}^{0}}}}}{\partial{\nu_{\circ}}}$,
for $y\in\Gamma$, $t>0$. Meanwhile, the approximants
$\\{({V}_{\text{D}}^{(m+1)},{V}_{\text{N}}^{(m+1)},{\tilde{\rho}}^{(m+1)})\\}_{m}$,
for $m=0,1,\ldots$, are defined as solutions to the problem

$\left\\{\begin{aligned}
&\Delta{V}_{\text{D}}^{(m+1)}+\boldsymbol{\mathbf{q}}_{\text{D}}\cdot\nabla{\tilde{\rho}}^{(m+1)}+{Q}_{\text{D}}{\tilde{\rho}}^{(m+1)}=\pazocal{{N}}[{V}_{\text{D}}^{(m)},{\tilde{\rho}}^{(m)}],\quad
y\in\Omega,\ t>0,\\\
&\qquad\quad\left.{{V}_{\text{D}}^{(m+1)}}\right|_{y\in\Sigma}=0,\qquad\left.{{V}_{\text{D}}^{(m+1)}}\right|_{y\in\Gamma}=0,\\\\[5.0pt]
&\Delta{V}_{\text{N}}^{(m+1)}+\boldsymbol{\mathbf{q}}_{\text{N}}\cdot\nabla{\tilde{\rho}}^{(m+1)}+{Q}_{\text{N}}{\tilde{\rho}}^{(m+1)}=\pazocal{{N}}[{V}_{\text{N}}^{(m)},{\tilde{\rho}}^{(m)}],\quad
y\in\Omega,\ t>0,\\\
&\qquad\quad\left.{\dfrac{\partial{}}{\partial{\nu_{\circ}}}{V}_{\text{N}}^{(m+1)}}\right|_{y\in\Sigma}=0,\qquad\left.{{V}_{\text{N}}^{(m+1)}}\right|_{y\in\Gamma}=0,\\\\[5.0pt]
&\Delta{\tilde{\rho}}^{(m+1)}=0,\quad y\in\Omega,\
t>0\qquad\left.{{\tilde{\rho}}^{(m+1)}}\right|_{y\in\Sigma}=0,\\\
&\left.{\left(\dfrac{\partial{{\tilde{\rho}}^{(m+1)}}}{\partial{t}}+b_{\perp}\dfrac{\partial{{\tilde{\rho}}^{(m+1)}}}{\partial{\nu_{\circ}}}+\boldsymbol{\mathbf{b}}_{\parallel}\cdot\nabla_{\Gamma}{\tilde{\rho}}^{(m+1)}+h{\tilde{\rho}}^{(m+1)}+B_{0}\dfrac{\partial{{V}_{\text{N}\text{D}}}}{\partial{\nu_{\circ}}}\right)}\right|_{y\in\Gamma}\\\
&\qquad\qquad\qquad\qquad\
=\pazocal{{B}}[{V}_{\text{D}}^{(m)},{V}_{\text{N}}^{(m)},{\tilde{\rho}}^{(m)}]-B_{0}\dfrac{\partial{{{U}_{\text{N}\text{D}}^{0}}}}{\partial{\nu_{\circ}}},\\\
&\qquad{{\tilde{\rho}}^{(m+1)}}\big{|}_{y\in\Gamma,\,t=0}=0.\end{aligned}\right.$
(69)

Our goal now is to show that all approximants are defined on some interval
${{I}^{\star}}$ and that the sequence of sub-approximants
$\\{{V}_{\text{D}}^{(m)}\\}$, $\\{{V}_{\text{N}}^{(m)}\\}$, and
$\\{{\tilde{\rho}}^{(m)}\\}$ converge.

By virtue of Theorem 7.1, the initial approximation
$({V}_{\text{D}}^{(0)},{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)})$ is defined
for all $t\in[0,T]$ and satisfies the estimates

${\mathbb{W}}[{V}_{\text{D}}^{(0)},{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)}](t)\leqslant\tilde{c}_{3}\left\|{(U^{0}_{\text{D}},U^{0}_{\text{N}})}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}=:\tilde{c}_{3}\left\|{{{U}_{{\text{D},\text{N}}}^{0}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}},\quad\text{for
all $t\leqslant T$},$ (70)

for some constant $\tilde{c}_{3}>0$. Meanwhile, to get an estimate for the
first approximant
$({V}_{\text{D}}^{(1)},{V}_{\text{N}}^{(1)},{\tilde{\rho}}^{(1)})$, we apply
Theorem 7.1 to (69) with $m=0$, to obtain

$\displaystyle{\mathbb{W}}[{V}_{\text{D}}^{(1)},{V}_{\text{N}}^{(1)},{\tilde{\rho}}^{(1)}](t)$
$\displaystyle\leqslant\tilde{c}_{3}\left\\{\left|{\pazocal{{N}}[{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)}]}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\pazocal{{N}}[{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)}]}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}\right.$
(71)
$\displaystyle\qquad\quad+\left.\left|{\pazocal{{B}}[{V}_{\text{D}}^{(0)},{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)}]}\right|_{\Gamma}^{(1+\alpha)}+\left\|{{{U}_{{\text{D},\text{N}}}^{0}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right\\}$
$\displaystyle\leqslant\tilde{c}_{3}\left\\{\tilde{c}_{1}t^{\eta}\left({\mathbb{W}}[{V}_{\text{D}}^{(0)},{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)}](t)\right)^{2}+\left\|{{{U}_{{\text{D},\text{N}}}^{0}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right\\}$
$\displaystyle\leqslant\tilde{c}_{3}^{3}\tilde{c}_{1}t^{\eta}\left(\left\|{{{U}_{{\text{D},\text{N}}}^{0}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right)^{2}+\tilde{c}_{3}\left\|{{{U}_{{\text{D},\text{N}}}^{0}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}$
$\displaystyle=:{{L}}(t),$

where the second inequality follows from (67), while the third one is due to
estimate (70).

We introduce the following notations:

${W}_{\text{D}}^{(m+1)}={V}_{\text{D}}^{(m+1)}-{V}_{\text{D}}^{(m)},\qquad{W}_{\text{N}}^{(m+1)}={V}_{\text{N}}^{(m+1)}-{V}_{\text{N}}^{(m)},\qquad{\tilde{\varrho}}^{(m+1)}={\tilde{\rho}}^{(m+1)}-{\tilde{\rho}}^{(m)}.$

Then, by telescoping sums, we see that

${V}_{\text{D}}^{(m+1)}={V}_{\text{D}}^{(1)}+\sum_{j=1}^{m}{W}_{\text{D}}^{(j+1)},\quad{V}_{\text{D}}^{(m+1)}={V}_{\text{D}}^{(1)}+\sum_{j=1}^{m}{W}_{\text{D}}^{(j+1)},\quad{\tilde{\rho}}^{(m+1)}={\tilde{\rho}}^{(1)}+\sum_{j=1}^{m}{\tilde{\varrho}}^{(j+1)}.$

Moreover, by (69), ${W}_{\text{D}}^{(m+1)}$, ${W}_{\text{N}}^{(m+1)}$, and
${\tilde{\varrho}}^{(m+1)}$ satisfy

$\left\\{\begin{aligned}
&\Delta{W}_{\text{D}}^{(m+1)}+\boldsymbol{\mathbf{q}}_{\text{D}}\cdot\nabla{\tilde{\rho}}^{(m+1)}+{Q}_{\text{D}}{\tilde{\rho}}^{(m+1)}\\\
&\qquad\qquad\qquad\qquad\ \
=\pazocal{{N}}[{W}_{\text{D}}^{(m)},{\tilde{\rho}}^{(m)}]-\pazocal{{N}}[{W}_{\text{D}}^{(m-1)},{\tilde{\rho}}^{(m-1)}],\quad
y\in\Omega,\ t>0,\\\
&\qquad\quad\left.{{W}_{\text{D}}^{(m+1)}}\right|_{y\in\Sigma}=0,\qquad\left.{{W}_{\text{D}}^{(m+1)}}\right|_{y\in\Gamma}=0,\\\\[10.00002pt]
&\Delta{W}_{\text{N}}^{(m+1)}+\boldsymbol{\mathbf{q}}_{\text{N}}\cdot\nabla{\tilde{\rho}}^{(m+1)}+{Q}_{\text{N}}{\tilde{\rho}}^{(m+1)}\\\
&\qquad\qquad\qquad\qquad\ \
=\pazocal{{N}}[{V}_{\text{N}}^{(m)},{\tilde{\rho}}^{(m)}]-\pazocal{{N}}[{V}_{\text{N}}^{(m-1)},{\tilde{\rho}}^{(m-1)}],\quad
y\in\Omega,\ t>0,\\\
&\qquad\quad\left.{\dfrac{\partial{}}{\partial{\nu_{\circ}}}{W}_{\text{N}}^{(m+1)}}\right|_{y\in\Sigma}=0,\qquad\left.{{W}_{\text{N}}^{(m+1)}}\right|_{y\in\Gamma}=0,\\\\[10.00002pt]
&\Delta{\tilde{\rho}}^{(m+1)}=0,\quad y\in\Omega,\
t>0\qquad\left.{{\tilde{\rho}}^{(m+1)}}\right|_{y\in\Sigma}=0,\\\
&\left.{\left(\dfrac{\partial{{\tilde{\rho}}^{(m+1)}}}{\partial{t}}+b_{\perp}\dfrac{\partial{{\tilde{\rho}}^{(m+1)}}}{\partial{\nu_{\circ}}}+\boldsymbol{\mathbf{b}}_{\parallel}\cdot\nabla_{\Gamma}{\tilde{\rho}}^{(m+1)}+h{\tilde{\rho}}^{(m+1)}+B_{0}\dfrac{\partial{{V}_{\text{N}\text{D}}}}{\partial{\nu_{\circ}}}\right)}\right|_{y\in\Gamma}\\\
&\qquad\qquad\qquad\qquad\
=\pazocal{{B}}[{W}_{\text{D}}^{(m)},{W}_{\text{N}}^{(m)},{\tilde{\rho}}^{(m)}]-\pazocal{{B}}[{W}_{\text{D}}^{(m-1)},{W}_{\text{N}}^{(m-1)},{\tilde{\rho}}^{(m-1)}],\\\
&\qquad{{\tilde{\rho}}^{(m+1)}}\big{|}_{y\in\Gamma,\,t=0}=0.\end{aligned}\right.$
(72)

We again utilize Theorem 7.1, and make use of the estimate (68), to obtain

$\displaystyle{\mathbb{W}}[{W}_{\text{D}}^{(m+1)},{W}_{\text{N}}^{(m+1)},{\tilde{\varrho}}^{(m+1)}](t)$
(73)
$\displaystyle\qquad\leqslant\tilde{c}_{2}\tilde{c}_{3}t^{\eta}\left(\sum_{j=m,m-1}{\mathbb{W}}[{V}_{\text{D}}^{(j)},{V}_{\text{N}}^{(j)},{\tilde{\rho}}^{(j)}](t)\right){\mathbb{W}}[{W}_{\text{D}}^{(m)},{W}_{\text{N}}^{(m)},{\tilde{\varrho}}^{(m)}](t).$

We are now in the position to determine the maximum time ${{t}^{\star}}$ so
that the claims we state in Theorem 4.4 are valid. First, we choose
${{t}^{\star}}$ so that our first approximant satisfy

$({{t}^{\star}})^{\eta}{\mathbb{W}}[{V}_{\text{D}}^{(1)},{V}_{\text{N}}^{(1)},{\tilde{\rho}}^{(1)}](t)\leqslant({{t}^{\star}})^{\eta}{{L}}(t)\leqslant{{M}_{0}}$

and such that ${{M}^{\star}}=2\tilde{c}_{2}\tilde{c}_{3}{{M}_{0}}<M_{1}<1$ for
some fixed constant $M_{1}>0$. Clearly, the initial approximant
${\mathbb{W}}[{V}_{\text{D}}^{(0)},{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)}](t)$
also satisfy the estimate
$({{t}^{\star}})^{\eta}{\mathbb{W}}[{V}_{\text{D}}^{(0)},{V}_{\text{N}}^{(0)},{\tilde{\rho}}^{(0)}](t)\leqslant{{M}_{0}}$.
We then assume that
$({V}_{\text{D}}^{(j)},{V}_{\text{N}}^{(j)},{\tilde{\rho}}^{(j)})$,
$j=0,1,\ldots,m$, are defined for $t\in{{I}^{\star}}$ and satisfy the estimate

$({{t}^{\star}})^{\eta}{\mathbb{W}}[{V}_{\text{D}}^{(j)},{V}_{\text{N}}^{(j)},{\tilde{\rho}}^{(j)}](t)\leqslant{{M}_{0}},\quad\text{for
all $j=0,1,\ldots,m$}.$ (74)

By Theorem 7.1, it follows that the $(m+1)$th approximant given by the triplet
$({V}_{\text{D}}^{(m+1)},{V}_{\text{N}}^{(m+1)},{\tilde{\rho}}^{(m+1)})$ is
also defined in ${{I}^{\star}}$. Combining estimates (73) and (74), we get

$\displaystyle{\mathbb{W}}[{W}_{\text{D}}^{(m+1)},{W}_{\text{N}}^{(m+1)},{\tilde{\varrho}}^{(m+1)}]({{t}^{\star}})$
$\displaystyle\leqslant
2\tilde{c}_{2}\tilde{c}_{3}{{M}_{0}}{\mathbb{W}}[{W}_{\text{D}}^{(m)},{W}_{\text{N}}^{(m)},{\tilde{\varrho}}^{(m)}]({{t}^{\star}})$
$\displaystyle\leqslant{{M}^{\star}}{\mathbb{W}}[{W}_{\text{D}}^{(m)},{W}_{\text{N}}^{(m)},{\tilde{\varrho}}^{(m)}]({{t}^{\star}}).$

Taking the sum of above inequality with respect to $j$ from $0$ to $m$ yields

$\displaystyle{{S}}_{m+1}$
$\displaystyle:=\sum_{j=0}^{m}{\mathbb{W}}[{W}_{\text{D}}^{(j+1)},{W}_{\text{N}}^{(j+1)},{\tilde{\varrho}}^{(j+1)}]({{{t}}})$
$\displaystyle\quad\leqslant{{M}^{\star}}\sum_{j=0}^{m}{\mathbb{W}}[{W}_{\text{D}}^{(j)},{W}_{\text{N}}^{(j)},{\tilde{\varrho}}^{(j)}]({{{t}}})$
$\displaystyle\quad={{M}^{\star}}{{S}}_{m+1}+{{M}^{\star}}\left({\mathbb{W}}[{W}_{\text{D}}^{(0)},{W}_{\text{N}}^{(0)},{\tilde{\varrho}}^{(0)}]({{{t}}})-{\mathbb{W}}[{W}_{\text{D}}^{(m+1)},{W}_{\text{N}}^{(m+1)},{\tilde{\varrho}}^{(m+1)}]({{{t}}})\right)$
$\displaystyle\quad\leqslant{{M}^{\star}}{{S}}_{m+1}+{{M}^{\star}}{\mathbb{W}}[{W}_{\text{D}}^{(0)},{W}_{\text{N}}^{(0)},{\tilde{\varrho}}^{(0)}]({{{t}}})$
$\displaystyle\quad\leqslant{{M}^{\star}}{{S}}_{m+1}+{{L}}(t),$

where the last inequality follows from (71) and the fact that
${{M}^{\star}}<1$. The above estimate clearly shows that the sums
${{S}}_{m+1}$, $m=0,1,\ldots$, are uniformly bounded by ${{L}}(t)$. This
information implies that the series sums
$\sum_{j=1}^{\infty}{\mathbb{W}}[{W}_{\text{D}}^{(j+1)},{W}_{\text{N}}^{(j+1)},{\tilde{\varrho}}^{(j+1)}](t)$
converges for all $t\in{{I}^{\star}}$ and that
$\\{({V}_{\text{D}}^{(m+1)},{V}_{\text{N}}^{(m+1)},{\tilde{\rho}}^{(m+1)})\\}_{m}$,
$m=0,1,\ldots$, also satisfy inequality condition (74).

Finally, passing to the limit $m\to\infty$, we conclude that

$({V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}})=\lim_{m\to\infty}({V}_{\text{D}}^{(m)},{V}_{\text{N}}^{(m)},{\tilde{\rho}}^{(m)})$

is a solution to (49) satisfying the estimate

${\mathbb{W}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}](t)\leqslant
c{{L}}(t)\leqslant
c\left\|{{{U}_{{\text{D},\text{N}}}^{0}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}},\quad\forall
t\in{{I}^{\star}},$ (75)

for some constant $c>0$.

To complete the proof of Theorem 4.4, we need to revert to the original
unknown functions that resolve the transformed problem (30). As stated in
Section 6, we have the relations
$U_{\text{D}}={V}_{\text{D}}+U^{0}_{\text{D}}$ and
$U_{\text{N}}={V}_{\text{N}}+U^{0}_{\text{N}}$. These functions, along with
${\tilde{\rho}}$, solve (30). Additionally, estimate (37) is derived from
(75), in conjunction with estimates (43) and (45). Furthermore, since the
difference between two solutions must also satisfy estimate (37) with a zero
right-hand side, the solution $(U_{\text{D}},U_{\text{N}},{\tilde{\rho}})$ is
unique. Finally, to obtain $\rho$, we set
$\rho={\tilde{\rho}}\big{|}_{\Gamma}$ and then determine the free boundary
using the description given in (18).

$\Gamma(t):=\left\\{x\in\mathbb{R}^{d}\mid
x=\xi+\rho(\xi,t)\boldsymbol{\mathbf{N}}(\xi),\
\xi\in\Gamma\right\\},\qquad\text{for $t\in{{I}^{\star}}$}.$

∎

To conclude this work, we remark – as alluded at an earlier part of this note
– that the main result given by Theorem 3.1 follows from Theorem 4.4 which
asserts the local-in-time solvability of problem (30). Indeed, recalling the
change of variables $Y^{-1}(y,t)$ from (28), we can transform the fixed domain
$\Omega$ to the moving domain $\Omega(t)$, for $t\in{{I}^{\star}}$, and
retrieve the functions
$u_{\text{D}}(x,t)=U_{\text{D}}(y,t)=u_{\text{D}}(Y^{-1}(y,t),t)$ and
$u_{\text{N}}(x,t)=U_{\text{N}}(y,t)=u_{\text{N}}(Y^{-1}(y,t),t)$ (cf. (C.1))
as the pair of solutions to (1). In addition, estimate (19) immediately
follows from (37), confirming our assertion. This verifies Theorem 3.1.

## Data availability

No data was used for the research described in the article.

## Acknowledgments

The work of JFTR is partially supported by JST CREST Grant Number JPMJCR2014
and JSPS KAKENHI Grant Number JP23K13012. The work of MK is supported by JSPS
KAKENHI Grant Numbers JP20KK0058 and JP20H01812.

##  Appendix A Comparison of normal derivatives

In this appendix, we examine the normal derivatives of two functions that
fulfill a pure Dirichlet problem and a mixed Dirichlet-Neumann problem on the
boundary of axisymmetric domains. Throughout the section, $\rho>0$ denotes the
radius of either a circle or a sphere. In this section, $f$ and $g$ are two
fixed positive-valued scalar functions defined on $\mathbb{R}^{2}$.

###  A.1 The case of concentric circles

For axisymmetric case in two spatial dimensions, the Laplace equation becomes

$\Delta
u(\rho)=\frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{d{u}}{d\rho}\right)=0\qquad\text{or}\qquad\Delta
u(\rho)=\frac{d}{d\rho}\left(\rho\frac{d{u}}{d\rho}\right)=0.$

Integrating once, we get $d{u(\rho)}/{d\rho}={a}/{\rho}$, and then again,
$u=a\log\rho+b$, for some unknowns $a$ and $b$.

Let $B_{r}:=B(\boldsymbol{\mathbf{0}},r)$ be a circle with radius $r>0$. Let
$r^{\star}\in(0,R)$ and consider the PDE system

$\Delta{{u}^{\star}}=0\quad\text{in}\
{\Omega}^{\star}:=B_{R}\setminus\overline{B}_{{r}^{\star}},\qquad{{u}^{\star}}=f>0\quad\text{on}\
\Sigma,\qquad{{u}^{\star}}=0\quad\text{on}\
{\Gamma}^{\star}:=\partial{B_{{r}^{\star}}}.$ (A.1)

Then, we can compute the exact solution as

${u}^{\star}(\rho)=\frac{f\log{(\rho/{r}^{\star})}}{\log{(R/{r}^{\star})}},\qquad\rho\in({r}^{\star},R).$
(A.2)

Note that for any $\rho\in({r}^{\star},R]$, $u(\rho)$ is positive. Now,
differentiating this function with respect to $\rho$, we get

$\frac{\partial}{\partial\nu_{R}}{u}^{\star}(\rho)=\frac{\partial}{\partial\left|{x}\right|}{u}^{\star}(\left|{x}\right|)=\frac{\partial}{\partial\rho}{u}^{\star}(\rho)=\frac{f}{\rho\log{(R/{r}^{\star})}},$
(A.3)

where $\nu_{R}$ is outward unit normal vector to $\partial B_{R}$. At
$\rho=R$, we have

$\frac{\partial}{\partial\nu_{R}}{u}^{\star}(R)=\frac{f}{R\log{(R/{r}^{\star})}}.$

We define

$g:=\frac{f}{R\log{(R/{r}^{\star})}}>0,\qquad 0<{r}^{\star}<R.$

Let us consider Equations (8) and (9) with
$\Omega=B_{R}\setminus\overline{B}_{r}$, $\Gamma=\partial{B_{R}}$, and
$\Sigma=\partial{B_{R}}$, where $0<r<R$. In view of (A.2), it can easily be
verified that the exact solutions to these systems of PDEs are respectively
given by

$u_{\text{D}}(\rho)=\frac{f\log{(\rho/r)}}{\log{(R/r)}}\qquad\text{and}\qquad
u_{\text{N}}(\rho)=\frac{f\log{(\rho/r)}}{\log{(R/{r}^{\star})}},\qquad\rho\in[r,R].$

Now, on $\Gamma$, straightforward computation of $V_{n}$ gives us

$\displaystyle V_{n}$
$\displaystyle=-\left(\dfrac{\partial{u_{\text{D}}}}{\partial\nu}-\dfrac{\partial{u_{\text{N}}}}{\partial\nu}\right)$
$\displaystyle=-\left.\left[\frac{f}{\log{(R/r)}}\left(\frac{1}{\rho}\right)-\frac{f}{\log{(R/{r}^{\star})}}\left(\frac{1}{\rho}\right)\right]\right|_{\rho=r}$
$\displaystyle=-\frac{f}{r{\log{(R/r)\log{(R/{r}^{\star})}}}}\left[\log{(R/{r}^{\star})}-\log{(R/r)}\right]$
$\displaystyle=-\underbrace{\frac{f}{r{\log{(R/r)\log{(R/{r}^{\star})}}}}}_{>0}\underbrace{\left[\log{\left(\dfrac{r}{{r}^{\star}}\right)}\right]}_{:=L_{1}(r)}.$

Observe that the sign of $V_{n}$ only depends on the relations between $r$ and
${r}^{\star}$. If we want $V_{n}$ to be negative, we need $L_{1}(r)>0$ to be
positive, and this happens when $r>{r}^{\star}$.

Now, motivated by the previous discussion, we examine the existence and
uniqueness of solution to the following initial value problem: for a given
$T>0$,

$\left\\{\begin{aligned}
r^{\prime}(t)&=-\dfrac{f\left(\log{r(t)}-\log{{r}^{\star}}\right)}{r(t)\log{(R/{r}^{\star})}\left(\log{R}-\log{r(t)}\right)}=:F_{1}(t,r(t)),&&\text{for
}0<t<T,\\\ r(0)&=r_{0},&&r_{0}>{r}^{\star},\end{aligned}\right.$ (A.4)

for some given $r_{0}\in({r}^{\star},R)$.

Given the techical assumption state in Remark 1.3, we show – although it is
trivial – that the initial value problem (A.4) admits a unique solution. We
prove the existence using Peano’s theorem [57, Thm. 2.1, p. 10] and its
uniqueness via Picard-Lindelöf theorem [57, Thm. 1.1, p. 8]. First, we observe
that for $t>0$, $0<{r}^{\star}\leqslant r(t)<R$ (see the proof of Proposition
A.2), and for $T>0$, $F_{1}$ is continuous for all $t\in[0,T]$ and
$\left|{r(t)-r_{0}}\right|\leqslant R-{r}^{\star}$. Moreover, $r$ is uniformly
bounded in $[0,T]$, and therefore, $F_{1}(t,r(t))$ is bounded in
$[0,T]\times[{r}^{\star},R_{0}]$ for some ${r}^{\star}<r_{0}\leqslant
R_{0}<R$. In fact, we see that

$\left|{F_{1}(t,r(t))}\right|\leqslant\dfrac{\|{f}\|_{L^{\infty}(\Sigma)}}{{r}^{\star}\log{(R/R_{0})}}=:K_{0},\qquad\text{for
all $(t,r(t))\in[0,T]\times[{r}^{\star},R_{0}]$}.$

By Peano’s theorem [57, Thm. 2.1, p. 10], (A.4) possesses at least one
solution $r=r(t)$ on $[0,T_{0}]$ where
$T_{0}=\min\\{T,(R-{r}^{\star})/K_{0}\\}$.

######  Remark A.1.

The choice $T_{0}=\min\\{T,(R-{r}^{\star})/K_{0}\\}$ for the existence of
solution to (A.4) is natural. On the one hand, the requirement $T_{0}\leqslant
T$ is necessary. On the other hand, the requirement
$T_{0}<(R-{r}^{\star})/K_{0}$ is due to the fact that if $r=r(t)$ is a
solution of (A.4) on $t\in[0,T_{0}]$, then
$\left|{r^{\prime}(t)}\right|\leqslant K_{0}$ implies that
$\left|{r(t)-r_{0}}\right|\leqslant R-{r}^{\star}\leqslant K_{0}(t-0)=K_{0}t$.
However, we note that $\left|{r(t)-r_{0}}\right|\leqslant R-{r}^{\star}$.
Thus, we require $t\leqslant(R-{r}^{\star})/K_{0}$ so that the previous
inequality holds.

To prove that the solution to (A.4) is unique, we will prove that $F_{1}$ is
uniformly Lipschitz continuous with respect to $r$. Let us consider

$F_{1}(r_{j}):=F_{1}(t,r_{j}(t))=-\frac{f\left(\log{r_{j}(t)}-\log{{r}^{\star}}\right)}{r_{j}(t)\log{(R/{r}^{\star})}\left(\log{R}-\log{r_{j}(t)}\right)},\quad\text{for
$j=1,2$, and $t>0$},$

where ${r}^{\star}<r_{j}\leqslant R_{0}$ for $j=1,2$, for some $0<R_{0}<R$
(refer to the third point in Remark 1.3). Hence,we have

$\displaystyle\left|{F_{1}(r_{1})-F_{1}(r_{2})}\right|$
$\displaystyle=\left|{\frac{f\left(\log{r_{1}}-\log{{r}^{\star}}\right)}{r_{1}\log{(R/{r}^{\star})}\left(\log{R}-\log{r_{1}}\right)}-\frac{f\left(\log{r_{2}}-\log{{r}^{\star}}\right)}{r_{2}\log{(R/{r}^{\star})}\left(\log{R}-\log{r_{2}}\right)}}\right|$
$\displaystyle\leqslant\dfrac{\|{f}\|_{L^{\infty}(\Sigma)}}{\log{(R/{r}^{\star})}}\left|{\frac{r_{2}\left(\log{R}-\log{r_{2}}\right)\left(\log{r_{1}}-\log{{r}^{\star}}\right)-r_{1}\left(\log{R}-\log{r_{1}}\right)\left(\log{r_{2}}-\log{{r}^{\star}}\right)}{r_{1}r_{2}\left(\log{R}-\log{r_{1}}\right)\left(\log{R}-\log{r_{2}}\right)}}\right|$
$\displaystyle=:\dfrac{\|{f}\|_{L^{\infty}(\Sigma)}}{\log{(R/{r}^{\star})}}\left|{\frac{\textsf{n}(r_{1},r_{2})}{\textsf{d}(r_{1},r_{2})}}\right|,$

where

$\displaystyle\textsf{n}(r_{1},r_{2})$
$\displaystyle=(\log{R})r_{2}\log{r_{1}}-r_{2}\log{r_{1}}\log{r_{2}}-\log{{r}^{\star}}(\log{R})r_{2}+(\log{{r}^{\star}})r_{2}\log{r_{2}}$
$\displaystyle\qquad-(\log{R})r_{1}\log{r_{2}}+r_{1}\log{r_{2}}\log{r_{1}}+\log{{r}^{\star}}(\log{R})r_{1}-(\log{{r}^{\star}})r_{1}\log{r_{1}}$
$\displaystyle=\log{R}(r_{2}\log{r_{1}}-r_{1}\log{r_{2}})-\log{r_{1}}\log{r_{2}}(r_{2}-r_{1})$
$\displaystyle\qquad+\log{{r}^{\star}}(r_{2}\log{r_{2}}-r_{1}\log{r_{1}})+\log{{r}^{\star}}\log{R}(r_{1}-r_{2}),$
$\displaystyle\textsf{d}(r_{1},r_{2})$
$\displaystyle=r_{1}r_{2}\left(\log{R}-\log{r_{1}}\right)\left(\log{R}-\log{r_{2}}\right).$

For $r_{j}>0$ such that ${r}^{\star}<r_{j}\leqslant R_{0}<R$ for $j=1,2$, the
denominator d is bounded below away from zero; i.e.,

$\left|{\textsf{d}(r_{1},r_{2})}\right|\geqslant\left[{r}^{\star}\left(\log{R}-\log{R_{0}}\right)\right]^{2}=:K_{1}.$

Now, on the other hand, we claim that $\textsf{n}(r_{1},r_{2})$ can be bounded
by $\left|{r_{1}-r_{2}}\right|$; that is, there exists a constant $K_{2}>0$
such that

$\left|{\textsf{n}(r_{1},r_{2})}\right|\leqslant
K_{2}\left|{r_{1}-r_{2}}\right|,$

where $K_{2}$ does not depends on $r_{1}$ and $r_{2}$.

Let us note that

$\displaystyle\left|{\textsf{n}(r_{1},r_{2})}\right|$
$\displaystyle=\left|{\log{R}}\right|\left|{r_{2}\log{r_{1}}-r_{1}\log{r_{2}}}\right|+\left|{\log{r_{1}}}\right|\left|{\log{r_{2}}}\right|\left|{r_{2}-r_{1}}\right|$
$\displaystyle\qquad+\left|{\log{{r}^{\star}}}\right|\left|{r_{2}\log{r_{2}}-r_{1}\log{r_{1}}}\right|+\left|{\log{{r}^{\star}}}\right|\left|{\log{R}}\right|\left|{r_{1}-r_{2}}\right|.$

We focus on proving that $\left|{r_{2}\log{r_{1}}-r_{1}\log{r_{2}}}\right|$
can be bounded by $\left|{r_{1}-r_{2}}\right|$. The rest can be shown in a
similar manner. So, let us consider the function $\textsf{f}(r)=r\log{r}$. We
note that f is continuous on $[r_{1},r_{2}]$ and differentiable on
$(r_{1},r_{2})$. Hence, by the mean value theorem, there exists
$\bar{r}\in(r_{1},r_{2})$ such that

$\log{\bar{r}}+1=\textsf{f}^{\prime}(\bar{r})=\dfrac{\textsf{f}(r_{2})-\textsf{f}(r_{1})}{r_{2}-r_{1}}=\dfrac{r_{2}\log{r_{2}}-r_{1}\log{r_{1}}}{r_{2}-r_{1}},\qquad(0<{r}^{\star}<r_{j}\leqslant
R_{0}<R,\ j=1,2).$

Because $r_{1},r_{2}>0$, $\left|{\log{\bar{r}}+1}\right|$ is clearly bounded
and does not depend on $r_{1}$ and $r_{2}$. Therefore,
$\left|{r_{2}\log{r_{2}}-r_{1}\log{r_{1}}}\right|=c\left|{r_{2}-r_{1}}\right|$,
for some constant of $c>0$.

Combining the bounds on the numerator and denominator, we get

$\left|{\textsf{f}(r_{2})-\textsf{f}(r_{1})}\right|\leq\frac{K_{2}}{K_{1}}\left|{u_{1}-u_{2}}\right|,$

which shows that $\textsf{f}(r)$ is uniformly Lipschitz continuous. Finally,
applying Picard-Lindelöf theorem [57, Thm. 1.1, p. 8], we deduce that the
solution to (A.4) is unique.

###  A.2 The case of concentric spheres

For axisymmetric case in three spatial dimensions, the Laplace equation
becomes

$\Delta
u(\rho)=\frac{1}{\rho^{2}}\frac{d}{d\rho}\left(\rho^{2}\frac{du}{d\rho}\right)=0\qquad\text{or}\qquad\Delta
u(\rho)=\frac{d}{d\rho}\left(\rho^{2}\frac{du}{d\rho}\right)=0.$

Integrating once, we get ${du(\rho)}/{d\rho}={a}/{\rho^{2}}$, and then again,
$u(\rho)=-{a}/{\rho}+b$, for some unknowns $a$ and $b$.

Let $B_{r}:=B(\boldsymbol{\mathbf{0}},r)$ be a sphere with radius $r>0$. Let
$r^{\star}\in(0,R)$ and consider again the PDE system (A.1) but with
${\Omega}^{\star}:=B_{R}\setminus\overline{B}_{{r}^{\star}}$ and
${\Gamma}^{\star}:=\partial{B_{{r}^{\star}}}$ described by the given spheres.
Then, we can compute the exact solution as

${u}^{\star}(\rho)=\frac{fR}{R-{r}^{\star}}\left(1-{r}^{\star}\rho^{-1}\right),\qquad\rho\in({r}^{\star},R).$
(A.5)

Observe, that for any $\rho\in({r}^{\star},R]$, $u(\rho)$ is positive.

Now, we differentiate ${u}^{\star}$ with respect to $\rho$, to obtain

$\frac{\partial}{\partial\nu_{R}}{u}^{\star}(\rho)=\frac{\partial}{\partial\left|{x}\right|}{u}^{\star}(\left|{x}\right|)=\frac{\partial}{\partial\rho}{u}^{\star}(\rho)=\frac{f}{\rho\log{(R/{r}^{\star})}},$
(A.6)

where $\nu_{R}$ is outward unit normal vector to $\partial B_{R}$. At
$\rho=R$, we have

$\frac{\partial}{\partial\nu_{R}}{u}^{\star}(R)=\frac{fR{r}^{\star}}{R-{r}^{\star}}\rho^{-2}\Big{|}_{\rho=R}=\frac{f{r}^{\star}}{R(R-{r}^{\star})}.$

We define

$g:=\frac{f{r}^{\star}}{R(R-{r}^{\star})}>0,\qquad 0<{r}^{\star}<R.$

Let us again consider Equations (8) and (9), but now in three dimensions, with
$\Omega=B_{R}\setminus\overline{B}_{r}$, $\Gamma=\partial{B_{R}}$, and
$\Sigma=\partial{B_{R}}$, where $0<r<R$. In view of (A.5), it can easily be
shown that the exact solutions to these systems of PDEs are respectively given
by

$u_{\text{D}}(\rho)=\frac{fR}{R-r}\left(1-r\rho^{-1}\right)\qquad\text{and}\qquad
u_{\text{N}}(\rho)=\frac{fR{r}^{\star}}{R-{r}^{\star}}\left(r^{-1}-\rho^{-1}\right),\qquad\rho\in[r,R].$

Now, computing $V_{n}$ on $\Gamma$ yields the following

$\displaystyle V_{n}$
$\displaystyle=-\left(\dfrac{\partial{u_{\text{D}}}}{\partial\nu}-\dfrac{\partial{u_{\text{N}}}}{\partial\nu}\right)$
$\displaystyle=-\left.\left[\frac{fR}{R-r}\left(\frac{r}{\rho^{2}}\right)-\frac{fR{r}^{\star}}{R-{r}^{\star}}\left(\frac{1}{\rho^{2}}\right)\right]\right|_{\rho=r}$
$\displaystyle=-\frac{fR}{r}\left[\frac{1}{R-r}-\frac{{r}^{\star}}{(R-{r}^{\star})r}\right]$
$\displaystyle=-\frac{fR}{r}\left[\frac{(R-{r}^{\star})r-(R-r){r}^{\star}}{r(R-r)(R-{r}^{\star})}\right]$
$\displaystyle=-\underbrace{\frac{fR^{2}}{r^{2}(R-r)(R-{r}^{\star})}}_{>0}\underbrace{(r-{r}^{\star})}_{=:L_{2}(r)}.$

Notice that, as expected, the sign of $V_{n}$ only depends on the relations
between $r$ and ${r}^{\star}$. If we want $V_{n}$ to be negative, we obviously
need $L_{2}(r)>0$ to be positive, and this occurs when $r>{r}^{\star}$.

Similar to the previous subsection, one can prove the existence and uniqueness
of solution to the following initial value problem

$\left\\{\begin{aligned}
r^{\prime}(t)&=-\frac{fR^{2}(r(t)-{r}^{\star})}{r(t)^{2}(R-r(t))(R-{r}^{\star})}=:F_{2}(t,r(t)),&&\text{for
}0<t<T,\\\ r(0)&=r_{0},&&r_{0}>{r}^{\star},\end{aligned}\right.$ (A.7)

for some given $r_{0}\in({r}^{\star},R)$. To prove that the right-hand side of
(A.7) is uniformly Lipshitz continuous with respect to $r$, we can again apply
the mean value theorem. For this purpose, we can utilize the derivative of $G$
with respect to $r$ given by

$F_{2}^{\prime}(\cdot,r)=-\dfrac{fR^{2}(2r^{2}+{r}^{\star}(2R-3r)-Rr)}{r^{3}(R-{r}^{\star})(R-r)^{2}},$

and show that this expression is bounded within the interval
$[{r}^{\star},R)$.

In view of the existence of unique solution to (A.4) and (A.7), we formally
have the following proposition.

######  Proposition A.2.

The unique solutions to (A.4) and (A.7) satisfy $r^{\prime}(t)<0$, for all
$t>0$.

###### Proof.

From previous discussion, we have deduced that $V_{n}$ remains negative for
$r>{r}^{\star}$. We want to prove then that $r(t)\geqslant{r}^{\star}$, for
all $t>0$, and we claim that $r(t)\to{r}^{\star}$ as $t\to T$ where
$T=\infty$. In other words, $r(t)$ will not reach ${r}^{\star}$ in finite
time. We will prove the latter claim via a contradiction. That is, we first
suppose that there exists a $T<\infty$ such that $r(T)={r}^{\star}$ and show
that this will lead to a contradiction. To do this, we first comment that
$F_{i}(\cdot,r)\in{{C}}^{\infty}(0,R)$, where $F_{i}$, $i=1,2$, are the
rational functions given in (A.4) and (A.7), respectively. Hence, via Taylor
expansion, we can write $r^{\prime}$ (after some transformation normalizing
the coefficient) as follows

$r^{\prime}(t)=({r}^{\star}-r(t))s_{i}(r(t)),$

where $s_{i}$ are some functions over $r$ such that $s_{i}({r}^{\star})>0$ for
$i=1,2$. The functions $s_{i}$, $i=1,2$, can be expressed in the form

$s_{i}(r)=\frac{r^{\prime}(t)}{{r}^{\star}-r(t)}=\log\left(\frac{1}{r(t)-{r}^{\star}}\right).$

Now, if $r\in{{C}}^{1}([0,T))\cap{{C}}^{0}([0,T])$ and $r(T)={r}^{\star}$ for
$0\leqslant a<T-\varepsilon<T<\infty$ for some constant $a$ and
$\varepsilon>0$ then

$\int_{a}^{T-\varepsilon}s_{i}(r(t))\
{d}{t}=\left.\log\left(\frac{1}{r(t)-{r}^{\star}}\right)\right|_{a}^{T-\varepsilon}=\log\left(\frac{r(a)-{r}^{\star}}{r(T-\varepsilon)-{r}^{\star}}\right).$

Notice that

$\lim_{\varepsilon\to 0}\int_{a}^{T-\varepsilon}s_{i}(r(t))\
{d}{t}=\lim_{\varepsilon\to
0}\log\left(\frac{r(a)-{r}^{\star}}{r(T-\varepsilon)-{r}^{\star}}\right)=\infty,$

which is a contradiction. ∎

##  Appendix B Mullins-Sekerka analysis for 2D axisymmetric domains

In this appendix, we examine the main system (1) on axisymmetric case in two
dimension through Mullins-Sekerka analysis given that $f$ and $g$ are
positive-valued functions. Let $r_{\Sigma}>0$ and $\rho^{0}>0$ be given. We
let the exterior (fixed) and the interior (free) boundary be defined by two
concentric circles
$\Sigma=\\{x\in\mathbb{R}^{2}\mid\left|{x}\right|=r_{\Sigma}\\}$ and
$\Gamma=\\{x\in\mathbb{R}^{2}\mid\left|{x}\right|=\rho^{0}\\}$. We define
$\Gamma^{0}(t)=\\{x\in\mathbb{R}^{2}\mid\left|{x}\right|=\rho^{0}(t)\\}$, and
consider, in polar coordinate $(r,\theta)$, its perturbation given by

$\Gamma^{\varepsilon}(t)=\\{(r,\theta)\mid r=\rho^{\varepsilon}(\theta,t)\\},$

where we suppose that $\rho^{\varepsilon}$ formally has the expansion

$\rho^{\varepsilon}(\theta,t)=\rho^{0}(t)+\varepsilon\rho^{1}(\theta,t)+O(\varepsilon^{2}).$

$\Sigma$$\Omega^{\varepsilon}(t)$$\Gamma^{\varepsilon}(t)$$\nu^{\varepsilon}$$\nu$
Figure 3: The axisymmetric annular domain

Here, in particular, we shall consider perturbation of the type
$\rho^{1}(\theta,t)=R(t)\cos k\theta$ or $R(t)\sin k\theta$. Our aim is to
derive a necessary condition of $R(t)$ for $\Gamma^{\varepsilon}(t)$
satisfying our main system (1) with suitable $u_{\text{D}}^{\varepsilon}(x,t)$
and $u_{\text{N}}^{\varepsilon}(x,t)$, except for the initial condition. For
simplicity, let us suppose that $f$ and $g$ are positive constants. We define

$\Omega^{\varepsilon}(t)=\\{x\in\mathbb{R}^{2}\mid\rho^{\varepsilon}(t)<\left|{x}\right|<r_{\Sigma}\\}.$

Of course, we set
$\Omega^{0}(t)=\\{x\in\mathbb{R}^{2}\mid\rho^{0}(t)<\left|{x}\right|<r_{\Sigma}\\}$.

We consider the following systems of PDEs:

$\displaystyle\left\\{\begin{aligned} \Delta
u_{\text{D}}^{\varepsilon}=0\quad\text{in}\ \Omega^{\varepsilon}(t),\qquad
u_{\text{D}}^{\varepsilon}=f>0\quad\text{on}\ \Sigma,\qquad
u_{\text{D}}^{\varepsilon}=0\quad\text{on}\ \Gamma^{\varepsilon}(t),\\\ \Delta
u_{\text{N}}^{\varepsilon}=0\quad\text{in}\
\Omega^{\varepsilon}(t),\qquad\dfrac{\partial{u_{\text{N}}^{\varepsilon}}}{\partial\nu}=g>0\quad\text{on}\
\Sigma,\qquad u_{\text{N}}^{\varepsilon}=0\quad\text{on}\
\Gamma^{\varepsilon}(t),\\\ \end{aligned}\right.$ (B.1)

where $\nu$ is the outward unit normal vector to $\Sigma$. Accordingly, the
normal velocity of the moving boundary $\Gamma^{\varepsilon}$ is given by

$v^{\varepsilon}(\theta,t)=-\left(\dfrac{\partial{u_{\text{D}}^{\varepsilon}}}{\partial\nu^{\varepsilon}}-\dfrac{\partial{u_{\text{N}}^{\varepsilon}}}{\partial\nu^{\varepsilon}}\right),\qquad
y\in\Gamma^{\varepsilon}(t),$

where $\nu^{\varepsilon}$ is the inward unit normal to
$\Gamma^{\varepsilon}(t)$ (see Figure 3 for illustration) given by

$\nu^{\varepsilon}(\theta,t)=\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}+\varepsilon\frac{\rho^{1}_{\theta}(\theta,t)}{\rho^{0}(t)}\begin{pmatrix}\sin\theta\\\
-\cos\theta\end{pmatrix}+O(\varepsilon^{2}).$

Here, the initial geometric profile of the perturbed boundary is
$\Gamma^{\varepsilon}(0)=\\{(r,\theta)\mid
r=\rho^{\varepsilon}_{0}(\theta)\\}$ where
$\rho^{\varepsilon}_{0}(\theta)=\rho^{0}_{0}+\varepsilon\rho^{1}_{0}(\theta)+O(\varepsilon^{2})$.
Now, we assume that

$u_{\text{p}}^{\varepsilon}(x,t)=u_{\text{p}}^{0}(x,t)+\varepsilon
u_{\text{p}}^{1}(x,t)+O(\varepsilon^{2}),\qquad\text{p}=\text{D},\text{N},$

in a neighborhood of $\Omega^{\varepsilon}(t)$. So, for
$x=\rho^{0}(t)\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\in\Gamma^{0}(t)$ and
$y=\rho^{\varepsilon}(\theta,t)\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\in\Gamma^{\varepsilon}(t)$ we have, after applying
Taylor’s expansion and combining terms with respect to $\varepsilon$, the
following identities

$\displaystyle u_{\text{p}}^{\varepsilon}(y,t)$
$\displaystyle=u_{\text{p}}^{0}(x,t)+\varepsilon\left(u_{\text{p}}^{1}(x,t)+\rho^{1}(\theta,t)\frac{\partial
u_{\text{p}}^{0}}{\partial{r}}(x,t)\right)+O(\varepsilon^{2}),$
$\displaystyle=u_{\text{p}}^{0}(x,t)+\varepsilon\left(u_{\text{p}}^{1}(x,t)+\rho^{1}(\theta,t)\frac{\partial
u_{\text{p}}^{0}}{\partial{r}}(x,t)\right)+O(\varepsilon^{2}),$
$\displaystyle\frac{\partial
u_{\text{p}}^{\varepsilon}}{\partial\nu^{\varepsilon}}(y,t)$
$\displaystyle=\frac{\partial
u_{\text{p}}^{0}}{\partial{r}}(x,t)+\varepsilon\left(\frac{\partial
u_{\text{p}}^{1}}{\partial{r}}(x,t)+\rho^{1}(\theta,t)\frac{\partial^{2}u_{\text{p}}^{0}}{\partial{r}^{2}}(x,t)\right)+O(\varepsilon^{2}).$

The equations and terms of order $O(\varepsilon^{0})$ are as follows:

$\displaystyle\left\\{\begin{aligned} \Delta u_{\text{D}}^{0}=0\quad\text{in}\
\Omega^{0}(0),\qquad u_{\text{D}}^{0}=f\quad\text{on}\ \Sigma,\qquad
u_{\text{D}}^{0}=0\quad\text{on}\ \Gamma^{0}(0),\\\ \Delta
u_{\text{N}}^{0}=0\quad\text{in}\
\Omega^{0}(0),\qquad\dfrac{\partial{u_{\text{N}}^{0}}}{\partial\nu}=g\quad\text{on}\
\Sigma,\qquad u_{\text{N}}^{0}=0\quad\text{on}\ \Gamma^{0}(0),\\\
\end{aligned}\right.$ (B.2)

with

$\rho_{t}^{0}(t)=-\left(\frac{\partial
u_{\text{D}}^{0}}{\partial{r}}(x,t)-\frac{\partial
u_{\text{N}}^{0}}{\partial{r}}(x,t)\right),\quad
x\in\Gamma^{0}(t),\qquad\text{and}\qquad\rho^{0}(0)=\rho_{0}^{0}.$ (B.3)

Meanwhile, the equations and terms of order $O(\varepsilon^{1})$ are as
follows:

$\displaystyle\left\\{\begin{aligned} \Delta
u_{\text{D}}^{1}&=0\quad\text{in}\ \Omega^{0}(t),\\\
u_{\text{D}}^{1}&=f\quad\text{on}\ \Sigma,\\\
u_{\text{D}}^{1}(x,t)+\rho^{1}(\theta,t)\frac{\partial
u_{\text{D}}^{0}}{\partial{r}}(x,t)&=0\quad\text{on}\ \Gamma^{0}(t),\\\ \Delta
u_{\text{N}}^{1}&=0\quad\text{in}\ \Omega^{0}(t),\\\
\dfrac{\partial{u_{\text{N}}^{1}}}{\partial\nu}&=g\quad\text{on}\ \Sigma,\\\
u_{\text{N}}^{1}(x,t)+\rho^{1}(\theta,t)\frac{\partial
u_{\text{N}}^{0}}{\partial{r}}(x,t)&=0\quad\text{on}\ \Gamma^{0}(t),\\\
\end{aligned}\right.$ (B.4)

with

$\rho_{t}^{1}(\theta,t)=-\left[\left(\frac{\partial
u_{\text{D}}^{1}}{\partial{r}}(x,t)-\frac{\partial
u_{\text{N}}^{1}}{\partial{r}}(x,t)\right)+\left(\frac{\partial^{2}u_{\text{D}}^{0}}{\partial{r}^{2}}(x,t)-\frac{\partial^{2}u_{\text{N}}^{0}}{\partial{r}^{2}}(x,t)\right)\rho^{1}(\theta,t)\right]$
(B.5)

and

$\rho^{1}(\theta,0)=\rho_{0}^{1}(\theta).$

Next, we compute the exact solutions to (B.2). We let

$u_{\text{p}}^{0}=C_{\text{p}}\log\dfrac{\left|{x}\right|}{\rho^{0}(t)},\quad\text{for
\text{p} = \text{D}, \text{N}, \quad where $\left|{x}\right|=r$}.$

Then, for $f(t)=C_{\text{D}}\log\dfrac{r_{\Sigma}}{\rho^{0}(t)}$, we get
$C_{\text{D}}(t)=\dfrac{f(t)}{\log(r_{\Sigma}/\rho^{0}(t))}$. Meanwhile, for
$\dfrac{\partial
u_{\text{p}}^{0}}{\partial\nu}\Big{|}_{\Sigma}=\dfrac{C_{\text{p}}}{r_{\Sigma}}=g(t)$
at $\text{p}=\text{N}$ (note here that $\nu$ is the outward unit normal to
$\Sigma$), we get $C_{\text{N}}(t)=g(t)r_{\Sigma}$. Therefore, we obtain the
following exact solutions to (B.2):

$\displaystyle u_{\text{D}}^{0}$
$\displaystyle=C_{\text{D}}(t)\log\left(\frac{\left|{x}\right|}{\rho^{0}(t)}\right)=\frac{f(t)}{\log\left(\frac{r_{\Sigma}}{\rho^{0}(t)}\right)}\log\left(\frac{\left|{x}\right|}{\rho^{0}(t)}\right),$
$\displaystyle u_{\text{N}}^{0}$
$\displaystyle=C_{\text{N}}(t)\log\left(\frac{\left|{x}\right|}{\rho^{0}(t)}\right)=g(t)r_{\Sigma}\log\left(\frac{\left|{x}\right|}{\rho^{0}(t)}\right),$

which gives us

$w^{0}:=u_{\text{N}}^{0}-u_{\text{D}}^{0}=(C_{\text{N}}(t)-C_{\text{D}}(t))\log\left(\frac{\left|{x}\right|}{\rho^{0}(t)}\right)=\left(g(t)r_{\Sigma}-\frac{f(t)}{\log\left(\frac{r_{\Sigma}}{\rho^{0}(t)}\right)}\right)\log\left(\frac{\left|{x}\right|}{\rho^{0}(t)}\right).$

Using the above notation, and from (B.3) and (B.5), we get

$\displaystyle O(\varepsilon^{0}):\ \left\\{\begin{aligned}
\rho_{t}^{0}(t)&=\frac{\partial w^{0}}{\partial{r}}(x,t),\\\ \frac{\partial
w^{0}}{\partial{r}}(x,t)\Big{|}_{x\in\Gamma^{0}(t)}&=\frac{C_{\text{N}}(t)-C_{\text{D}}(t)}{\rho^{0}(t)},\qquad
x\in\Gamma^{0}(t),\end{aligned}\right.$ (B.6) $\displaystyle
O(\varepsilon^{1}):\ \left\\{\begin{aligned}
\rho_{t}^{1}(\theta,t)&=\frac{\partial
w^{1}}{\partial{r}}(x,t)+\frac{\partial^{2}w^{1}}{\partial{r}^{2}}(x,t)\rho^{1}(\theta,t),\\\
\frac{\partial^{2}w^{0}}{\partial{r}^{2}}(x,t)\Big{|}_{x\in\Gamma^{0}(t)}&=-\frac{C_{\text{N}}(t)-C_{\text{D}}(t)}{\rho^{0}(t)^{2}},\quad
x=\rho^{0}(t)\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\in\Gamma^{0}(t),\end{aligned}\right.$ (B.7)

where

$C_{\text{D}}(t)=\frac{f(t)}{\log\left(\frac{r_{\Sigma}}{\rho^{0}(t)}\right)}>0\qquad\text{and}\qquad
C_{\text{N}}(t)=g(t)r_{\Sigma}>0.$ (B.8)

Let us now suppose, as mentioned earlier, that $\rho^{1}(\theta,t)=R(t)\cos
k\theta$, $k\in\mathbb{N}$. Then, we consider the ansatz

$u_{\text{p}}^{1}=a_{\text{p}}(t)r^{k}\cos k\theta+b_{\text{p}}(t)r^{-k}\cos
k\theta,\qquad r=\left|{x}\right|,\qquad\text{p}=\text{D},\text{N},$

which leads to the equivalence

$u_{\text{p}}^{1}+\rho^{1}\frac{\partial
u_{\text{p}}^{0}}{\partial{r}}=0\qquad\Longleftrightarrow\qquad
a_{\text{p}}(t)\rho^{0}(t)^{k}+b_{\text{p}}(t)\rho^{0}(t)^{-k}+R(t)\frac{C_{\text{p}}(t)}{\rho^{0}(t)}=0.$
(B.9)

So, for the first-order term, we have

$\displaystyle w^{1}$
$\displaystyle:=u_{\text{N}}^{1}-u_{\text{D}}^{1}=\left\\{(a_{\text{N}}(t)-a_{\text{D}}(t))r^{k}+(b_{\text{N}}(t)-b_{\text{D}}(t))r^{-k}\right\\}\cos
k\theta,$ $\displaystyle\frac{\partial
w^{1}}{\partial{r}}\Big{|}_{\Gamma^{0}(t)}$
$\displaystyle=k\left\\{(a_{\text{N}}(t)-a_{\text{D}}(t))\rho^{0}(t)^{k-1}-(b_{\text{N}}(t)-b_{\text{D}}(t))\rho^{0}(t)^{-(k+1)}\right\\}\cos
k\theta.$

From (B.7)1 and the fact that $\rho^{1}(\theta,t)=R(t)\cos k\theta$, we get

$R^{\prime}(t)=\lambda_{k}(t)R(t),$

where

$\lambda_{k}(t)=k\left\\{(\tilde{a}_{\text{N}}(t)-\tilde{a}_{\text{D}}(t))\rho^{0}(t)^{k-1}-(\tilde{b}_{\text{N}}(t)-\tilde{b}_{\text{D}}(t))\rho^{0}(t)^{-(k+1)}\right\\}-\frac{C_{\text{N}}(t)-C_{\text{D}}(t)}{\rho^{0}(t)^{2}},$

where $\tilde{a}_{\text{p}}=a_{\text{p}}/R(t)$ and
$\tilde{b}_{\text{p}}=b_{\text{p}}/R(t)$.

To get the exact form of $\lambda_{k}(t)$, we need to compute for the
coefficients $a_{\text{D}}(t)$, $b_{\text{D}}(t)$, $a_{\text{N}}(t)$, and
$b_{\text{N}}(t)$. To this end, we note the following implications of the
boundary conditions

$\displaystyle u_{\text{D}}^{1}=0\quad\text{on $\Sigma$}$
$\displaystyle\Longleftrightarrow\qquad
a_{\text{D}}(t)r_{\Sigma}^{k}+b_{\text{D}}(t)r_{\Sigma}^{-k}=0,$ (B.10)
$\displaystyle\frac{\partial u_{\text{N}}^{1}}{\partial\nu}=0\quad\text{on
$\Sigma$}$ $\displaystyle\Longleftrightarrow\qquad
a_{\text{D}}(t)r_{\Sigma}^{k-1}-b_{\text{D}}(t)r_{\Sigma}^{-(k+1)}=0.$

Equations (B.9) and (B.10) respectively lead to a system of equations which,
upon solving, provides the forms of $a_{\text{D}}$ and $b_{\text{D}}$ as well
as $a_{\text{N}}$ and $b_{\text{N}}$. That is, we have

$\begin{pmatrix}\rho^{0}(t)^{k}&\rho^{0}(t)^{-k}\\\
r_{\Sigma}^{k}&r_{\Sigma}^{-k}\end{pmatrix}\begin{pmatrix}\tilde{a}_{\text{D}}(t)\\\
\tilde{b}_{\text{D}}(t)\end{pmatrix}=\begin{pmatrix}-\frac{C_{\text{D}}(t)}{\rho^{0}(t)}\\\
0\end{pmatrix}$

and

$\begin{pmatrix}\rho^{0}(t)^{k}&\rho^{0}(t)^{-k}\\\
r_{\Sigma}^{k-1}&-r_{\Sigma}^{-(k+1)}\end{pmatrix}\begin{pmatrix}\tilde{a}_{\text{N}}(t)\\\
\tilde{b}_{\text{N}}(t)\end{pmatrix}=\begin{pmatrix}-\frac{C_{\text{N}}(t)}{\rho^{0}(t)}\\\
0\end{pmatrix},$

which would give us

$\displaystyle\begin{pmatrix}\tilde{a}_{\text{D}}(t)\\\
\tilde{b}_{\text{D}}(t)\end{pmatrix}$
$\displaystyle=\frac{1}{\rho^{0}(t)^{k}r_{\Sigma}^{-k}-\rho^{0}(t)^{-k}r_{\Sigma}^{k}}\begin{pmatrix}r_{\Sigma}^{-k}&-\rho^{0}(t)^{-k}\\\
-r_{\Sigma}^{k}&\rho^{0}(t)^{k}\end{pmatrix}\begin{pmatrix}-\frac{C_{\text{D}}(t)}{\rho^{0}(t)}\\\
0\end{pmatrix}$
$\displaystyle=\frac{1}{\operatorname{{det}}^{k}_{\text{D}}(t)}\left(\frac{C_{\text{D}}(t)}{\rho^{0}(t)}\right)\begin{pmatrix}r_{\Sigma}^{-k}\\\\[5.0pt]
-r_{\Sigma}^{k}\\\ \end{pmatrix},$
$\displaystyle\begin{pmatrix}\tilde{a}_{\text{N}}(t)\\\
\tilde{b}_{\text{N}}(t)\end{pmatrix}$
$\displaystyle=\frac{1}{-\rho^{0}(t)^{k}r_{\Sigma}^{-(k+1)}-\rho^{0}(t)^{-k}r_{\Sigma}^{k-1}}\begin{pmatrix}-r_{\Sigma}^{-(k+1)}&-\rho^{0}(t)^{-k}\\\
-r_{\Sigma}^{k-1}&\rho^{0}(t)^{k}\end{pmatrix}\begin{pmatrix}-\frac{C_{\text{N}}(t)}{\rho^{0}(t)}\\\
0\end{pmatrix}$
$\displaystyle=\frac{1}{\operatorname{{det}}^{k}_{\text{N}}(t)}\left(\frac{g(t)}{\rho^{0}(t)}\right)\begin{pmatrix}r_{\Sigma}^{-k}\\\\[5.0pt]
r_{\Sigma}^{k}\\\ \end{pmatrix},$

where

$\displaystyle\operatorname{{det}}^{k}_{\text{D}}(t)$
$\displaystyle=-\rho^{0}(t)^{k}r_{\Sigma}^{-k}+\rho^{0}(t)^{-k}r_{\Sigma}^{k}=\left(\frac{r_{\Sigma}}{\rho^{0}(t)}\right)^{k}-\left(\frac{\rho^{0}(t)}{r_{\Sigma}}\right)^{k},$
$\displaystyle\operatorname{{det}}^{k}_{\text{N}}(t)$
$\displaystyle=\rho^{0}(t)^{k}r_{\Sigma}^{-k}+\rho^{0}(t)^{-k}r_{\Sigma}^{k}=\left(\frac{\rho^{0}(t)}{r_{\Sigma}}\right)^{k}+\left(\frac{r_{\Sigma}}{\rho^{0}(t)}\right)^{k}.$

Note that both of these determinants are positive because
$\rho^{0}(t)<r_{\Sigma}$ for all $t>0$. Moreover, we observe that

$0<\frac{\operatorname{{det}}^{k}_{\text{N}}(t)}{\operatorname{{det}}^{k}_{\text{D}}(t)}\longrightarrow
1,\qquad
0<\frac{\operatorname{{det}}^{k}_{\text{D}}(t)}{\operatorname{{det}}^{k}_{\text{N}}(t)}\longrightarrow
1,\qquad\text{as $k\to\infty$}.$

Then, finally, we have

$\displaystyle\lambda_{k}(t)$
$\displaystyle=\frac{k}{\rho^{0}(t)}\left\\{\left[-\tilde{a}_{\text{D}}(t)\rho^{0}(t)^{k}+\tilde{b}_{\text{D}}(t)\rho^{0}(t)^{-k}\right]+\left[\tilde{a}_{\text{N}}(t)\rho^{0}(t)^{k}-\tilde{b}_{\text{N}}(t)\rho^{0}(t)^{-k}\right]\right\\}$
$\displaystyle\qquad+\frac{C_{\text{D}}(t)-C_{\text{N}}(t)}{\rho^{0}(t)^{2}}$
$\displaystyle=\frac{k}{\rho^{0}(t)^{2}}\left\\{\frac{-C_{\text{D}}(t)}{\operatorname{{det}}^{k}_{\text{D}}(t)}\left(r_{\Sigma}^{-k}\rho^{0}(t)^{k}+r_{\Sigma}^{k}\rho^{0}(t)^{-k}\right)+\frac{g(t)}{\operatorname{{det}}^{k}_{\text{N}}(t)}\left(r_{\Sigma}^{-k}\rho^{0}(t)^{k}-r_{\Sigma}^{k}\rho^{0}(t)^{-k}\right)\right\\}$
$\displaystyle\qquad+\frac{C_{\text{D}}(t)-C_{\text{N}}(t)}{\rho^{0}(t)^{2}}$
$\displaystyle=-\underbrace{\frac{k}{\rho^{0}(t)^{2}}\left(C_{\text{D}}(t)\frac{\operatorname{{det}}^{k}_{\text{N}}(t)}{\operatorname{{det}}^{k}_{\text{D}}(t)}+g(t)\frac{\operatorname{{det}}^{k}_{\text{D}}(t)}{\operatorname{{det}}^{k}_{\text{N}}(t)}\right)}_{>0}+\frac{C_{\text{D}}(t)-C_{\text{N}}(t)}{\rho^{0}(t)^{2}}.$

Evidently, for large enough $k$, $\lambda_{k}(t)$ is negative, for any
$t\geqslant 0$. This implies, formally, that the main system (1) is well-posed
in the case of annular domains formed by concentric radially symmetric shapes.

##  Appendix C Change of variables

For the benefit of the reader, we provide here the computation of the
transformed domain. We focus on rewriting the main equations for the variable
$u_{\text{N}}$ (the fourth to sixth equations in (1)). For the case
$u_{\text{D}}$, the computations are similar. So, consider a smooth function
$u_{\text{N}}(x,t)$, $x\in\Omega(t)$, and $U_{\text{N}}(y,t)$, $y\in\Omega$,
where $y=Y(x,t)$. Then, $x=Y^{-1}(y,t)$), and we have

$\displaystyle
u_{\text{N}}(x,t)=U_{\text{N}}(y,t)=u_{\text{N}}(Y^{-1}(y,t),t)=:u_{\text{N}}\circ
Y^{-1}(y,t).$ (C.1)

Let us introduce the notation
$x:=(x_{1},\ldots,x_{d})^{\top}\in\mathbb{R}^{d}$,
$y:=(y_{1},\ldots,y_{d})^{\top}\in\mathbb{R}^{d}$, so that
$Y^{-1}(y,t):=(x_{1},\ldots,x_{d})^{\top}(y,t)\in\mathbb{R}^{d}$. For clarity,
we shortly verify here identity (29).

By generalized chain rule, the $m$th ($m=1,\ldots,n$) column entry of
$\nabla_{y}^{\top}$ is computed as follows (dropping $t$)

$\displaystyle\frac{\partial}{\partial
x_{m}}u_{\text{N}}(x)=\sum_{k=1}^{d}\frac{\partial U_{\text{N}}}{\partial
y_{k}}(y)\frac{\partial Y_{k}}{\partial x_{m}}(x)=\sum_{k=1}^{d}\frac{\partial
U_{\text{N}}}{\partial
y_{k}}(y)\textsf{J}^{km}=\sum_{k=1}^{d}\textsf{J}^{km}\frac{\partial
U_{\text{N}}}{\partial y_{k}}(y).$

We next focus on the transformation of the Laplace equation. We recall that

$\nabla_{x}\cdot\nabla_{x}u=\Delta_{x}u=\left(\frac{\partial^{2}}{\partial
x_{1}^{2}}+\cdots+\frac{\partial^{2}}{\partial x_{d}^{2}}\right)u.$

From the relation $u_{\text{N}}(x,t)=U_{\text{N}}(y,t)$, we get
$\Delta_{x}u_{\text{N}}(x,t)=\Delta_{y}U_{\text{N}}(y,t)=\nabla_{y}^{\top}\cdot\nabla_{y}^{\top}U_{\text{N}}(y,t)$.
Each summand $\partial^{2}u/\partial x_{k}^{2}$, $k=1,\ldots,n$, is computed
by performing the following calculation of partial derivatives

$\displaystyle\frac{\partial}{\partial x_{k}}\left(\frac{\partial}{\partial
x_{k}}u_{\text{N}}(x)\right)$ $\displaystyle=\frac{\partial}{\partial
x_{k}}\left(\sum_{p=1}^{d}\frac{\partial U_{\text{N}}}{\partial
y_{p}}(y)\frac{\partial Y_{p}}{\partial x_{k}}(x)\right)$
$\displaystyle=\sum_{p=1}^{d}\left\\{\underbrace{\frac{\partial}{\partial
x_{k}}\left(\frac{\partial U_{\text{N}}}{\partial
y_{p}}(y)\right)}_{=:D_{1}}\frac{\partial Y_{p}}{\partial
x_{k}}(x)+\frac{\partial U_{\text{N}}}{\partial
y_{p}}(y)\underbrace{\frac{\partial}{\partial x_{k}}\left(\frac{\partial
Y_{p}}{\partial x_{k}}(x)\right)}_{=:D_{2}}\right\\}.$

The partial derivative $D_{1}$ is computed as follows:

$D_{1}=\frac{\partial}{\partial x_{k}}\left(\frac{\partial
U_{\text{N}}}{\partial y_{p}}(y)\right)=\sum_{m=1}^{d}\frac{\partial
U_{\text{N}}}{\partial y_{m}y_{p}}(y)\frac{\partial Y_{m}}{\partial
x_{k}}(x)=\sum_{m=1}^{d}\frac{\partial U_{\text{N}}}{\partial
y_{m}y_{p}}(y)\textsf{J}^{mk}.$

Hence, the first summand where $D_{1}$ appears can be written as follows:

$\displaystyle\sum_{p=1}^{d}\frac{\partial}{\partial
x_{k}}\left(\frac{\partial U_{\text{N}}}{\partial
y_{p}}(y)\right)\frac{\partial Y_{p}}{\partial x_{k}}(x)$
$\displaystyle=\sum_{p=1}^{d}\left(\sum_{m=1}^{d}\frac{\partial
U_{\text{N}}}{\partial y_{m}y_{p}}(y)\textsf{J}^{mk}\right)\textsf{J}^{pk}.$

Now, summing $k$ from $1$ to $d$, and after a few rearrangements, we get

$\displaystyle\sum_{k=1}^{d}\sum_{p=1}^{d}\sum_{m=1}^{d}\textsf{J}^{mk}\textsf{J}^{pk}\frac{\partial^{2}U_{\text{N}}}{\partial
y_{m}\partial y_{p}}$
$\displaystyle=\sum_{m,p=1}^{d}\left(\sum_{k=1}^{d}\textsf{J}^{mk}\textsf{J}^{pk}\right)\frac{\partial^{2}U_{\text{N}}}{\partial
y_{m}\partial y_{p}}$
$\displaystyle=\sum_{m,p=1}^{d}{{\textsf{A}}}_{mp}\frac{\partial^{2}U_{\text{N}}}{\partial
y_{m}\partial y_{p}},$

where ${{\textsf{A}}}_{mp}=\sum_{k=1}^{d}\textsf{J}^{mk}\textsf{J}^{pk}$ are
the entries of the matrix $A=\textsf{J}^{-1}\textsf{J}^{-\top}$. For the
partial derivative $D_{2}$, we have the form

$\sum_{k,m,p=1}^{d}\textsf{J}^{mp}\frac{\partial\textsf{J}^{pk}}{\partial
y_{m}}\frac{\partial U_{\text{N}}}{\partial
y_{p}}=\sum_{k,m,p=1}^{d}\frac{\partial Y_{m}}{\partial
x_{k}}(x)\frac{\partial}{\partial y_{m}}\left(\frac{\partial Y_{p}}{\partial
x_{k}}(x)\right)\frac{\partial U_{\text{N}}}{\partial y_{p}}(y).$

The following computations yield the above expression:

$\displaystyle D_{2}$ $\displaystyle=\frac{\partial}{\partial
x_{k}}\left(\frac{\partial Y_{p}}{\partial x_{k}}(x)\right)$
$\displaystyle=\frac{\partial}{\partial x_{k}}\left(\varphi\circ
Y(x)\right)\quad\left(\varphi(y):=\left(\frac{\partial Y^{-1}_{p}}{\partial
y_{k}}(y)\right)^{-1},\ \Omega\ni y=Y(x,t),\ x\in\Omega(t)\right)$
$\displaystyle=\sum_{m=1}^{d}\frac{\partial\varphi}{\partial
y_{m}}(Y(x))\frac{\partial Y_{m}}{\partial x_{k}}(x)$
$\displaystyle=\sum_{m=1}^{d}\textsf{J}^{mk}\frac{\partial\varphi}{\partial
y_{m}}(y)\quad(\text{$m$-$k$: row-column entry})$
$\displaystyle=\sum_{m=1}^{d}\textsf{J}^{mk}\frac{\partial}{\partial
y_{m}}\textsf{J}^{pk}.$

Inserting the above expression for $D_{2}$ in the second summand earlier
above, and then summing $k$ from $1$ to $d$, we get

$\displaystyle\sum_{k=1}^{d}\sum_{p=1}^{d}\frac{\partial
U_{\text{N}}}{\partial y_{p}}(y)\frac{\partial}{\partial
x_{k}}\left(\frac{\partial Y_{p}}{\partial x_{k}}(x)\right)$
$\displaystyle=\sum_{k=1}^{d}\sum_{m=1}^{d}\sum_{p=1}^{d}\frac{\partial
Y_{m}}{\partial x_{k}}(x)\frac{\partial\varphi}{\partial
y_{m}}(y)\frac{\partial U_{\text{N}}}{\partial y_{p}}(y)$
$\displaystyle=\sum_{k,m,p=1}^{d}\textsf{J}^{mk}\frac{\partial\textsf{J}^{pk}}{\partial
y_{m}}\frac{\partial U_{\text{N}}}{\partial y_{p}},$

where $\varphi(y)$ is again defined as before. In summary, the Laplace
equation $\Delta_{x}u_{\text{N}}(x,t)=0$, $x\in\Omega(t)$, $t>0$, when
expressed over the fixed domain $\Omega$, has the following form:

$\sum_{m,p=1}^{d}{{\textsf{A}}}_{mp}\frac{\partial^{2}U_{\text{N}}}{\partial
y_{m}\partial
y_{p}}+\sum_{k,m,p=1}^{d}\textsf{J}^{mk}\frac{\partial\textsf{J}^{pk}}{\partial
y_{m}}\frac{\partial U_{\text{N}}}{\partial y_{p}}=0,$

where ${{\textsf{A}}}_{mp}$ is as mentioned earlier. For the Laplace equation
$\Delta_{x}u_{\text{D}}(x,t)=0$, $x\in\Omega(t)$, $t>0$, simply replace
$U_{\text{N}}$ by $U_{\text{D}}$ above.

We next perform the transformations of equations on the boundary. First, it is
easy to check that $U_{\text{N}}(y,t)=u_{\text{N}}(Y^{-1}(y,t),t)=0$,
$y\in\Gamma$, $t>0$. Similarly, on $\Sigma$, we can write $g(x,t)$
($x\in\Sigma$) as $g(Y^{-1}(y,t),t)$. Since ${\tilde{\rho}}(x,t)\equiv 0$ on
$\Sigma$, we simply write the function $g(Y^{-1}(y,t),t)$ as $g(y,t)$. Now,
let us rewrite the expression $\nabla_{x}u_{\text{N}}(x,t)\cdot\nu(x,t)$,
where $x\in\partial\Omega(t)$, in terms of $U_{\text{N}}(y,t)$ on the fixed
boundary $\partial\Omega$. This is done as follows (hereon, we occasionally
drop $t$ and write $\nu_{\circ}:=\nu_{\circ}(y)$ for convenience):

$\displaystyle\nabla_{x}u_{\text{N}}(x)\cdot\nu(x)$
$\displaystyle=\textsf{J}^{-\top}\nabla_{y}^{\top}U_{\text{N}}(y)\cdot\frac{\textsf{J}^{-\top}\nu_{\circ}}{|\textsf{J}^{-\top}\nu_{\circ}|}(y).$

Clearly, on the fixed (interior) boundary $\Sigma$, we have

$\dfrac{\partial{}}{\partial{\nu_{\circ}}}U_{\text{N}}(y,t)=g(y,t),\qquad
y\in\Sigma,\ t>0.$

In the same manner, we obtain the following transformation of boundary
conditions for the variable $u_{\text{D}}$

$U_{\text{D}}(y,t)=u_{\text{D}}(Y^{-1}(y,t),t)=0,\qquad y\in\Gamma,\ t>0,$

and, on $\Sigma$, we can write $f(x,t)$ ($x\in\Sigma$) as $f(Y^{-1}(y,t),t)$.
Again, because ${\tilde{\rho}}(x,t)\equiv 0$ on $\Sigma$, we simply write the
function $f(Y^{-1}(y,t),t)$ as $f(y,t)$. Moreover, from $\Gamma(t)$ to
$\Gamma$, we have the following sequence of computations and transformations
with regards to the equation for the normal velocity

$\displaystyle V_{n}(x,t)$
$\displaystyle=\boldsymbol{\mathbf{V}}(x,t)\cdot\nu(x,t)=\frac{d}{dt}x(t)\cdot\nu(x,t),$
$\displaystyle=\frac{d}{dt}\left(y+\boldsymbol{\mathbf{N}}(y){\tilde{\rho}}(y,t)\right)\cdot\nu(x,t)$
$\displaystyle=\dfrac{\partial{}}{\partial{t}}{\tilde{\rho}}(y,t)\boldsymbol{\mathbf{N}}(y)\cdot\frac{\textsf{J}^{-\top}\nu_{\circ}}{|\textsf{J}^{-\top}\nu_{\circ}|}(y),\quad\left({\tilde{\rho}_{t}}(y,t):=\frac{d}{dt}{\tilde{\rho}}(y,t),\
y\in\Gamma,\ t>0\right).$

and

$\displaystyle\left(\nabla u_{\text{D}}(x,t)-\nabla
u_{\text{N}}(x,t)\right)\cdot\nu(x,t)=\textsf{J}^{-\top}\left(\nabla_{y}^{\top}U_{\text{D}}(y)-\nabla_{y}^{\top}U_{\text{N}}(y)\right)\cdot\frac{\textsf{J}^{-\top}\nu_{\circ}}{|\textsf{J}^{-\top}\nu_{\circ}|}(y),$

for $x=x(t)\in\Gamma(t)$, $y\in\Gamma$, and $t>0$. These identities, with
${\textsf{A}}:=\textsf{J}^{-1}\textsf{J}^{-\top}$, gives us

$\dfrac{\partial{}}{\partial{t}}{\tilde{\rho}}(y,t)=-\left(\boldsymbol{\mathbf{N}}(y)\cdot\textsf{J}^{-\top}\nu_{\circ}(y)\right)^{-1}\left(\nabla_{y}^{\top}U_{\text{D}}(y)-\nabla_{y}^{\top}U_{\text{N}}(y)\right)\cdot{{\textsf{A}}}\nu_{\circ},$

for $y\in\Gamma$, $t>0$. Taking into account the fact that
$U_{\text{N}}(y,t)=U_{\text{D}}(y,t)=0$ on $\Gamma$, $t>0$, we can
equivalently write the above equation as follows:

$\displaystyle\dfrac{\partial{}}{\partial{t}}{\tilde{\rho}}(y,t)+B_{{\tilde{\rho}}}\dfrac{\partial{}}{\partial{\nu_{\circ}}}\left(U_{\text{D}}(y,t)-U_{\text{N}}(y,t)\right)=0,\qquad(y\in\Gamma,\
t>0),$

where
$B_{\tilde{\rho}}:=B_{{\tilde{\rho}}}(y)=\left(\boldsymbol{\mathbf{N}}(y)\cdot\textsf{J}^{-\top}\nu_{\circ}(y)\right)^{-1}\left(\nu_{\circ}\cdot{{\textsf{A}}}\nu_{\circ}\right)$.
We remind that throughout the paper we shall assume that
${\tilde{\rho}}(y,t)$, $y\in\Gamma$, $t\in[0,T]$, is small enough so that
$\boldsymbol{\mathbf{N}}(y)\cdot\textsf{J}^{-\top}\nu_{\circ}(y)$ and
$|\textsf{J}^{-1}\nu_{\circ}(y)|$ are strictly positive for $y\in\Gamma$, and
that the map $Y^{-1}(y,t)$ is invertible.

In the rest of the section, we look at how the transformed Laplace equation
splits into the operators given in (32). The form of
$\pazocal{{L}}_{\tilde{\rho}}$ is already clear, so we focus on rewriting the
expression

$\sum_{k,m,p=1}^{d}\textsf{J}^{mk}\frac{\partial\textsf{J}^{pk}}{\partial
y_{m}}\frac{\partial U_{\text{N}}}{\partial y_{p}}(y),\qquad(y\in\Omega),$
(C.2)

in terms of $\boldsymbol{\mathbf{N}}(y)$ and ${\tilde{\rho}}(y)$ (again, we
drop the dependence of ${\tilde{\rho}}$ to $t$ for convenience) to confirm the
forms of $\pazocal{{K}}_{\tilde{\rho}}$ and $\pazocal{{M}}_{\tilde{\rho}}$.
From the identity $\textsf{J}\textsf{J}^{-1}=\textsf{I}$, we have
$\partial_{y}\textsf{J}\textsf{J}^{-1}+\textsf{J}\partial_{y}\textsf{J}^{-1}=0$,
or equivalently,
$\partial_{y}\textsf{J}^{-1}=-\textsf{J}^{-1}\left(\partial_{y}\textsf{J}\right)\textsf{J}^{-1}$.
We compute $\partial_{y}\textsf{J}$ by differentiating the equation
$\textsf{J}=\textsf{I}+\boldsymbol{\mathbf{N}}\nabla_{y}^{\top}{\tilde{\rho}}+(\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}){\tilde{\rho}}$
with respect to $y$:

$\partial_{y}\textsf{J}=\boldsymbol{\mathbf{N}}\operatorname{Hess}_{y}{\\!({\tilde{\rho}})}+\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}\left(\nabla_{y}^{\top}{\tilde{\rho}}\right)+\operatorname{Hess}_{y}{\\!(\boldsymbol{\mathbf{N}})}{\tilde{\rho}}+(\nabla_{y}^{\top}\boldsymbol{\mathbf{N}})\nabla_{y}^{\top}{\tilde{\rho}},$

where $\operatorname{Hess}_{y}{\\!(\cdot)}$ denotes the bilinear form
associated with the Hessian matrix, i.e.,
$\operatorname{Hess}_{y}{\\!({\tilde{\rho}})}:=\operatorname{Hess}_{y}{\\!({\tilde{\rho}}(y))}=\partial_{y}(\nabla^{\top}{\tilde{\rho}}(y))$,
$y\in\Omega$. Thus, we have

$\displaystyle\partial_{y}\textsf{J}^{-1}$
$\displaystyle=-\textsf{J}^{-1}\boldsymbol{\mathbf{N}}\operatorname{Hess}_{y}{\\!({\tilde{\rho}})}\textsf{J}^{-1}-\textsf{J}^{-1}\left[\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}\left(\nabla_{y}^{\top}{\tilde{\rho}}\right)+(\nabla_{y}^{\top}\boldsymbol{\mathbf{N}})\nabla_{y}^{\top}{\tilde{\rho}}+\operatorname{Hess}_{y}{\\!(\boldsymbol{\mathbf{N}})}{\tilde{\rho}}\right]\textsf{J}^{-1}.$

The entries of the above matrix can be computed without difficulty. For
instance, the matrix
$-\textsf{J}^{-1}\boldsymbol{\mathbf{N}}\operatorname{Hess}_{y}{\\!({\tilde{\rho}})}\textsf{J}^{-1}$
has entries (in Einstein’s notation)

$-\textsf{J}^{pj}\left(\boldsymbol{\mathbf{N}}\partial_{y_{m}}{(\nabla^{\top}{\tilde{\rho}})}\right)_{jq}\textsf{J}^{qk}\\\
=-\textsf{J}^{pj}\left(N_{j}\frac{\partial}{\partial
y_{m}}\left(\frac{\partial{\tilde{\rho}}}{\partial
y_{q}}\right)\right)\textsf{J}^{qk}.$

For the sake of clarity, we provide below a detailed computation of the matrix
$\partial_{y}\textsf{J}^{-1}$ in terms of its entries. To start, we apply the
matrix inversion lemma or the so-called Sherman–Morrison–Woodbury formula [58,
59, 60] to obtain the following identity

$\displaystyle D_{x}Y(x)$ $\displaystyle=\left[D_{y}Y^{-1}(y)\right]^{-1}$
$\displaystyle=\left[\textsf{I}+(\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}){\tilde{\rho}}+\boldsymbol{\mathbf{N}}\nabla_{y}^{\top}{\tilde{\rho}}\right]^{-1},\qquad(\boldsymbol{\mathbf{N}}:=\boldsymbol{\mathbf{N}}(y),\
{\tilde{\rho}}:={\tilde{\rho}}(y),\ y\in\Omega),$
$\displaystyle=\textsf{I}-\left[I+(\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}){\tilde{\rho}}+\boldsymbol{\mathbf{N}}\nabla_{y}^{\top}{\tilde{\rho}}\right]^{-1}\left((\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}){\tilde{\rho}}+\boldsymbol{\mathbf{N}}\nabla_{y}^{\top}{\tilde{\rho}}\right),$
$\displaystyle=\textsf{I}-\left[D_{y}Y^{-1}(y)\right]^{-1}\left((\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}){\tilde{\rho}}+\boldsymbol{\mathbf{N}}\nabla_{y}^{\top}{\tilde{\rho}}\right),$

where $\textsf{I}:=(\delta_{ij})\in\mathbb{R}^{d\times d}$ ($i,j=1,\ldots,d$)
denotes the $d\times d$ identity matrix. The above equations mean that

$\textsf{J}^{-1}=\textsf{I}-\textsf{J}^{-1}\left((\nabla_{y}^{\top}\boldsymbol{\mathbf{N}}){\tilde{\rho}}+\boldsymbol{\mathbf{N}}\nabla_{y}^{\top}{\tilde{\rho}}\right)=:\textsf{I}-\textsf{J}^{-1}\textsf{M},$

where, of course,

$\textsf{J}^{-1}=\left(\textsf{J}^{ij}\right)=(\textsf{J}_{ij})^{-1}=\left(\delta_{ij}+N_{i}\frac{\partial{\tilde{\rho}}}{\partial
y_{j}}+\frac{\partial N_{i}}{\partial y_{j}}{\tilde{\rho}}\right)^{-1}.$

Differentiating the above equation with respect to $y$, we get

$\displaystyle\partial_{y}\textsf{J}^{-1}=-\partial_{y}\textsf{J}^{-1}\textsf{M}-\textsf{J}^{-1}\partial_{y}\textsf{M}$
$\displaystyle\Longleftrightarrow$
$\displaystyle\partial_{y}\textsf{J}^{-1}(\textsf{I}+\textsf{M})$
$\displaystyle=-\textsf{J}^{-1}\partial_{y}\textsf{M}$
$\displaystyle\Longleftrightarrow$ $\displaystyle\partial_{y}\textsf{J}^{-1}$
$\displaystyle=-\textsf{J}^{-1}(\partial_{y}\textsf{M})(\textsf{I}+\textsf{M})^{-1}$
$\displaystyle\Longleftrightarrow$ $\displaystyle\partial_{y}\textsf{J}^{-1}$
$\displaystyle=-\textsf{J}^{-1}(\partial_{y}\textsf{M})\textsf{J}^{-1}.$

We recall that the product of three matrices, say
$\textsf{A}=(\textsf{A}_{ij})$, $\textsf{B}=(\textsf{B}_{ij})$, and
$\textsf{C}=(\textsf{C}_{ij})$, in terms of its $jm$-th entry, is given by
(using Einstein’s notation)
$(\textsf{A}\textsf{B}\textsf{C})_{ij}=\textsf{A}_{ik}\textsf{B}_{kl}\textsf{C}_{lj}$.
Therefore,
$\partial_{y_{i}}\textsf{J}^{-1}=-\textsf{J}^{-1}(\partial_{y_{i}}\textsf{M})\textsf{J}^{-1}$,
$i=1,\ldots,d$, in terms of its entries, is equivalent to

$\displaystyle\frac{\partial}{\partial y_{m}}\textsf{J}^{pk}$
$\displaystyle=-\sum_{j,q=1}^{d}\textsf{J}^{pj}(\partial_{y_{m}}\textsf{M})_{jq}\textsf{J}^{qk},\qquad(\text{$p$-$k$:
row-column entry}),$
$\displaystyle=-\sum_{j,q=1}^{d}\textsf{J}^{pj}\frac{\partial}{\partial
y_{m}}\left(\frac{\partial{N}_{j}}{\partial
y_{q}}{\tilde{\rho}}+{N}_{j}\frac{\partial{\tilde{\rho}}}{\partial
y_{q}}\right)\textsf{J}^{qk}$
$\displaystyle=-\sum_{j,q=1}^{d}\textsf{J}^{pj}\left(\frac{\partial^{2}{N}_{j}}{\partial
y_{m}\partial y_{q}}{\tilde{\rho}}+\frac{\partial{N}_{j}}{\partial
y_{q}}\frac{\partial{\tilde{\rho}}}{\partial
y_{m}}+\frac{\partial{N}_{j}}{\partial
y_{m}}\frac{\partial{\tilde{\rho}}}{\partial
y_{q}}+{N}_{j}\frac{\partial^{2}{\tilde{\rho}}}{\partial y_{m}\partial
y_{q}}\right)\textsf{J}^{qk}.$

Inserting the above expression to (C.2), we get the equivalent sum

$\displaystyle-\sum_{k,m,p=1}^{d}\textsf{J}^{mk}\left\\{\sum_{j,q=1}^{d}\textsf{J}^{pj}\left(\frac{\partial^{2}{N}_{j}}{\partial
y_{m}\partial y_{q}}{\tilde{\rho}}+\frac{\partial{N}_{j}}{\partial
y_{q}}\frac{\partial{\tilde{\rho}}}{\partial
y_{m}}+\frac{\partial{N}_{j}}{\partial
y_{m}}\frac{\partial{\tilde{\rho}}}{\partial
y_{q}}\right)\textsf{J}^{qk}\right\\}\frac{\partial U_{\text{N}}}{\partial
y_{p}}$
$\displaystyle\qquad-\sum_{k,m,p=1}^{d}\textsf{J}^{mk}\left(\sum_{j,q=1}^{d}\textsf{J}^{pj}{N}_{j}\frac{\partial^{2}{\tilde{\rho}}}{\partial
y_{m}\partial y_{q}}\textsf{J}^{qk}\right)\frac{\partial
U_{\text{N}}}{\partial y_{p}}\ (=:S_{1}+S_{2}).$

To get the desired expression (C.2) (with the same notation on indices), we
simply interchange the indices $k$ and $m$ in $S_{1}$, and for $S_{2}$, we
apply the change of notations on indices: $k\leftarrow j$, $j\leftarrow p$,
$p\leftarrow q$, and note that
${{\textsf{A}}}_{mp}=\sum_{k=1}^{d}\textsf{J}^{mk}\textsf{J}^{pk}$. The
resulting expression with these changes in notations finally provides the
desired expansion of (C.2) given in (32) from which the forms of the operators
$\pazocal{{K}}_{\tilde{\rho}}$ and $\pazocal{{M}}_{\tilde{\rho}}$ are made
clear.

##  Appendix D Computations of $\delta\textsf{A}^{(0)}$ and $\delta B_{0}$

In this appendix, we compute the variations $\delta\textsf{A}^{(0)}$ and
$\delta B_{0}$ for the more general variation

$\delta\pazocal{{F}}_{\rho_{0}}=\left.\frac{d}{d\lambda}\pazocal{{F}}_{\rho_{0}+\lambda({\tilde{\rho}}-\rho_{0})}\right|_{\lambda=0},$

where, basically, $\rho_{0}(y,0)=0$, $y\in\Sigma$, and $\rho_{0}(y,0)=0$,
$y\in\Gamma$ (see [14, Sec. 2]). In our case, $\rho_{0}\equiv 0$. For
convenience, let us first recall that the Jacobi matrix J has entries given by

$\textsf{J}_{km}=\delta_{km}+N_{k}\frac{\partial{\tilde{\rho}}}{\partial
y_{m}}+\frac{\partial N_{k}}{\partial y_{m}}{\tilde{\rho}},$

and that $\textsf{J}^{km}$ are, on the other hand, the entries of the inverse
matrix $\textsf{J}^{-1}$ which is the Jacobi matrix of the transform $Y(x,t)$,
and that $\textsf{J}^{-\top}=(\textsf{J}^{-1})^{\top}$. Here,
$\delta\textsf{A}^{(0)}=\textsf{J}_{0}^{-1}\textsf{J}_{0}^{-\top}$,
$\textsf{J}_{0}=\textsf{J}\big{|}_{{\tilde{\rho}}={\tilde{\rho}}_{0}}$, where
${\tilde{\rho}}_{0}$ satisfies the PDE system given by (27) with
${\tilde{\rho}}_{0}(y,0)=0$. Moreover, we recall the entries
${{\textsf{A}}}_{mp}:=\sum_{k=1}^{d}\textsf{J}^{mk}\textsf{J}^{pk}$ of the
matrix ${{\textsf{A}}}:=\textsf{J}^{-1}\textsf{J}^{-\top}$. We shall express
$\delta A$ in terms of $\delta\textsf{J}$, and the computation of
$\delta\textsf{A}^{(0)}$ goes as follows. First, we note that, from the
identity $\textsf{J}_{0}\textsf{J}_{0}^{-1}=\textsf{I}$, we get
$\delta(\textsf{J}_{0}\textsf{J}_{0}^{-1})=\textsf{J}_{0}(\delta\textsf{J}_{0}^{-1})+(\delta\textsf{J}_{0})\textsf{J}_{0}^{-1}=\boldsymbol{\mathbf{0}}$.
So,
$\delta\textsf{J}_{0}^{-1}=-\textsf{J}_{0}^{-1}(\delta\textsf{J}_{0})\textsf{J}_{0}^{-1}$.
Hence, we have the following sequence of identities (the index $0$ is dropped)

$\displaystyle\delta{\textsf{A}}=\delta(\textsf{J}^{-1}\textsf{J}^{-\top})$
$\displaystyle=\textsf{J}^{-1}(\delta\textsf{J}^{-1})^{\top}+(\delta\textsf{J}^{-1})\textsf{J}^{-\top}$
$\displaystyle=\textsf{J}^{-1}[-\textsf{J}^{-\top}(\delta\textsf{J})^{\top}\textsf{J}^{-\top}]+[-\textsf{J}^{-1}(\delta\textsf{J})\textsf{J}^{-1}]\textsf{J}^{-\top}$
$\displaystyle=-{\textsf{A}}(\delta\textsf{J})^{\top}\textsf{J}^{-\top}-\textsf{J}^{-1}(\delta\textsf{J}){\textsf{A}}.$

Here,

$\delta\textsf{J}=\left.\frac{d}{d\lambda}\left\\{\delta_{km}+{N}_{k}\frac{\partial(\lambda{\tilde{\rho}})}{\partial
y_{m}}+\frac{\partial{N}_{k}}{\partial
y_{m}}(\lambda{\tilde{\rho}})\right\\}\right|_{\lambda=0}={N}_{k}\frac{\partial{\tilde{\rho}}}{\partial
y_{m}}+\frac{\partial{N}_{k}}{\partial y_{m}}{\tilde{\rho}}.$

So,

$\delta{\textsf{A}}=-{\textsf{A}}(\nabla\otimes({\boldsymbol{\mathbf{N}}}{\tilde{\rho}}))\textsf{J}^{-\top}-\textsf{J}^{-1}(\nabla\otimes({\boldsymbol{\mathbf{N}}}{\tilde{\rho}}))^{\top}{\textsf{A}},$

or, in terms of its entries,

$\delta{\textsf{A}}_{ij}^{(0)}=-\sum_{k,m=1}\left({\textsf{A}}_{ik}^{(0)}\frac{\partial({N}_{m}{\tilde{\rho}})}{\partial
y_{k}}\textsf{J}_{0}^{jm}+\textsf{J}_{0}^{im}\frac{\partial({N}_{m}{\tilde{\rho}})}{\partial
y_{k}}{\textsf{A}}_{mj}^{(0)}\right).$

Next, let us compute the variation $\delta B_{0}$. To do so, let us first
compute the variations
$\delta(\boldsymbol{\mathbf{N}}\cdot\textsf{J}^{-\top}\nu_{\circ})$ and
$\delta\left(\textsf{A}^{(0)}\nu_{\circ}/(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ})\right)$,
and note that $\delta\boldsymbol{\mathbf{N}}\perp\nu_{\circ}$ and
$\boldsymbol{\mathbf{N}}\perp\delta\nu_{\circ}$:

$\displaystyle\delta(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ})$
$\displaystyle=\delta\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}+\boldsymbol{\mathbf{N}}\cdot\left(\delta\textsf{J}_{0}^{-\top}\nu_{\circ}+\textsf{J}_{0}^{-\top}\delta\nu_{\circ}\right)$
$\displaystyle=\left(-\textsf{J}_{0}^{-1}(\delta\textsf{J})\textsf{J}_{0}^{-1}\right)\nu_{\circ}^{\top}\boldsymbol{\mathbf{N}}$
$\displaystyle=-\sum_{i,k,m,j=1}\textsf{J}_{0}^{ki}\nu_{0k}\frac{\partial(N_{i}{\tilde{\rho}})}{\partial
y_{j}}\textsf{J}_{0}^{jm}N_{m}$
$\displaystyle=-\textsf{J}_{0}^{-\top}\nu_{\circ}\cdot(\nabla\otimes({\boldsymbol{\mathbf{N}}}{\tilde{\rho}}))^{\top}\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}.$

Therefore, we have the following computations

$\displaystyle\delta\left(\frac{\textsf{A}^{(0)}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}\right)$
$\displaystyle=\frac{\delta\textsf{A}^{(0)}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}-\textsf{A}^{(0)}\nu_{\circ}\frac{\delta(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ})}{(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ})^{2}}$
$\displaystyle=\frac{\left[-\textsf{A}^{(0)}(\nabla\otimes({\boldsymbol{\mathbf{N}}}{\tilde{\rho}}))\textsf{J}_{0}^{-\top}-\textsf{J}_{0}^{-1}(\nabla\otimes({\boldsymbol{\mathbf{N}}}{\tilde{\rho}}))^{\top}\textsf{A}^{(0)}\right]\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}$
$\displaystyle\qquad+\textsf{A}^{(0)}\nu_{\circ}\frac{\textsf{J}_{0}^{-\top}\nu_{\circ}\cdot(\nabla\otimes({\boldsymbol{\mathbf{N}}}{\tilde{\rho}}))^{\top}\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}}{(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ})^{2}}$
$\displaystyle=-\textsf{A}^{(0)}\nabla{\tilde{\rho}}-\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}\frac{\textsf{A}^{(0)}\nu_{\circ}\cdot\nabla{\tilde{\rho}}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}+\textsf{A}^{(0)}\nu_{\circ}\frac{\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}\cdot\nabla{\tilde{\rho}}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}$
$\displaystyle\qquad+\Bigg{[}-\frac{\textsf{A}^{(0)}(\nabla\otimes\boldsymbol{\mathbf{N}})\textsf{J}_{0}^{-\top}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}-\frac{\textsf{J}_{0}^{-1}(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{A}^{(0)}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}$
$\displaystyle\qquad\qquad+\frac{\textsf{A}^{(0)}\nu_{\circ}}{\left(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}\right)^{2}}\left(\textsf{J}_{0}^{-\top}\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}\right)\Bigg{]}{\tilde{\rho}}.$

Now, from (31), we recall that

$B_{\tilde{\rho}}:=B_{{\tilde{\rho}}}(y)=\left(\boldsymbol{\mathbf{N}}(y)\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}(y)\right)^{-1}\left(\nu_{\circ}(y)\cdot{{\textsf{A}}}\nu_{\circ}(y)\right),\qquad
y\in\Gamma.$

So, we have the following calculations

$\displaystyle\delta\left(\frac{\nu_{\circ}\cdot\textsf{A}^{(0)}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}\right)$
$\displaystyle=-\nu_{\circ}\cdot\textsf{A}^{(0)}\nabla{\tilde{\rho}}-\nu_{\circ}\cdot\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}\frac{\textsf{A}^{(0)}\nu_{\circ}\cdot\nabla{\tilde{\rho}}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}+\nu_{\circ}\cdot\textsf{A}^{(0)}\nu_{\circ}\frac{\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}\cdot\nabla{\tilde{\rho}}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}$
$\displaystyle\qquad+\Bigg{[}-\frac{\textsf{A}^{(0)}\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{J}_{0}^{-\top}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}-\frac{\textsf{J}_{0}^{-1}\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{A}^{(0)}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}$
$\displaystyle\qquad\qquad+\frac{\nu_{\circ}\cdot\textsf{A}^{(0)}\nu_{\circ}}{\left(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}\right)^{2}}\left(\textsf{J}_{0}^{-\top}\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}\right)\Bigg{]}{\tilde{\rho}}$
$\displaystyle=-2\textsf{A}^{(0)}\nu_{\circ}\cdot\nabla{\tilde{\rho}}-\frac{(\nu_{\circ}\cdot\textsf{A}^{(0)}\nu_{\circ})}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nabla{\tilde{\rho}}$
$\displaystyle\qquad+\Bigg{[}-\frac{\textsf{A}^{(0)}\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{J}_{0}^{-\top}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}-\frac{\textsf{J}_{0}^{-1}\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{A}^{(0)}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}$
$\displaystyle\qquad\qquad+\frac{\nu_{\circ}\cdot\textsf{A}^{(0)}\nu_{\circ}}{\left(\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}\right)^{2}}\left(\textsf{J}_{0}^{-\top}\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\textsf{J}_{0}^{-1}\boldsymbol{\mathbf{N}}\right)\Bigg{]}{\tilde{\rho}}$
$\displaystyle=:-2\textsf{A}^{(0)}\nu_{\circ}\cdot\nabla{\tilde{\rho}}-\frac{(\nu_{\circ}\cdot\textsf{A}^{(0)}\nu_{\circ})}{\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nu_{\circ}}\boldsymbol{\mathbf{N}}\cdot\textsf{J}_{0}^{-\top}\nabla{\tilde{\rho}}+h(\boldsymbol{\mathbf{N}},\partial_{y}\boldsymbol{\mathbf{N}}){\tilde{\rho}}.$

##  Appendix E Lemmata Proofs

For the benefit of the reader, we prove here some of the lemmas used in the
paper.

###  E.1 Proof of Lemma 2.5

In the proof we require a simple matrix inequality. For
$\textsf{A}=(a_{ij})\in\mathbb{R}^{d\times d}$, we define
$\left\|{\textsf{A}}\right\|=\sqrt{\sum_{i,j=1}^{d}{a_{ij}^{2}}}$ (i.e.,
$\left\|{\cdot}\right\|$ denotes the usual entry-wise matrix norm). Then, for
$\textsf{A},\textsf{B}\in\mathbb{R}^{d\times d}$, we have

$\left\|{\textsf{A}\textsf{B}}\right\|^{2}=\sum_{i=1}^{d}\sum_{j=1}^{d}\left(\sum_{k=1}^{d}a_{ik}b_{kj}\right)^{2}\leqslant\sum_{i=1}^{d}\sum_{j=1}^{d}\left(\sum_{k=1}^{d}a_{ik}^{2}\right)\left(\sum_{k=1}^{d}b_{kj}^{2}\right)=\left\|{\textsf{A}}\right\|^{2}\left\|{\textsf{B}}\right\|^{2}.$

So,
$\left\|{\textsf{A}\textsf{B}}\right\|\leqslant\left\|{\textsf{A}}\right\|\left\|{\textsf{B}}\right\|$.

###### Proof of Lemma 2.5.

Let $k\in\mathbb{N}$, $\alpha\in[0,1)$, and
$\phi\in{{C}}_{0}^{k+\alpha}(\Omega)$. Since
$\phi\in{{C}}_{0}^{k+\alpha}(\mathbb{R}^{d})$, we consider the case
$\Omega=\mathbb{R}^{d}$ without loss of generality. To show that
$\varphi\in\operatorname{Diffeo}^{k+\alpha}(\Omega,\Omega)$, we need to prove
the following:

  * (i)

$\varphi$ is injective;

  * (ii)

$\varphi$ is surjective;

  * (iii)

$\operatorname{{det}}(\nabla^{\top}\varphi)\neq 0$ for
$x\in\overline{\Omega}$.

Let us prove that $\varphi$ is injective. That is, for every $x,y\in\Omega$
such that $\varphi(x)=\varphi(y)$, we have $x=y$. So, let us suppose that
$x+\phi(x)=y+\phi(y)$. Then, we have the following sequence of computations:

$\displaystyle\left|{x-y}\right|=\left|{\phi(y)-\phi(x)}\right|$
$\displaystyle=\left|{\int_{0}^{1}\frac{d}{dt}\left[\phi(x+t(y-x))\right]\,dt}\right|$
$\displaystyle=\left|{\int_{0}^{1}\left[\nabla^{\top}\phi(x+t(y-x))\right](y-x)\,dt}\right|$
$\displaystyle\leqslant\int_{0}^{1}\left\|{\nabla^{\top}\phi(x+t(y-x))}\right\|\left|{y-x}\right|\,dt.$

Because $\max_{x\in\overline{\Omega}}\left\|{\nabla^{\top}\phi(x)}\right\|<1$,
then there exists ${\varepsilon_{\star}}>0$ such that
$\left\|{\nabla^{\top}\phi(x+t(y-x))}\right\|\leqslant{\varepsilon_{\star}}<1$.
Hence,

$\displaystyle\left|{x-y}\right|=\left|{\phi(y)-\phi(x)}\right|$
$\displaystyle\leqslant{\varepsilon_{\star}}\int_{0}^{1}\left|{y-x}\right|\,dt={\varepsilon_{\star}}\left|{y-x}\right|.$

This implies that $(1-{\varepsilon_{\star}})\left|{y-x}\right|=0$, or
equivalently, $y=x$. This proves the injectivity of $\varphi$.

Next let us prove that $\varphi$ is surjective. That is, for every
$y\in\mathbb{R}^{d}$, there exists $x\in\mathbb{R}^{d}$ such that
$\varphi(x)=y$. Let us define the map $T_{y}:\mathbb{R}^{d}\to\mathbb{R}^{d}$
where $T_{y}(x):=y-\phi(x)$. We want to solve the equation $T_{y}(x)=x$. Note
that

$\left|{T_{y}(x)-T_{y}(z)}\right|=\left|{\phi(x)-\phi(z)}\right|\leqslant{\varepsilon_{\star}}\left|{x-z}\right|.$

Because ${\varepsilon_{\star}}<1$, then by Contraction Mapping Theorem, there
exists an $x\in\mathbb{R}^{d}$ such that $x=T_{y}(x)$. This holds for
arbitrary $y\in\mathbb{R}^{d}$, proving the surjectivity of $\varphi$.

Lastly, let us show that $\operatorname{{det}}(\nabla^{\top}\varphi)\neq 0$ in
$\mathbb{R}^{d}$. Note that
$\nabla^{\top}\varphi=\textsf{I}+\nabla^{\top}\phi$, where
$\textsf{I}=\textsf{I}(x)=x$ is the identity matrix/operator. If
$\left\|{\nabla^{\top}\phi(x)}\right\|<1$, then
$\operatorname{{det}}(\textsf{I}+\nabla^{\top}\phi)\neq 0$. Let
$\textsf{A}(x):=\nabla^{\top}\phi(x)$,
$\textsf{A}\in{{C}}^{0}(\overline{\Omega},\mathbb{R}^{d})$,
$\left\|{\textsf{A}(x)}\right\|\leqslant{\varepsilon_{\star}}<1$, for all
$x\in\overline{\Omega}$, and
$\textsf{B}_{n}(x):=\sum_{k=0}^{n}(-\textsf{A}(x))^{k}$,
$\textsf{B}_{n}\in{{C}}^{0}(\overline{\Omega};\mathbb{R}^{d\times d})$ and
$\textsf{B}(x):=\sum_{k=0}^{\infty}(-\textsf{A}(x))^{k}$. Hence,

$\left\|{\textsf{B}_{n}(x)}\right\|\leqslant\sum_{k=0}^{n}\left\|{\textsf{A}(x)}\right\|^{k}\leqslant\sum_{k=0}^{n}\varepsilon_{\star}^{k}=\frac{1-\varepsilon_{\star}^{n+1}}{1-{\varepsilon_{\star}}}\leqslant\frac{1}{1-{\varepsilon_{\star}}},$

and $\\{\textsf{B}_{n}\\}$ is a Cauchy sequence in
${{C}}^{0}(\overline{\Omega};\mathbb{R}^{d\times d})$ and
$\textsf{B}(x)=\lim_{n\to\infty}\textsf{B}_{n}(x)$. Therefore,
$\textsf{B}\in{{C}}^{0}(\overline{\Omega};\mathbb{R}^{d\times d})$ and
$\lim_{n\to\infty}\textsf{B}_{n}=\textsf{B}$ in
${{C}}^{0}(\overline{\Omega};\mathbb{R}^{d\times d})$. These arguments show
that $\textsf{B}(x)=\sum_{n=0}^{\infty}(-\textsf{A}(x))^{n}$ uniformly
converges in $\overline{\Omega}$ where
$\textsf{B}\in{{C}}^{0}(\overline{\Omega};\mathbb{R}^{d\times d})$ and
$\textsf{B}(x)=(\textsf{I}+\textsf{A}(x))^{-1}$.

Now, it follows that,

$(\textsf{I}+\textsf{A}(x))\textsf{B}(x)=\lim_{n\to\infty}(\textsf{I}+\textsf{A}(x))\textsf{B}_{n}(x)=\lim_{n\to\infty}(\textsf{I}-(-\textsf{A}(x))^{n+1})=\textsf{I},$

or equivalently, $\textsf{I}+\textsf{A}(x)=\textsf{B}^{-1}(x)$, for all
$x\in\overline{\Omega}$. Taking the determinant of both sides of this
equation, together with the matrix inequality, we deduce that
$\operatorname{{det}}(\nabla^{\top}\varphi)\neq 0$ in $\mathbb{R}^{d}$.

To conclude that $\operatorname{{det}}(\nabla^{\top}\varphi)>0$ in
$\mathbb{R}^{d}$. We argue as follows. Consider the set
$\mathscr{O}:=\\{\textsf{A}\in\mathbb{R}^{d\times
d}\mid\left\|{\textsf{A}}\right\|<1\\}$. Clearly, $\mathscr{O}$ is a connected
convex open set. Let us define
$F(\textsf{A}):=\operatorname{{det}}(\textsf{I}+\textsf{A})$. We already prove
that $F(\textsf{A})\neq 0$ for all $\textsf{A}\in\mathscr{O}$ and
$F\in{{C}}^{0}(\mathscr{O};\mathbb{R})$. Evidently, the zero matrix
$\textsf{O}\in\mathscr{O}$ and
$F(\textsf{O})=\operatorname{{det}}\textsf{I}=1$. Hence, it follows that
$F(\textsf{A})>0$, for all $\textsf{A}\in\mathscr{O}$, as desired.

To finish the proof, let us note that, for every $x\in\Omega$, we have
$\varphi(x)\in\Omega$. Indeed, if $y=\varphi(x)\in\Omega^{c}$, we have
$\Omega\ni x=\varphi^{-1}(y)=y\in\Omega^{c}$, which is a contradiction. This
concludes the proof. ∎

###  E.2 Proof of Lemma 4.2

###### Proof of Lemma 4.2.

We first show that for any $y\in\overline{\Omega}$, we have that
$Z(y)\in\overline{\Omega(\rho)}$. We prove this by contradiction. So, suppose
the statement is not true, then there is a point $y_{a}\in\Omega$ such that
$Z(y_{a})\in\Gamma(\rho)\cup\Sigma$, $\Gamma(\rho)={\mathscr{S}}(\rho)$ for
sufficiently small $\rho>0$ (i.e.,
$\left\|{\rho}\right\|_{C^{2+\alpha}(\Gamma)}<\varepsilon$ where
$0<\varepsilon\ll 1$). We have the following situations:

  1. (i)

there is a point $z_{a}\in\Omega$ such that $Z(z_{a})\in\Omega(\rho)$;

  2. (ii)

there is a point $z_{b}\in\overline{\Omega}$ such that
$Z(z_{b})\not\in\overline{\Omega(\rho)}$;

  3. (iii)

there exists a curve $\mathscr{C}(z_{a},z_{b})$ such that
$\mathscr{C}(z_{a},z_{b})\setminus\\{z_{b}\\}\subset\Omega$;

  4. (iv)

there is a point $y_{a}\in\Omega\cap\mathscr{C}(z_{a},z_{b})$ such that
$Z(y_{a})\in\Gamma(\rho)\cup\Sigma$.

However, we can find a point $y_{b}\in\Gamma\cup\Sigma$ such that
$Z(y_{b})=Z(y_{a})$ because $Z$ is a ${{C}}^{2+\alpha}$-diffeomorphism from
$\partial\Omega$ onto $\partial\Omega(\rho)$. But, as we shall show next, the
map is injective, i.e., we have that $y_{b}=y_{a}$, a contradiction.

Now, let us verify that the map is one-to-one. We need to show that for any
point $y_{a},y_{b}\in\overline{\Omega}$ such that $Z(y_{b})=Z(y_{a})$, it must
be that $y_{b}=y_{a}$. Let us note that for any pair of points
$y_{a},y_{b}\in\overline{\Omega}$, there exists a curve
$\mathscr{C}(y_{a},y_{b})\subset\overline{\Omega}$ of class ${{C}}^{1}$ such
that $\left|{\mathscr{C}(y_{a},y_{b})}\right|\leqslant
c(\Omega)|y_{b}-y_{a}|$, for some constant $c(\Omega)>0$. Hence, the equation
$Z(y_{b})=Z(y_{a})$ implies that
$y_{b}+\boldsymbol{\mathbf{N}}(y_{b})\tilde{\rho}(y_{b})=y_{a}+\boldsymbol{\mathbf{N}}(y_{a})\tilde{\rho}(y_{a})$,
and we get the inequality

$\left|{y_{b}-y_{a}}\right|=\left|{\boldsymbol{\mathbf{N}}(y_{b})\tilde{\rho}(y_{b})-\boldsymbol{\mathbf{N}}(y_{a})\tilde{\rho}(y_{a})}\right|\leqslant
c(\Omega)\max_{y\in\overline{\Omega}}\left|{\nabla(\boldsymbol{\mathbf{N}}(y)\tilde{\rho}(y))}\right|\left|{y_{b}-y_{a}}\right|.$

By successively applying the inequality, we will get
$\left|{y_{b}-y_{a}}\right|=0$, or equivalently, $y_{b}=y_{a}$. This completes
the proof of the lemma. ∎

###  E.3 Proof of Lemma 7.2

###### Proof of Lemma 7.2.

The existence result follows from [54, Sec. 3, Thm. 4, p. 1270], where the
general linear problem with the time-derivative in the boundary condition was
studied. For the proof of the a priori estimate (53) we apply Schauder’s
method. That is, using the estimate for the solution to the corresponding
model problem in a half-space [42, Chap. 3, Sec. 2], one first obtain an
estimate of the form

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Theta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}\leqslant
c\left(\left|{\Theta}\right|^{(2)}_{{\overline{\Omega};\,[0,t]}}+\left|{\psi_{2}}\right|^{(2+\alpha)}_{{\Sigma;\,[0,t]}}+\left|{\Theta}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left|{\psi_{1}}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}\right),$

for some constant $c>0$. Here we apply the interpolation inequality (see,
e.g., [42, Eq. (2.1), p. 117] or [44, Chap. 6, Sec. 8])

$\left|{u}\right|_{\overline{\Omega}}^{(2)}\leqslant\varepsilon_{1}\sum_{\left|{k}\right|=2}|{{D}}^{k}u|_{\overline{\Omega}}^{(\alpha)}+c_{\varepsilon_{1}}{{\max_{\overline{\Omega}}}}\left|{u}\right|,$
(E.1)

where $\varepsilon_{1}>0$ is an arbitrary small number and
$c_{\varepsilon_{1}}\to\infty$ as $\varepsilon_{1}\to 0$, to
$\left|{\Theta}\right|^{(2)}_{{{{\overline{\Omega}}}}}$. We take
$\varepsilon_{1}$ small enough so that we will get an estimate of the form

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Theta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}\leqslant
c\left(\left|{\Theta}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left|{\psi_{1}}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left|{\psi_{2}}\right|^{(2+\alpha)}_{{\Sigma;\,[0,t]}}+\left|{\Theta}\right|^{\infty}_{\overline{\Omega};\,[0,t]}\right).$
(E.2)

Now, to get further with our computation, we note the following equalities for
$u\in{{C}}^{2+\alpha}(\overline{\Omega})$:

$\left\\{\begin{aligned}
\left|{u}\right|^{(1+\alpha)}_{\overline{\Omega}}&=\sum_{\left|{j}\right|<1+\alpha}{{\max_{\overline{\Omega}}}}\left|{{D}^{j}u(x)}\right|+[u]^{(1+\alpha)}_{\overline{\Omega}}\\\
&={{\max_{\overline{\Omega}}}}\left|{\Theta}\right|+\sum_{\left|{j}\right|=1}{{\max_{\overline{\Omega}}}}\left|{{D}^{j}u(x)}\right|+\sum_{\left|{j}\right|=1}\max_{x,{{\hat{x}}}\in\overline{\Omega}}\frac{\left|{{D}^{j}u(x)-{D}^{j}u(y)}\right|}{\left|{x-y}\right|^{\alpha}}\\\
[{D}u]^{(\alpha)}_{\overline{\Omega}}&=\sum_{\left|{j}\right|=[\alpha]}\max_{x,{{\hat{x}}}\in\overline{\Omega}}\frac{\left|{{D}^{j+1}u(x)-{D}^{j+1}u(y)}\right|}{\left|{x-y}\right|^{\alpha-[\alpha]}}\\\
&=\sum_{\left|{j}\right|=[1+\alpha]}\max_{x,{{\hat{x}}}\in\overline{\Omega}}\frac{\left|{{D}^{j}u(x)-{D}^{j}u(y)}\right|}{\left|{x-y}\right|^{1+\alpha-[1+\alpha]}}\\\
&\equiv[u]^{(1+\alpha)}_{\overline{\Omega}}.\end{aligned}\right.$ (E.3)

In above equalities, we apply the interpolation inequalities (see [42, Eq.
(2.1), p. 117])

$[{D}{u}]_{\overline{\Omega}}^{(\alpha)}\leqslant\varepsilon_{2}\left|{{u}}\right|_{\overline{\Omega}}^{(2)}+c_{\varepsilon_{2}}\left|{{u}}\right|_{\overline{\Omega}}^{(0)}\qquad\text{and}\qquad[{u}]_{\overline{\Omega}}^{(\alpha)}\leqslant\varepsilon_{3}\left|{{u}}\right|_{\overline{\Omega}}^{(1)}+c_{\varepsilon_{3}}\left|{{u}}\right|_{\overline{\Omega}}^{(0)},$
(E.4)

together with (E.1) to get an estimate for
$\left|{\Theta}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}$. We choose small
enough values for $\varepsilon_{2}>0$ and $\varepsilon_{3}>0$ such that the
right-hand side of the inequality (E.2) contains only the norms of known
functions and the quantity ${{\max_{\overline{\Omega}}}}|\Theta(\cdot,\tau)|$.
Since $\Theta$ is harmonic in $\Omega$, we can then apply the maximum
principle (see, e.g., [46, Chap. I.3, III, p. 7]) on the aforesaid norm to
obtain

$\left|{\Theta}\right|^{\infty}_{\overline{\Omega};\,[0,t]}\leqslant\left|{\psi_{2}}\right|^{\infty}_{\Sigma;\,[0,t]}+\left|{\Theta}\right|^{\infty}_{\Gamma;\,[0,t]}.$

We next estimate $\left|{\Theta}\right|^{\infty}_{\Gamma;\,[0,t]}$. Applying
the fundamental theorem of calculus, together with the homogenous initial
condition for the function $\Theta\big{|}_{\Gamma}$, we get

$\left|{\Theta}\right|^{\infty}_{\Gamma;\,[0,t]}=\max_{0\leqslant\tau\leqslant
t}\max_{\Gamma}\left|{\int_{0}^{\tau}{\frac{\partial}{\partial
s}}\Theta(\cdot,s){{{d}}s}}\right|\leqslant\int_{0}^{t}{{\max_{\Gamma}}}\left|{\frac{\partial}{\partial
s}}\Theta(\cdot,s)\right|{{{d}}s}.$ (E.5)

Inserting this estimate to (E.2), and letting
${{{\Psi}}(t)}=\left|{\psi_{1}}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left|{\psi_{2}}\right|^{(2+\alpha)}_{{\Sigma;\,[0,t]}}$,
leads to

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Theta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}\leqslant
c\left({{{\Psi}}(t)}+\int_{0}^{t}{{\max_{\Gamma}}}\left|{{\frac{\partial}{\partial
s}}\Theta(\cdot,s)}\right|{{{d}}s}\right).$ (E.6)

Because ${{{\Psi}}(t)}$ is non-negative, then we may apply Grönwall’s
inequality (see, e.g., [61, Appx. 5.1, p. 498]) to (E.6) deducing the desired
estimate (53). ∎

###  E.4 Proof of Lemma 8.1

We provide here the proof of Lemma 8.1.

###### Proof of Lemma 8.1.

For some $\alpha\in(0,1)$, let $\Omega\subset\mathbb{R}^{d}$ be an open,
connected, bounded set of class ${{C}}^{2+\alpha}$ and
$u\in{{C}}^{2+\alpha}(\overline{\Omega})$. Then, $\Omega$ satisfies the cone
condition and, in particular, the (interior/exterior) ball condition. Let us
suppose that, in specific, $\Omega$ has cone property with the cone having the
opening angle $\hat{\theta}:=\hat{\theta}(\Omega)$ and height $h$. To prove
estimates (64) and (65), it is enough to show that the following inequalities
are true:

$\displaystyle[u]_{\overline{\Omega}}^{(1)}$
$\displaystyle\leqslant\varepsilon^{1+\alpha}[u]_{\overline{\Omega}}^{(2+\alpha)}+\frac{c_{7}}{\varepsilon}\left|{u}\right|_{\overline{\Omega}}^{(0)},$
(E.7) $\displaystyle[u]_{\overline{\Omega}}^{(2)}$
$\displaystyle\leqslant\varepsilon^{\alpha}[u]_{\overline{\Omega}}^{(2+\alpha)}+\frac{c_{8}}{\varepsilon^{2}}\left|{u}\right|_{\overline{\Omega}}^{(0)}.$
(E.8)

We will only verify in detail inequality (E.8), but (E.7) can be shown in a
similar fashion. To this end, let us consider another cone
${\pazocal{{C}}}_{1}:={\pazocal{{C}}}_{1}(\hat{\theta},1)$ with the same
opening angle, but with cone height equal to $1$. Using the interpolation
inequality (E.15) (with $\varepsilon=1$) in Lemma E.1, we get

$[u]_{{\pazocal{{C}}}_{1}}^{(2)}\leqslant[u]_{{\pazocal{{C}}}_{1}}^{(2+\alpha)}+c(d,\hat{\theta})\left|{u}\right|_{{\pazocal{{C}}}_{1}}^{(0)},$
(E.9)

for any function $u\in{{C}}^{2+\alpha}(\overline{{\pazocal{{C}}}}_{1})$, where
$c(d,\hat{\theta})>0$ is a constant that depends on $d$ and $\hat{\theta}$. To
get the desired coefficients in the inequality condition, we observe that for
any constant $\tilde{\mu}>0$ and function
$\tilde{u}(\tilde{x})=u(\tilde{\mu}x)$ we have the equality

$[\tilde{u}]_{\overline{\Omega}}^{(l)}=\tilde{\mu}^{l}[u]_{\overline{\Omega}}^{(l)},$
(E.10)

for all $0\leqslant l\leqslant 2+\alpha$. So, for
$u\in{{C}}^{2+\alpha}({\overline{\pazocal{{C}}}}_{\varepsilon})$, where the
cone
${{\pazocal{{C}}}_{\varepsilon}}:={{\pazocal{{C}}}_{\varepsilon}}(\hat{\theta},\varepsilon)$
with its vertex at the origin, we apply the change of variables
$\tilde{x}=x/\varepsilon$, $0<\varepsilon\leqslant h$, and consider the
function $\tilde{u}(\tilde{x})=u(\varepsilon\tilde{x})=u(x)$. Obviously,
$\tilde{u}\in{{C}}^{2+\alpha}(\overline{{\pazocal{{C}}}}_{1})$, and (E.9) also
holds for $\tilde{u}$. Applying identity (E.10) to our case, we get

$[\tilde{u}]_{{\pazocal{{C}}}_{\varepsilon}}^{(2)}\leqslant[\tilde{u}]_{{\pazocal{{C}}}_{\varepsilon}}^{(2+\alpha)}+\tilde{c}(d,\hat{\theta})|\tilde{u}|_{{\pazocal{{C}}}_{\varepsilon}}^{(0)}\quad\Longleftrightarrow\quad[u]_{{\pazocal{{C}}}_{\varepsilon}}^{(2)}\leqslant\varepsilon^{\alpha}[u]_{{\pazocal{{C}}}_{\varepsilon}}^{(2+\alpha)}+\frac{\tilde{c}(d,\hat{\theta})}{\varepsilon^{2}}\left|{u}\right|_{{\pazocal{{C}}}_{\varepsilon}}^{(0)},$
(E.11)

for some constant $\tilde{c}(d,\hat{\theta})>0$. Now, note that for any
$x\in\Omega\in{{C}}^{2+\alpha}$, we can construct a cone
${{\pazocal{{C}}}_{\varepsilon}}$ with vertex $x$ and such that
${{\pazocal{{C}}}_{\varepsilon}}\subset\Omega$. Hence, using (E.11), we can
further get the estimate

$[u]_{{\pazocal{{C}}}_{\varepsilon}}^{(2)}\leqslant\varepsilon^{\alpha}[u]_{\overline{\Omega}}^{(2+\alpha)}+\frac{\tilde{c}}{\varepsilon^{2}}\left|{u}\right|_{\overline{\Omega}}^{(0)}.$

Because $x\in\Omega$ is taken arbitrarily, we conclude that

$[u]_{\overline{\Omega}}^{(2)}\leqslant\varepsilon^{\alpha}[u]_{\overline{\Omega}}^{(2+\alpha)}+\frac{c_{7}}{\varepsilon^{2}}\left|{u}\right|_{\overline{\Omega}}^{(0)},$

for some constant $c_{7}>0$, proving (E.8). Using the definition of the norm
$\left|{\cdot}\right|_{\overline{\Omega}}^{(l)}$, we finally recover
inequality (65). The other interpolation inequality (E.7) from which (64)
relies can be shown in a similar manner as above with the help of the
interpolation inequality (E.14) in Lemma E.1 (stated in the subsection E.5
below). ∎

###  E.5 Interpolation inequalities

We justify in this subsection the interpolation inequalities used to proved
Lemma 8.1. Similar and more general interpolation inequalities are given and
proved in [44, Sec. 6.8, pp. 130–136] (see also [49, Cor. 5.3.4, p. 144] for
analogous results).

######  Lemma E.1.

Let $B_{r}$ be a ball in $\mathbb{R}^{d}$ with radius $r$ and consider a
function $u\in{{C}}^{1+\alpha}(\overline{B}_{r})$, $\alpha\in(0,1)$. For any
real number $\varepsilon\in(0,r]$, the following inequalities hold:

$\displaystyle[u]_{B_{r}}^{(1)}$
$\displaystyle\leqslant\varepsilon^{\alpha}[u]_{B_{r}}^{(1+\alpha)}+\frac{{{c}_{d}}}{\varepsilon}\left|{u}\right|_{B_{r}}^{(0)},$
(E.12) $\displaystyle[u]_{B_{r}}^{(\alpha)}$
$\displaystyle\leqslant\varepsilon[u]_{B_{r}}^{(1+\alpha)}+\frac{{{c}_{d}}}{\varepsilon^{\alpha}}\left|{u}\right|_{B_{r}}^{(0)},$
(E.13)

where ${{c}_{d}}$ is a positive constant that depends only on the dimension
$d$. Moreover, if $u\in{{C}}^{2+\alpha}(\overline{B}_{r})$, then we further
have

$\displaystyle[u]_{B_{r}}^{(1)}$
$\displaystyle\leqslant\varepsilon^{(1+\alpha)}[u]_{B_{r}}^{(2+\alpha)}+\frac{{{c}_{d}}}{\varepsilon}\left|{u}\right|_{B_{r}}^{(0)},$
(E.14) $\displaystyle[u]_{B_{r}}^{(2)}$
$\displaystyle\leqslant\varepsilon^{\alpha}[u]_{B_{r}}^{(2+\alpha)}+\frac{{{c}_{d}}}{\varepsilon^{2}}\left|{u}\right|_{B_{r}}^{(0)}.$
(E.15)

###### Proof.

Consider a ball $B_{r}\subset\mathbb{R}^{d}$ with radius $r$ and let $x\in
B_{r}$. Let us choose a point in $B_{r}$, say $x_{0}$, such that $x\in
B_{\varepsilon/{2}}(x_{0})\subset B_{r}$. Over $B_{\varepsilon/{2}}(x_{0})$,
we integrate $\partial_{x_{i}}u$ and apply Green’s formula to obtain

$\int_{B_{\varepsilon/{2}}(x_{0})}{\partial_{x_{i}}u}{\,dx}=\int_{\partial
B_{\varepsilon/{2}}(x_{0})}{u\cos(\nu_{\partial
B_{\varepsilon/{2}}},x_{i})}{\,ds},$ (E.16)

where $\nu_{\partial B_{\varepsilon/{2}}}$ is the outward unit normal vector
to the boundary $\partial B_{\varepsilon/{2}}$. Because
$u\in{{C}}^{1+\alpha}(\overline{B}_{r})$, then the mean value theorem can be
applied and it tells us that there is a point $\bar{x}\in B_{\varepsilon/{2}}$
such that

$\partial_{x_{i}}u(\bar{x})=\frac{1}{\left|{B_{\varepsilon/{2}}}\right|}\int_{B_{\varepsilon/{2}}(x_{0})}{\partial_{x_{i}}u}{\,dx},$

from which, together with (E.16), we can obtain the estimate (cf. [44, Eq.
(2.1), p. 22])

$\displaystyle\left|{\partial_{x_{i}}u(\bar{x})}\right|=\dfrac{1}{\left|{B_{\varepsilon/{2}}}\right|}\left|{\int_{B_{\varepsilon/{2}}(x_{0})}{\partial_{x_{i}}u}{\,dx}}\right|\leqslant\dfrac{\left|{\partial
B_{\varepsilon/{2}}}\right|}{\left|{B_{\varepsilon/{2}}}\right|}\leqslant\dfrac{\dfrac{2\pi^{d/2}}{\Gamma_{f}({d}/{2})}\left(\dfrac{\varepsilon}{2}\right)^{d-1}}{\dfrac{\pi^{d/2}}{\Gamma_{f}({d}/{2}+1)}\left(\dfrac{\varepsilon}{2}\right)^{d}}\left|{u}\right|_{B_{r}}^{(0)}=\frac{2d}{\varepsilon}\left|{u}\right|_{B_{r}}^{(0)},$

where $\Gamma_{f}$ denotes the gamma-function. For $\bar{x}\neq x$, we can
then make the following estimation of $\left|{\partial_{x_{i}}u(x)}\right|$:

$\displaystyle\left|{\partial_{x_{i}}u(x)}\right|$
$\displaystyle\leqslant\left|{\partial_{x_{i}}u(x)-\partial_{x_{i}}u(\bar{x})}\right|+\left|{\partial_{x_{i}}u(\bar{x})}\right|\leqslant\frac{|\partial_{x_{i}}u(x)-\partial_{x_{i}}u(\bar{x})|}{\left|{x-\bar{x}}\right|^{\alpha}}\left|{x-\bar{x}}\right|^{\alpha}+\frac{2d}{\varepsilon}\left|{u}\right|_{B_{r}}^{(0)}$
$\displaystyle\leqslant\varepsilon^{\alpha}[\partial_{x_{i}}u]_{B_{r}}^{(\alpha)}+\frac{2d}{\varepsilon}\left|{u}\right|_{B_{r}}^{(0)}.$

This verifies estimate (E.12).

Next, we validate (E.13) using (E.12). For any
$x,{{\hat{x}}}\in\overline{B}_{r}$, $x\neq{{\hat{x}}}$, we have either
$\left|{x-{{\hat{x}}}}\right|<\varepsilon$ or
$\left|{x-{{\hat{x}}}}\right|\geqslant\varepsilon$ which would give us the
inequalities

$\frac{|u(x)-u({{\hat{x}}})|}{\left|{x-{{\hat{x}}}}\right|^{\alpha}}=\left|{x-{{\hat{x}}}}\right|^{1-\alpha}\frac{|u(x)-u({{\hat{x}}})|}{\left|{x-{{\hat{x}}}}\right|}<\varepsilon^{1-\alpha}\left|{u}\right|_{B_{r}}^{(1)},$

and

$\frac{|u(x)-u({{\hat{x}}})|}{\left|{x-{{\hat{x}}}}\right|^{\alpha}}\leqslant\frac{1}{\varepsilon^{\alpha}}\left(\max_{x\in\Omega}|u(x)|+\max_{{{\hat{x}}}\in\Omega}|u({{\hat{x}}})|\right)=\frac{2}{\varepsilon^{\alpha}}\left|{u}\right|_{B_{r}}^{(0)},$

respectively. In either case, we get the estimate

$\frac{|u(x)-u({{\hat{x}}})|}{\left|{x-{{\hat{x}}}}\right|^{\alpha}}\leqslant\varepsilon^{1-\alpha}\left|{u}\right|_{B_{r}}^{(1)}+\frac{2}{\varepsilon^{\alpha}}\left|{u}\right|_{B_{r}}^{(0)}.$

Combining the above estimate with (E.12), we obtain (E.13).

We next prove (E.14) and (E.15) using (E.12) and (E.13), respectively. To this
end, let us assume that $u\in{{C}}^{2+\alpha}(\overline{B}_{r})$. Then, (E.13)
implies that

$[u]_{B_{r}}^{(1+\alpha)}\leqslant\varepsilon_{1}[u]_{B_{r}}^{(2+\alpha)}+\frac{{{c}_{d}}}{\varepsilon_{1}^{\alpha}}\left|{u}\right|_{B_{r}}^{(1)},$
(E.17)

for any $\varepsilon_{1}\in(0,r]$. Let us denote $\varepsilon=\varepsilon_{2}$
in (E.12) and use the resulting inequality to estimate the quantity
$\left|{u}\right|_{B_{r}}^{(1)}$ in (E.17). These lead us to

$\left(1-{{c}_{d}}\left(\frac{\varepsilon_{2}}{\varepsilon_{1}}\right)^{\alpha}\right)[u]_{B_{r}}^{(1+\alpha)}\leqslant\varepsilon_{1}[u]_{B_{r}}^{(2+\alpha)}+\frac{{{c}_{d}}^{2}}{\varepsilon_{1}^{\alpha}\varepsilon_{2}}\left|{u}\right|_{B_{r}}^{(1)}.$
(E.18)

To force a single parameter $\varepsilon\in(0,r]$ in the above inequality, we
apply a scaling technique. To this end, we introduce the scaling
$\varepsilon_{1}=t_{1}\varepsilon$ and $\varepsilon_{2}=t_{2}\varepsilon$,
where $t_{1},t_{2}\in(0,1)$ and substitute these quantities to (E.18) to
obtain

$\frac{1}{t_{1}}\left(1-{{c}_{d}}\left(\frac{t_{2}}{t_{1}}\right)^{\alpha}\right)[u]_{B_{r}}^{(1+\alpha)}\leqslant\varepsilon[u]_{B_{r}}^{(2+\alpha)}+\frac{{{c}_{d}}^{2}}{t_{1}^{1+\alpha}t_{2}\varepsilon^{1+\alpha}}\left|{u}\right|_{B_{r}}^{(0)}.$
(E.19)

Now, we only need to choose the values of $t_{1}$ and $t_{2}$ such that
$t_{1}^{-1}\left(1-{{c}_{d}}\left(t_{2}t_{1}^{-1}\right)^{\alpha}\right)=1$
and denote by
${{\tilde{c}}_{d}}:={{c}_{d}}^{2}({t_{1}^{1+\alpha}t_{2}})^{-1}>0$ to obtain

$\varepsilon^{\alpha}[u]_{B_{r}}^{(1+\alpha)}\leqslant\varepsilon^{1+\alpha}[u]_{B_{r}}^{(2+\alpha)}+{{\tilde{c}}_{d}}\left|{u}\right|_{B_{r}}^{(0)}.$

Using the above estimate for the quantity
$\varepsilon^{\alpha}[u]_{B_{r}}^{(1+\alpha)}$ in (E.12), and after some
change of notations, we finally get (E.14). Of course, the inequality holds
for any $\varepsilon\in(0,r]$ since we choose $t_{1},t_{2}\in(0,1)$.

We apply the same technique to prove (E.15) using (E.12) and (E.14). That is,
using (E.12) with $u\in{{C}}^{2+\alpha}(\overline{B}_{r})$, and then applying
(E.14) to estimate $\varepsilon\left|{u}\right|_{B_{r}}^{(1)}$, we obtain

$\displaystyle\varepsilon^{2}[u]_{B_{r}}^{(2)}\leqslant\varepsilon^{2+\alpha}[u]_{B_{r}}^{(2+\alpha)}+{{c}_{d}}\varepsilon\left|{u}\right|_{B_{r}}^{(1)}\leqslant(1+{{c}_{d}})\varepsilon^{2+\alpha}[u]_{B_{r}}^{(2+\alpha)}+{{c}_{d}}^{2}\left|{u}\right|_{B_{r}}^{(0)}.$

Again, to force the constant, which is dependent only on the dimension $d$, to
appear only as a coefficient of $\left|{u}\right|_{B_{r}}^{(0)}$, we consider
the scaled free parameter $t\tilde{\varepsilon}=\varepsilon\in(0,r]$, where
$t\in(0,1)$, to obtain

$\tilde{\varepsilon}^{2}[u]_{B_{r}}^{(2)}\leqslant(1+{{c}_{d}})t^{\alpha}\tilde{\varepsilon}^{2+\alpha}[u]_{B_{r}}^{(2+\alpha)}+\left(\frac{{{c}_{d}}}{t}\right)^{2}\left|{u}\right|_{B_{r}}^{(0)}.$

We choose $t\in(0,1)$ such that $(1+{{c}_{d}})t^{\alpha}=1$ and let
${{\tilde{c}}_{d}}:=\left({{c}_{d}}t^{-1}\right)^{2}$ in order to finally get

$[u]_{B_{r}}^{(2)}\leqslant\tilde{\varepsilon}^{\alpha}[u]_{B_{r}}^{(2+\alpha)}+\frac{{{\tilde{c}}_{d}}}{\tilde{\varepsilon}^{2}}\left|{u}\right|_{B_{r}}^{(0)},$

which holds for all $\tilde{\varepsilon}\in(0,r]$. This verifies the last
estimate (E.15), completing the proof of the lemma. ∎

## References

  * [1] M. C. Delfour, J.-P. Zolésio, Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization, 2nd Edition, Vol. 22 of Adv. Des. Control, SIAM, PA, 2011.
  * [2] A. Henrot, M. Pierre, Shape Variation and Optimization: A Geometrical Analysis, Vol. 28 of Tracts in Mathematics, European Mathematical Society, Zürich, 2018\.
  * [3] J. Sokołowski, J.-P. Zolésio, Introduction to Shape Optimization, Vol. 16 of Springer Series in Computational Mathematics, Springer, Berlin, Heidelberg, 1992.
  * [4] K. Eppler, H. Harbrecht, A regularized newton method in electrical impedance tomography using shape Hessian information, Control Cybernet. 34 (1) (2005) 203–225.
  * [5] J. R. Roche, J. Sokołowski, Numerical methods for shape identification problems, Control Cybernet. 25 (5) (1996) 867–895.
  * [6] L. Afraites, A new coupled complex boundary method (CCBM) for an inverse obstacle problem, Discrete Contin. Dyn. Syst. Ser. S 15 (1) (2022) 23 – 40.
  * [7] I. Akduman, R. Kress, Electrostatic imaging via conformal mapping, Inverse Problems 18 (2002).
  * [8] G. Alessandrini, V. Isakov, J. Powell, Local uniqueness in the inverse problem with one measurement, Trans. Am. Math. Soc. 347 (1995) 3031–3041.
  * [9] G. Alessandrini, A. D. Valenzuela, Unique determination of multiple cracks by two measurements, SIAM J. Control Optim. 34 (1996) 913–921.
  * [10] L. Bourgeois, J. Dardé, A quasi-reversibility approach to solve the inverse obstacle problem, Inverse Prob. Imaging 4 (2010) 351–377.
  * [11] R. Chapko, R. Kress, A hybrid method for inverse boundary value problems in potential theory, J. Inv. Ill-Posed Problems 13 (1) (2005) 27–40.
  * [12] F. Hettlich, W. Rundell, The determination of a discontinuity in a conductivity from a single boundary measurement, Inverse Problems 14 (1998) 67–82.
  * [13] V. Isakov, Inverse Problems for Partial Differential Equations, Vol. 127, Springer Business & Media, 2066.
  * [14] G. I. Bizhanova, V. A. Solonnikov, On free boundary problems for the second order parabolic equationson free boundary problems for the second order parabolic equations, Algebra Anal. 12 (6) (2000) 98–139.
  * [15] V. A. Solonnikov, Lectures on evolution free boundary problems: classical solutions, Lect. Notes Math., Springer, 2003, pp. 123–175.
  * [16] J. Escher, G. Simonett, Classical solutions of multidimensional hele–shaw models, SIAM J. Math. Anal. 28 (5) (1997) 1028–1047.
  * [17] P. Cardaliaguet, O. Ley, Some flows in shape optimization, Arch. Rational Mech. Anal. 183 (2007) 21–58.
  * [18] P. Cardaliaguet, E. Rouy, Viscosity solutions of increasing flows of sets. application of the Hele-Shaw problem for power-law fluids, SIAM J. Math. Anal. 38 (1) (2006) 143–165.
  * [19] Y. Giga, Surface Evolution Equations A Level Set Approach, Vol. 99 of Monographs in Mathematics, Birkhäuser, Basel, 2006.
  * [20] P. I. Plotnikov, J. Sokołowski, Geometric aspects of shape optimization, J. Geom. Anal. 33 (2023) Article 206.
  * [21] P. I. Plotnikov, J. Sokołowski, Gradient flow for kohn–vogelius functional, Siberian Electron. Math. Rep. (2023) to appear.
  * [22] L. Hörmander, The Analysis of Linear Partial Differential Operators, Springer, New York, 1983–1985.
  * [23] V. Šverák, On optimal shape design, J. Math. Pures Appl. 72 (1993) 537–551.
  * [24] A. Henrot, W. Horn, J. Sokołowski, Domain optimization problem for stationary heat equation. , 6, 2, ., Appl. Math. & Comp. Sci. 6 (2) (1996) 353–374.
  * [25] D. Bucur, J.-P. Zolésio, $n$-dimensional shape optimization under the capacitary constraints, J. Differ. Equ. 123 (2) (1995) 504–522.
  * [26] L. Afraites, J. F. T. Rabago, Boundary shape reconstruction with robin condition: existence result, stability analysis, and inversion via multiple measurements, Comput. Appl. Math. 43 (2024) Art. 270.
  * [27] J. Crank, Free and Moving Boundary Problems, Clarendon press, 1984.
  * [28] C. M. Elliot, J. R. Ockendon, Weak and Variational Methods for Moving Boundary Problems, Pitman, Boston, 1982.
  * [29] S. Richardson, Hele-Shaw flows with a free boundary produced by the injection of the fluid into a narrow channel, J. Fluid Mech. 56 (1972) 609–618.
  * [30] J. Neuberger, Sobolev Gradients and Differential Equations, 2nd Edition, Vol. 1670 of Lecture Notes in Mathematics, Springer-Verlag, Berlin, Heidelberg, 2010\.
  * [31] H. Azegami, Shape Optimization Problems, Springer Optimization and its Applications, Springer, Singapore, 2020.
  * [32] G. Dogan, P. Morin, R. H. Nochetto, M. Verani, Discrete gradient flows for shape optimization and applications, Comput. Methods Appl. Mech. Engrg. 196 (37–40) (2007) 3898–3914.
  * [33] B. Mohammadi, O. Pironneau, Applied Shape Optimization for Fluids, Clarendon press, Oxford, 2001.
  * [34] L. Afraites, M. Dambrine, K. Eppler, D. Kateb, Detecting perfectly insulated obstacles by shape optimization techniques of order two, Discrete Contin. Dyn. Syst. B 8 (2) (2007) 389–416.
  * [35] J. F. T. Rabago, H. Azegami, Shape optimization approach to defect-shape identification with convective boundary condition via partial boundary measurement, Japan J. Indust. Appl. Math 36 (1) (2019) 131–176.
  * [36] Y. Sunayama, M. Kimura, J. F. T. Rabago, Comoving mesh method for certain classes of moving boundary problems, Japan J. Indust. Appl. Math 39 (2022) 973–1001.
  * [37] Y. Sunayama, J. F. T. Rabago, M. Kimura, Comoving mesh method for multi-dimensional moving boundary problems: Mean-curvature flow and Stefan problems, Math. Comput. Simul. 221 (2024) 589–605.
  * [38] J. F. T. Rabago, H. Azegami, A second-order shape optimization algorithm for solving the exterior Bernoulli free boundary problem using a new boundary cost functional, Comput. Optim. Appl. 77 (1) (2020) 251–305.
  * [39] D. Colton, R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, 4th Edition, Springer-Verlag, New York, 2019.
  * [40] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (2012) 251–265.
  * [41] N. V. Krylov, Lectures on Elliptic and Parabolic Equations in Holder Spaces, Vol. 12 of Graduate Studies in Mathematics, American Mathematical Society, Providence, 1997.
  * [42] O. A. Ladyženskaja, N. N. Ural’ceva, Equations aux Dérivées Partielles de Type Elliptique, Dunod, Paris, 1968.
  * [43] E. V. Frolova, Quasistationary approximation for the stefan problem, J. Math. Sci. 132 (4) (2006) 562–575.
  * [44] D. Gilbarg, N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, 2nd Edition, Springer-Verlag, Berlin, Heidelberg, 2001.
  * [45] L. C. Evans, Partial Differential Equations, 1st Edition, Vol. 19 of Graduate Series in Mathematics, AMS, 1998.
  * [46] C. Miranda, Partial Differential Equations of Elliptic Type, 2nd Edition, Ergebnisse der Mathematik und ihrer Grenzgebiete 2, Springer-Verlag, Berlin, Heidelberg, 1970.
  * [47] O. D. Kellog, On the derivates of harmonic functions on the boundary, Trans. Mat. Soc. 33 (1931) 486–510.
  * [48] J. Schauder, Über lineare elliptische differentialgleichungen zweiter ordnung, Math. Z. 38 (1934) 257–282.
  * [49] R. Fiorenza, Hölder and locally Hölder continuous functions, and open sets of class $C^{k}$, $C^{k,\lambda}$, Birkhaüser, 2017.
  * [50] G. Nardi, Schauder estimate for solutions of Poisson’s equation with Neumann boundary condition, L’Enseign. Math. 60 (2014) 421–435.
  * [51] B. K. Driver, Analysis with Applications Part 2, Springer, 2003.
  * [52] G. M. Troianiello, Elliptic Differential Equations and Obstacle Problems, Plenum Press, New York, 1987.
  * [53] B. K. Driver, Analysis with Applications Part 1, Springer, 2003.
  * [54] S. N. Antontsev, C. R. Gonçalves, A. M. Meirmanov, Exact estimates for the classical solutions to the free boundary problem in the hele-shaw cell, Adv. Diff. Equ. 8 (10) (2003) 1259–1280.
  * [55] V. Volpert, Elliptic Partial Differential Equations Volume 2: Reaction-Diffusion Equations, Birkhäuser, 2014.
  * [56] M. R. Ebert, M. Reissig, Methods for Partial Differential Equations. Qualitative Properties of Solutions, Phase Space Analysis, Semilinear Models, Birkhäuser, Cham, 2018.
  * [57] P. Hartman, Ordinary Differential Equations, 2nd Edition, SIAM, 2002.
  * [58] W. W. Hager, Updating the inverse of a matrix, SIAM Review 31 (2) (1989) 221–239.
  * [59] J. Sherman, W. J. Morrison, Adjustment of an inverse matrix corresponding to changes in the elements of a given column or a given row of the original matrix, Ann. Math. Stat. 20 (621) (1949).
  * [60] M. A. Woodbury, The stability of out-input matrices, Chicago III (1949) 5 pp.
  * [61] S. N. Ethier, T. G. Kurtz, Markov processes, characterization and convergence, Probability and Statistics, Wiley, 2005.