Higher-order Riemannian splines and the interpolation problem: an approach of gradient flows

Chun-Chi Lin Department of Mathematics, National Taiwan Normal University, Taipei, 116 Taiwan<EMAIL_ADDRESS>Dung The Tran Department of Mathematics, National Taiwan Normal University, Taipei, 116 Taiwan;
Current Address: Faculty of Mathematics-Mechanics-Informatics, Vietnam National University, University of Science, Hanoi, Vietnam<EMAIL_ADDRESS>
Keywords: higher-order Riemannian spline,
geometric $k$-spline, Riemannian cubic, spline interpolation

MSC(2020): primary  35K52, 49J20;  secondary  41A15, 35K25

In this article, we demonstrate that spline interpolation problems on smooth Riemannian manifolds can be resolved by applying the method of gradient flows proposed here. Our approach appears to be novel in addressing such problems and can be extended to the spline interpolation problem on Lie groups, which are commonly implicated in formulations of mechanical optimal control theory. This approach seems innovative in both optimal and geometric control theory. We provide the existence of global solutions in Hölder spaces for the gradient flow, whose asymptotic limits establish the existence of Riemannian splines, thereby offering solutions to spline interpolation problems on smooth Riemannian manifolds.

=13.5 pt plus .5pt minus .5pt

§ INTRODUCTION

The significance of data-fitting problems in curved spaces emerged in the late 1980s, spurred by advancements in mechanical and robotics industries. Following groundbreaking work on spherical smoothing splines, there has been a growing recognition of the challenges associated with understanding functions defined on non-Euclidean spaces for accurate interpretation of given data points. This challenge has reverberated across diverse research areas, including image processing, trajectory planning, and data analysis. Over the past few decades, researchers have delved into alternative approaches to extend interpolation and approximation methods to the Riemannian setting. This interdisciplinary exploration spans calculus of variations, optimal control, geophysics, geometric mechanics, robotics, quantum information processing, computer graphics, medical imaging, and data science. The complexity of these problems serves as a catalyst for investigating interpolation challenges in curved spaces. Notably, in contrast to Euclidean spaces, the primary hurdle in generalizing classical fitting or interpolation problems to curved spaces lies in the absence of a straightforward understanding of polynomial curve expressions in such spaces.

Geometric polynomials are conventionally defined using a variational approach, with recent emphasis on Riemannian cubics as a means to extend spline-based methods to manifolds in physics and engineering. However, constructing splines using Riemannian cubics poses a substantial challenge.
Consequently, proving the existence of Riemannian polynomial or cubic splines has garnered interest and remains an ongoing research issue. Riemannian splines find extensive use in physical and applied sciences, such as Riemannian cubics for rigid body control [19, 20] and geometric $k$-splines for computational anatomy [9, 10].
These splines generalize polynomial splines in Euclidean spaces to splines in Lie groups or Riemannian manifolds. The generalization was first achieved by Gabriel and Kajiya in [8],
and by Noakes, Heinzinger, and Paden in [20]
for cubic splines ($k=2$), and by Camarinha, Silva Leite, and Crouch in [3] for polynomial splines ($k\ge2$).
In their approach, splines are defined as critical points of the higher-order energy functional

\sum_{l=1}^{q}\int_{x_{l-1}}^{x_l} \, \bigg\langle \frac{D^{k-1}}{dx^{k-1}} \gamma_x, \frac{D^{k-1}}{dx^{k-1}} \gamma_x \bigg\rangle \, dx,

where $\gamma:[x_{0},x_{q}]\rightarrow M$ is a sufficiently smooth map into a Riemannian manifold $M$ such that $\gamma(x_j)=p_j$ for all $j$,
$\{p_0,\cdots,p_q\}\subset M$ is a set of prescribed ordered points,
$\gamma_x$ is the tangent vector of the map $\gamma$,
$\frac{D}{dx}$ is the (directional) covariant differentiation on $M$, and
$\langle \cdot,\cdot\rangle$ is the inner product induced from the metric of $M$.

In this article, we demonstrate that spline interpolation problems on smooth Riemannian manifolds can be resolved by applying the method of gradient flows.
Our approach appears to be novel in addressing such problems and can be extended to spline interpolation problems in Lie groups, which are commonly implicated in mechanical optimal control theory formulations.
This approach seems innovative in both optimal and geometric control theory.

§.§ Previous results

In optimal control theory, the existence of solutions is usually obtained by working with the Euler-Lagrange equations and applying Pontryagin's maximum principle. The result on the existence of global solutions for Riemannian splines in the literature is limited to cases with certain assumptions. For the case where $q=1$ and $M$ is a complete Riemannian manifold, the proofs on the existence of (global) weak solutions for the Riemannian splines using the min-max method are provided by Giambò, Giannoni, and Piccione in [11] for $k=2$ and in [12] for $k\geq 2$. When $q=1$ and $M$ is a Riemannian manifold, Camarinha, Silva Leite, and Crouch in [4] use biexponential maps to prove the local existence and uniqueness of Riemannian cubics subject to boundary conditions on positions and velocities at a point. The approaches mentioned above do not seem easy to extend to the case when $q\ge 2$ or $M$ is not a Lie group with $k\ge 2$.

It seems that there are not many results on the case when $q\ge 2$ and $k\ge 2$ in the literature. In 2019, Heeren, Rumpf, and Wirth [13] proved the existence of weak solutions of Riemannian cubics in the Sobolev space $W^{2,2}$ through energy-minimizing sequences when $q\ge 2$, $k=2$, and $M$ is a Riemannian manifold. They also provided an example in <cit.> showing the nonexistence of Riemannian splines when $\lambda=0$.
Note that, as $\lambda=0$, the existence of classical Riemannian splines seems to be still an open problem.

§.§ Our results

To the best of the author's knowledge, it seems that there is no rigorous proof on the existence of classical solutions to the problem of higher-order spline interpolation problem in the literature.
We provide the existence of global solutions in the Hölder spaces for the gradient flow, whose asymptotic limits establish the existence of Riemannian splines, thereby offering solutions to spline interpolation problems on smooth Riemannian manifolds.

In this article, let $M$ be an $m$-dimensional smooth Riemannian manifold, and $T_{p}M$ be the tangent space of $M$ at $p\in M$ with the inner product $\langle \cdot , \cdot \rangle_{p}$ on $T_{p}M$.
Denote by $D$ the Levi-Civita connection associated with $M$.
Suppose $X$, $Y$, and $Z$ are smooth vector fields in $M$,
and define the so-called covariant derivative of $Y$ with respect to $X$ by $D_X Y$.
The Riemannian curvature tensor
$R(\boldsymbol{\cdot}, \boldsymbol{\cdot}) \boldsymbol{\cdot}: TM \times TM \times TM \rightarrow TM$ is defined by the identity, $R(X,Y) Z=D_X D_Y Z- D_Y D_X Z- D_{[X,Y]} Z$, where $[X,Y]$ denotes the Lie bracket of the vector fields $X$ and $Y$. Thanks to Nash's embedding theorem, we may assume that $M \subset\mathbb{R}^n$ is isometrically embedded for some $n>m$. Denote by $T_p M \subset T_p\mathbb{R}^n$ the $m$-dimensional tangent space to $M$ at a point $p \in M$.
We further assume that the Riemannian curvature tensor $R_M$ and the second fundamental form $\Pi_M$ of the immersed Riemannian manifold $M\subset\mathbb{R}^n$ are smooth and uniformly bounded in a $C^\infty$-norm.

Given a set of ordered (knot) points in a Lie group or a Riemannian manifold $M$,

$$\mathcal{P}=\{p_0, \cdots, p_q\}\subset M, $$

define the set of admissible maps by

\begin{align*}
\Phi_{\mathcal{P}}
\gamma:[x_0,x_q]\rightarrow M:
\gamma\in C^{2k-2}([x_0,x_q]), \gamma_l \in C^{2k,\alpha}([x_{l-1},x_l]), \gamma(x_l)=p_l, x_l=l, \forall\, l
\right\},
\end{align*}

$k \in \mathbb{N}\setminus \{1\}$, $I_{l}=(x_{l-1},x_{l})$, and $\gamma_l=\gamma_{|_{\bar I_l}}$,  for all $ l \in \{1,\ldots, q\}$.
Note that $\gamma \in C^m(\bar I)$ means
$\partial_{x}^\ell \gamma$ are continuous up to the boundary $\partial I$, $\forall\, \ell \in\{0,1,\cdots, m\}$,
where $I=(x_0, x_q)$.
The induced tangent map
$\partial_x \gamma: I \subset \mathbb{R}\rightarrow T_\gamma M$
and $\partial_x \gamma_l: I_l \subset \mathbb{R}\rightarrow T_\gamma M$ may be denoted as $\gamma_x$ and $\gamma_{l,x}$ respectively.
Define the functional
$\mathcal{E}_{k}^{\lambda}:  \Phi_{\mathcal{P}} \rightarrow [0,\infty)$ by

\begin{equation}
\label{eq:E_s-energy}
\mathcal{E}_{k}^\lambda [\gamma] =\mathcal{E}_{k} [\gamma]+\lambda\cdot \mathcal{T} [\gamma],
\end{equation}

where $\lambda\ge 0$ is a constant,

\begin{align}
\label{eq:E-energy}
\mathcal{E}_{k} [\gamma] =&
\sum\limits_{l=1}^{q} \frac{1}{2}
\int_{I_l} |D_{x}^{k-1} \gamma_{l,x}(x)|^2 \, dx, \\
\label{eq:T-energy}
\mathcal{T} [\gamma] =&
\sum\limits_{l=1}^{q} \frac{1}{2}
\int_{I_l} | \gamma_{l,x}(x)|^2 \, dx,
\end{align}

and $|\cdot|=\sqrt{\langle\cdot,\cdot\rangle}$ represents the norm of the fiber space in $TM$.

When $\gamma\in \Phi_{\mathcal{P}}$ is a critical point of the functional $\mathcal{E}_{k}^\lambda$, then,
by applying the first variation formulae (see Lemma <ref> in <ref> Appendix) and the boundary conditions,
we obtain

\begin{align*}
\frac{d}{d\varepsilon}_{|\varepsilon=0} \mathcal{E}_{k}^{\lambda}[\gamma]
= \sum_{l=1}^{q}\int_{I_l}\, \langle -\mathcal{L}_x^{2k}(\gamma_l), W\rangle \, dx =0
\text{,}
\end{align*}

for any suitable smooth vector field $W:I \rightarrow T_{\gamma}M$ along $\gamma_l(I)\subset M$,

\begin{align}
\label{Euler-eq-RPT}
\mathcal{L}_x^{2k}(\gamma_l)
:=(-1)^{k+1} D_{x}^{2k-1} \gamma_{l,x}
+ \sum^{k}_{\mu=2} (-1)^{\mu+k+1} R\left(D^{2k-\mu-1}_{x}  \gamma_{l,x}, D^{\mu-2}_{x}  \gamma_{l,x} \right) \gamma_{l,x} +\lambda \cdot D_{x} \gamma_{l,x}.
\end{align}

Thus, if $\gamma_l $ is the critical point, then

\begin{align}\label{1.1}
\mathcal{L}_x^{2k}(\gamma_l)=0
~~ \text{(intrinsic formulation)}
~~ \text{ or }~~
\mathcal{L}_x^{2k}(\gamma_l) \perp T_{\gamma_l} M
~~ \text{(extrinsic formulation)}.
\end{align}

For any fixed $l$, $\gamma_l$ is called the geometric $k$-spline or Riemannian polynomial spline. As $k=2$, it is called a cubic spline or Riemannian cubic. The critical points of $\mathcal{E}_k^\lambda$ could also be viewed as poly-harmonic maps (as $k\ge2$) or biharmonic maps (as $k=2$) with tension energy. Moreover, we may call all of these critical points the Riemannian splines in the rest of the article (and Riemannian cubics as $k=2$). Thus, one could view the so-called Riemannian spline interpolations as the problem of finding the poly-harmonic maps or biharmonic maps $\gamma:[a,b]\rightarrow M$ such that $\gamma(x_l)=p_l\in M$, $\forall\, l\in \{0, \ldots, q\}$, and $a=x_0<x_1<\cdots<x_q=b$. We show in <ref> that the spline interpolation problem on smooth Riemannian manifolds can be resolved by the gradient flow method. In Theorem <ref>, we prove the existence of global classical solutions to the gradient flows, whose asymptotic limits provide classical solutions to the interpolation problems on smooth Riemannian manifolds, called higher-order Riemannian splines or geometric $k$-splines. Moreover, we provide $C^{2k-2}$-smoothness for the geometric $k$-splines at the inner-boundary points, which was expected to be only $C^1$ at the inner-boundary points as $k=2$ in <cit.>.

Let $\lambda\in (0,\infty)$,
Assume that $M \subset\mathbb{R}^n$ is an $m$-dimensional Riemannian manifold, for some $n>m$,
and the second fundamental form $\Pi_M$ of $M$ is smooth and uniformly bounded with respect to a $C^\infty$-norm.
Suppose that $\gamma_0\in \Phi_{\mathcal{P}}$ is an initial datum
and fulfills the compatibility conditions of order $0$ in Definition <ref>.
Then, there exists a global solution $\gamma$ to
with the regularity

\begin{equation*}
\begin{cases}
\gamma_{l}
\in
\left([0,\infty)\times [x_{l-1}, x_{l}]\right) \bigcap C^{\infty}
\left((0,\infty)\times [x_{l-1}, x_{l}]\right),
~~~~~~~~\forall\, l\in\{1,\ldots,q\},\\
\gamma(t,\cdot) \in C^{2k-2}([x_{0}, x_{q}]), ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\forall\, t\in[0,\infty).
\end{cases}
\end{equation*}

Moreover, there exists  $\gamma_{\infty} \in \Phi_{\mathcal{P}}$
with the regularity $\gamma_{\infty} \in C^{2k-2,1}([x_0,x_q])$
such that
$\gamma_{\infty}=\underset{t_j\rightarrow\infty}\lim \gamma(t_j,\cdot)$,
for some convergent subsequences
and $\gamma_{\infty, l}$ is a Riemannian spline fulfilling (<ref>),
for each $l$.

§.§ The method and structure of article

The proof of the main theorems is based on Solonnikov's theory of parabolic partial differential systems and $L^2$-estimates of the curvature using inequalities of the Gagliardo-Nirenberg type.

The short-time existence of solutions is obtained by applying Solonnikov's theory to the linearized parabolic systems and employing a fixed-point argument in a suitable Banach space. Since the induced mathematical objects are connected curves or networks (represented as a collection of maps), we need to reformulate the problem by rearranging the domain (e.g., see (<ref>)) and initial-boundary conditions so that we can apply Solonnikov's theory to linear parabolic systems. Comparing with [15, 14], where the sufficient conditions for applying Solonnikov's theory are required to be verified, the complementary conditions at the boundaries in this article are much more complicated due to the higher order of the gradient flow.

To establish the long-time existence of solutions, it is necessary to derive uniform bounds for a specific sequence of higher-order $L^2$-norms, allowing for the application of a contradiction argument.
In the proof of Theorem <ref>, we aim to obtain the following estimate:

\begin{align*}
\sum_{l=1}^{q}\left( \frac{d}{dt}| D_{t}^{\mu} \gamma_{l}|^2_{L^2(I_{l})}
+| D_{t}^{\mu} \gamma_{l}|^2_{L^2(I_{l})}
+| D_{x}^{2k\mu+k-1}\partial_x \gamma_{l} |^2_{L^2(I_{l})} \right)
=\text{lower order terms}.
\end{align*}

Then we apply Gagliardo-Nirenberg interpolation inequalities and the Gronwall inequalities to derive a uniform bound for $| D_{x}^{2k\mu-1}\partial_x \gamma_{l} |^2_{L^2(I_{l})}$. The challenge lies in estimating the boundary terms generated from integration by parts in the higher-order energy estimates. The key is to apply the boundary conditions to reduce the order of derivatives, allowing for the application of Gagliardo-Nirenberg interpolation inequalities and subsequently obtaining the Gronwall inequalities.

The rest of this article is organized as following.
In <ref>, we give the existence of Riemannian splines to the interpolation problem on Riemannian manifolds
by gradient flow method
(see Theorem <ref>).
To keep the article self-contained, we collect some notations and technical results from the literature in
<ref> Appendix.

§ THE SPLINE INTERPOLATION PROBLEM

Let $\gamma:[0,T]\times I\rightarrow M$ fulfill $\gamma(t,\cdot)\in \Phi_{\mathcal{P}}$, for any fixed $t\in[0,T]$.
From the first variation formulae
(see Lemma
<ref> in <ref> Appendix),
the gradient flow of $\mathcal{E}_{k}^{\lambda}$ can be set up as,

\begin{equation}
\label{eq:E_s-flow}
\partial_t \gamma_{l}= \mathcal{L}_x^{2k}(\gamma_l), ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\forall ~ l\in\{1,\cdots,q\}
\text{,}
\end{equation}

with the initial-boundary conditions,

\begin{align}
&\text{(initial datum)}
\nonumber
\\
\gamma(0,x)=\gamma_{0}(x)\in \Phi_{\mathcal{P}},
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ x\in [x_0, x_q],
\\
&\text{(Dirichlet boundary conditions)}
\nonumber
\\
\gamma(t,x_{l})=p_{l}, ~~~~~~~~~~~~~~~~~~~~~~ l \in \{0,\ldots,q\}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ t\in[0,T],
\\
D_{x}^{\mu-1}\partial_x \gamma_l (t,x^{\ast})=v^{\mu}_{x^{\ast}} ~~~~~~~~~~  (l,x^{\ast})\in \{(1, x_0), (q, x_q)\}, ~ 1\leq \mu \leq k-1, ~ t\in[0,T],
\\
&\text{(regularity conditions at the knot points)}
\nonumber
\\
\Delta_{l} D_{x}^{\mu-1}\partial_x \gamma (t)=0,
~~~~~~~~~~~~~~ l \in \{1,\ldots,q-1\}, ~~~~~~~~ 1 \leq \mu \leq 2k-2, ~ t\in[0,T],
\end{align}

where the boundary datum $\{v^{\mu}_{x^{\ast}}\}$ in (<ref>) are constants and

\begin{align}\label{def:Delta_l}
\Delta_{l} D_{x}^{\mu-1}\partial_x \gamma (t)
:=D_{x}^{\mu-1}\partial_x \gamma_{l+1}(t,x_{l}^{+})
-D_{x}^{\mu-1}\partial_x \gamma_{l}(t,x_{l}^{-}).
\end{align}

The conditions, (<ref>) and (<ref>), aim to
keep $\gamma(t,\cdot)\in \Phi_{\mathcal{P}}$, for any fixed $t$.

From (<ref>), (<ref>),
and the boundary conditions
we obtain the energy identity of $\mathcal{E}^\lambda_{k}$, i.e.,

\begin{equation*}
\frac{d}{dt} \mathcal{E}^\lambda_{k}[\gamma]= - \int_{I} |\partial_t \gamma|^2 \,dt .
\end{equation*}

As $\lambda>0$, we obtain from the energy identity that

\begin{align}\label{ineq:f_x-higher-order}
\sum\limits^{q}_{l=1}\|  \gamma_{l,x} \|_{L^2(I_{l})}^2 (t) \le \frac{1}{\lambda} \mathcal{E}^\lambda_{k}[\gamma(0,\cdot)],
~~~~~~~~~~~~~~ \forall\, t\ge0.
\end{align}

[The geometric compatibility conditions of order $0$ to (<ref>)$\sim$(<ref>)]
We say that the initial datum $\gamma_0$ with the smoothness $\gamma_0\in C^{2k-2}([x_0, x_q])$ and $\gamma_{0, l} \in C^{2k,\alpha}([x_{l-1}, x_{l}])$,
$\forall\, l$,
fulfills the compatibility conditions of order $0$ to the parabolic system (<ref>)$\sim$(<ref>) if

\begin{equation}
\label{eq:cc-order-0}
\begin{cases}
\mathcal{L}_x^{2k}(\gamma_{0, l})(x^{\ast}) =0, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ x^{\ast} \in \{x_{l-1}^{+}, x_{l}^{-}\},
~~~~~ l \in \{1,\ldots,q\},
\\
\gamma_{0}(x_{l})=p_{l}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
l \in \{0,\ldots,q\},
\\
D_{x}^{\mu-1}\partial_x\gamma_{0, l} (x^{\ast})=v^{\mu}_{x^{\ast}} ~~~~~~~~~~~~~~~1\leq \mu \leq k-1, ~~ (l, x^{\ast}) \in \{(1, x_0), (q, x_q)\},
\\
\Delta_{l} D_{x}^{\mu-1}\partial_x\gamma_{0}=0, ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1 \leq \mu \leq 2k-2,  ~~~ l \in \{1,\ldots,q-1\},
\end{cases}
\end{equation}

In this section, our goal is to establish the existence of local solutions to equations
To achieve this, we will transform the system
(<ref>)$\sim$(<ref>) into a parabolic system describing a map into some Euclidean space, denoted as $g:[0,t_0)\times \bar I \rightarrow \mathbb{R}^{n\cdot q}$.
This reformulation facilitates the application of Solonnikov's theory and the Banach fixed-point argument, enabling us to obtain short-time solutions.

Suppose that $I\subset\mathbb{R}$ is an open interval,
$\gamma:I \rightarrow M\subset\mathbb{R}^n$ a smooth map,
and $V: I\rightarrow T_\gamma M\subset T_\gamma \mathbb{R}^n$
is a smooth vector field along $\gamma$.
Then, the (directional) covariant differentiation, denoted as $D_{\gamma_x}$ or simply as $D_x$, associated to the Levi-Civita connection $D$, can be written as

\begin{align}\label{eq:co-deri}
D_x V=
\partial_x V- \Pi_\gamma(V, \partial_x\gamma),
\end{align}

$\Pi_\gamma(\cdot,\cdot):T_\gamma M\times T_\gamma M\rightarrow(T_\gamma M)^\perp$
represents the second fundamental form at $\gamma\in M\subset\mathbb{R}^n$.
As $M$ is assumed to be isometrically embedded in $\mathbb{R}^n$ for some $n\in\mathbb{N}$, equation (<ref>) allows us to convert covariant differentiation $D_x$ into partial differentiation $\partial_x$.
From (<ref>), we can inductively deduce

\begin{align*}
D_x^i\gamma_{l,x} = \partial_x^{i+1}\gamma_{l} + W_i(\partial_x^{i}\gamma_{l}, \ldots, \partial_x\gamma_{l}, \gamma_{l}),
\end{align*}

$W_i:\mathbb{R}^n \times \cdots \times \mathbb{R}^n \times M \rightarrow \mathbb{R}^n$ fulfills

\begin{align}\label{est:W_i-2nd-fund-form}
\|W_i\|_{C^{k-i+1}(\mathbb{R}^n \times \cdots \times \mathbb{R}^n \times M)} \le C \cdot \|\Pi_M\|_{C^{k}(M)}
\text{,}
\end{align}

\|\Pi_M\|_{C^{k}(M)}:=\sum_{0\le |\sigma|\le k}\, \|\nabla^{\sigma} \Pi_M(\gamma_l)\|_{C^{0}(M)},
\nabla^{\sigma}:=\partial_{u_1}^{\sigma_1} \cdots \partial_{u_n}^{\sigma_n}
~~\text{ as } ~~
\sigma=(\sigma_1, \cdots, \sigma_n) .

Consequently, equation (<ref>) and
the boundary conditions can be written as partial differential equations after converting the covariant differentiation.
Thus, we rewrite the equations
(<ref>)$\sim$(<ref>) in an extrinsic formulation as,

\begin{align}
\label{eq:E_s-flow-extrinsic}
&\partial_t \gamma_{l}= (-1)^{k+1} \partial_x^{2k}\gamma_{l}
+ F_{l} (\partial_x^{2k-1}\gamma_{l}, \ldots, \partial_x\gamma_{l}, \gamma_{l})~~~~~~\text{ in } (0,T) \times I_{l}, ~ l \in \{1,\ldots,q\},
\\
\gamma_{l}(0,x)=\gamma_{0, l}(x), ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x \in I_{l}, ~ l \in \{1,\ldots,q\},
\\
\gamma_{l}(t,x_{l-j})=p_{l-j},
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  t \in [0,T], ~ j \in \{0,1\}, ~ l \in \{1,\ldots,q\},
\\
\partial_{x}^{\mu}\gamma_{l}(t,x^{\ast})=b^{\mu}_{l, x^{\ast}},
~~~~~~~~~ t \in [0,T], ~\mu \in \{1,\ldots, k-1\}, ~ (l, x^{\ast}) \in \{(1, x_0), (q, x_q)\},
\\
\partial_{x}^{\mu}\gamma_{l+1}(t, x^{+}_{l})-\partial_{x}^{\mu}\gamma_{l}(t, x^{-}_{l})=0,
t \in [0,T], ~ \mu \in \{1,\ldots, 2k-2\}, ~ l \in \{1,\ldots,q-1\},
\end{align}

\begin{align}
\label{eq:F_l}
F_l(\partial_x^{2k-1}\gamma_l, \ldots, \gamma_l)
&= (-1)^{k+1} W_{2k-1}(\partial_x^{2k-1}\gamma_l, \ldots, \gamma_l)
\\ \nonumber
&+ \sum^{k}_{\mu=2} (-1)^{\mu+k+1}
\partial_x^{2k-\mu}\gamma_l
+W_{2k-\mu-1}(\partial_x^{2k-\mu-1}\gamma_l, \ldots, \gamma_l)\right. ,
\\ \nonumber
& ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\left. \partial_x^{\mu-1}\gamma_l
+W_{\mu-2}(\partial_x^{\mu-2}\gamma_l, \ldots, \gamma_l)
\right)
\partial_x \gamma_l
\\ \nonumber
&+\lambda \cdot \left( \partial_x^2\gamma_l+ W_{1}(\partial_x\gamma_l, \gamma_l) \right)
\text{.}
\end{align}

§.§ The nonlinear parabolic systems for the splines

(60, 45)

(-45.5,38) *2
(36.5,-33.5) *2

(0,0) circle(1);
[dashed] (0,0) ellipse (1 and .2);
[black, xshift=1cm] plot [smooth, tension=2] coordinates  (-0.2,0) (-0.5,-0.4) (-1.3,0.4) (-1.8,-0.4);

The orientation of $g$ as $q=4$.

To tackle the nonlinear interpolation problem,
we transform the formulation of problem
(<ref>)$\sim$(<ref>) into a parabolic system with appropriate boundary conditions. Utilizing the assumption $x_l=l$ for all $l\in\mathbb{Z}$, we define the parametrization as follows,

\begin{align}
\label{eq:parametrization}
g_{l}(t,y) = \gamma_{l}\left(t, (-1)^{l-1}y + x_{2\left[\frac{l}{2}\right]}\right),
\quad g_{0,l}(y) = \gamma_{0,l}\left((-1)^{l-1}y + x_{2\left[\frac{l}{2}\right]}\right).
\end{align}

Here, $t\in [0,T]$, $y\in [0,1]$, and $[A]$ denotes the greatest integer less than or equal to $A$.

$\bullet$ The case when $q$ is even.

Note that the orientation of $g_{2l}$ will be reversed, and the boundary conditions will change sign if the order of derivatives is odd.
Specifically, the boundary conditions in
(<ref>) at $x_{l}$ become

\begin{align*}
\partial_{y}^{\mu}g_{2l-y^{\ast}}(t, y^{\ast})-(-1)^{\mu}\partial_{y}^{\mu} g_{2l-y^{\ast}+1}(t, y^{\ast})=0,
\end{align*}

where $(t,y^{\ast}) \in [0,T]\times \{0,1\}$,
$l \in \{1,\ldots, \frac{q}{2}-1+y^{\ast}\}$,
$\mu \in \{1, \ldots,2k-2\}$.
Together with a straightforward computation, the formulation of problem
can be converted into the parabolic system,

\begin{equation}\label{parabolic-equation}
\partial_t g_{l}= (-1)^{k+1}\cdot \partial_y^{2k} g_{l}
+G_l(\partial_y^{2k-1} g_{l}, \ldots, \partial_y g_{l}, g_{l})
~~~~ \text{ in } (0,T) \times (0,1), ~ l \in\{1,\cdots,q\},
\end{equation}

with the initial-boundary conditions

\begin{equation}\label{BC-q-even}
\begin{cases}
\text{(initial conditions)} \\
g_{l}(0,y)= g_{0,l}(y),  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   y \in [0,1], ~~ l\in \{1,\ldots, q\}, \\
\text{(boundary conditions)}
\\
g_{1}(t,0)=p_{0}, ~~~~ g_{q}(t,0)=p_{q}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ t \in [0,T],
\\
g_{2l-y^{\ast}}(t, y^{\ast})=g_{2l-y^{\ast}+1}(t, y^{\ast})=p_{2l-y^{\ast}},
~~ (t, y^{\ast}) \in [0,T]\times \{0,1\}, ~l \in \{1,\ldots, \frac{q}{2}-1+y^{\ast}\},
\\
\partial_{y}^{\mu}g_{l}(t,0)=b^{\mu}_{l, 0}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ t \in [0,T], ~ \mu \in \{1,\ldots, k-1\}, ~ l \in \{1, q\},
\\
\partial_{y}^{\mu}g_{2l-y^{\ast}}(t,y^{\ast})+(-1)^{\mu-1}\partial_{y}^{\mu} g_{2l-y^{\ast}+1}(t,y^{\ast})=0, \\
~~~~~~~~~~~~~~~~~~~~~~~~~~~(t, y^{\ast}) \in [0,T]\times \{0,1\}, ~ \mu \in \{1, \ldots,2k-2\}, ~ l \in \{1,\ldots, \frac{q}{2}-1+y^{\ast}\},
\end{cases}
\end{equation}

\begin{align}\label{eq:G_l-F_l}
G_l(\partial_y^{2k-1} g_{l}, \ldots, g_{l}):
=F_{l} \left((-1)^{(l-1)(2k-1)}\partial_y^{2k-1} g_{l}, \ldots, (-1)^{(l-1)\mu}\partial_y^{\mu} g_{l}, \ldots, g_{l}\right),
\end{align}

$b^{\mu}_{1, 0}= b^{\mu}_{1, x_0}$, $b^{\mu}_{q, 0}=(-1)^{(q-1)\mu} b^{\mu}_{q, x_q}$ are constants.

$\bullet$ The case when $q$ is odd.

As $q$ is odd, the formulation of problem
can be converted into (<ref>)
with the initial-boundary conditions

\begin{equation}
\label{BC-q-odd}
\begin{cases}
\text{(initial conditions)}
\\
g_{l}(0,y)= g_{0,l}(y),  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ y \in [0,1], ~~ l\in \{1,\ldots, q\},
\\
\text{(boundary conditions)}
\\
g_{1}(t,0)=p_{0}, ~~~~ g_{q}(t,1)=p_{q}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ t \in [0,T],
\\
g_{2l-y^{\ast}}(t, y^{\ast})=g_{2l-y^{\ast}+1}(t, y^{\ast})=p_{2l-y^{\ast}}, ~ (t, y^{\ast}) \in [0,T]\times \{0,1\}, ~l \in \{1,\ldots, \frac{q-1}{2}\},
\\
\partial_{y}^{\mu}g_{1}(t, 0)=b^{\mu}_{1, 0}, ~ \partial_{y}^{\mu}g_{q}(t, 1)=b^{\mu}_{q, 1},
~~~~~~~~~~~~~~~~~~~~~ t \in [0,T], ~~ \mu \in \{1,\ldots, k-1\}, \\
\partial_{y}^{\mu}g_{2l-y^{\ast}}(t,y^{\ast})+(-1)^{\mu-1}\partial_{y}^{\mu} g_{2l-y^{\ast}+1}(t,y^{\ast})=0,
\\
~~~~~~~~~~~~~~~~~~~~~~~~ (t, y^{\ast}) \in [0,T]\times \{0,1\}, ~~ \mu \in \{1, \ldots,2k-2\},
~~ l \in \{1,\ldots, \frac{q-1}{2}\},
\end{cases}
\end{equation}

$b^{\mu}_{1, 0}=b^{\mu}_{1, x_0}$, and $b^{\mu}_{q, 1}=(-1)^{(q-1)\mu}b^{\mu}_{1, x_q}$.

Introducing the boundary operators $\mathcal{B}$ with $\mathcal{B} g=\Phi$, where $\Phi$ represents the boundary datum, allows for a simplification of the boundary conditions in
(<ref>) and (<ref>).
Specifically, let $g=(g_{1}, \ldots, g_{q}) \in\mathbb{R}^{nq}$, where $g_{\nu}=(g^{1}_{\nu}, \ldots, g^{n}_{\nu})$,
$\forall\, \nu\in\{1,\cdots,q \}$,
and denote by $g^{t}$ the transpose of the vector $g$.
The boundary conditions in (<ref>) and (<ref>) can then be expressed as

\begin{align*}
\sum\limits_{\substack{1 \leq \ell \leq n \\1 \leq \nu \leq q}} \mathcal{B}_{i \, (\nu-1)n+\ell} \left( y, t, \partial_{y},  \partial_{t} \right) g_{\nu}^{\ell}(t, y^{\ast})=\Phi_{ i }(t, y^{\ast}) ,~~i \in \{1, \ldots, knq\}, ~~ y^{\ast} \in \{0,1\}
\end{align*}

or simply as,

\begin{align*}
\mathcal{B}(y^{\ast}, t, \partial_{y}, \partial_{t})
g^{t}(t, y^{\ast})=\Phi(t, y^{\ast}),
~~~y^{\ast} \in \{0,1\}.
\end{align*}

Here, the matrices
$\mathcal{B}(y^{\ast},t,\partial_{y},\partial_{t})$ are given by

\begin{align*}
\mathcal{B}(y^\ast=0, t, \partial_{y}, \partial_{t})=
\text{diag}(B_0, \overbrace{B_1, \ldots, B_1}^{\frac{q}{2}-1}, B_0),
\mathcal{B}(y^\ast=1, t, \partial_{y}, \partial_{t})=\text{diag}(\overbrace{B_1, \ldots, B_1}^{\frac{q}{2}}),
\end{align*}

when $q$ is even, or

\begin{align*}
\mathcal{B}(y^{\ast}=0, t, \partial_{y}, \partial_{t})=\text{diag}(B_0, \overbrace{B_1, \ldots, B_1}^{\frac{q-1}{2}}),
\mathcal{B}( y^{\ast}=1, t, \partial_{y}, \partial_{t})=\text{diag}(\overbrace{B_1, \ldots, B_1}^{\frac{q-1}{2}}, B_0),
\end{align*}

when $q$ is odd.

\begin{align*}
=&B_{0}(y^\ast, t,\partial_{y}, \partial_{t})
Id_{n \times n}\\
\partial_{y} \cdot Id_{n \times n} \\
\vdots  \\
\partial_{y} ^{k-1} \cdot Id_{n \times n}\\
\end{pmatrix}\text{,}\\
=&B_1(y^{\ast}, t, \partial_{y}, \partial_{t})
Id_{n \times n}&0\\
0&Id_{n \times n}\\
\partial_{y}\cdot Id_{n \times n}&  \partial_{y} \cdot Id_{n \times n}\\
\vdots & \vdots\\
\partial_{y}^{2k-3}\cdot Id_{n \times n}& \partial_{y}^{2k-3} \cdot Id_{n \times n}\\
\partial_{y}^{2k-2}\cdot Id_{n \times n}& -\partial_{y}^{2k-2} \cdot Id_{n \times n}\\
\end{pmatrix}.
\end{align*}

Hence, $\mathcal{B}(y^{\ast},t,\partial_{y},\partial_{t})$ are diagonal block matrices.

(60, 45)

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_0$;(-0.4,-0.2) rectangle(0.4,0.4);          0

0        $B_0$;(-0.4,-0.2) rectangle(0.4,0.4);

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_0$;(-0.4,-0.2) rectangle(0.4,0.4);          0

$B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);

0       $B_0$;(-0.4,-0.2) rectangle(0.4,0.4);

The matrix $\mathcal{B}(y^\ast=0, t, \partial_{y}, \partial_{t})$ when $q=2$ or $q=4$.

(60, 45)

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);           0

0        $B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);

The matrix $\mathcal{B}(y^\ast=1, t, \partial_{y}, \partial_{t})$ when $q=2$ or $q=4$.

(60, 45)

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_0$;(-0.4,-0.2) rectangle(0.4,0.4);

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_0$;(-0.4,-0.2) rectangle(0.4,0.4);       0

0    $B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_0$;(-0.4,-0.2) rectangle(0.4,0.4);          0

$B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);

0        $B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);

The matrix $\mathcal{B}(y^\ast=0, t, \partial_{y}, \partial_{t})$ when $q=1$, $q=3$, or $q=5$.

(60, 45)

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_0$;(-0.4,-0.2) rectangle(0.4,0.4);

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);        0

0     $B_0$;(-0.4,-0.2) rectangle(0.4,0.4);

(magic) [matrix of nodes,left delimiter=(,right delimiter=)]

$B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);           0

$B_1$;(-0.8,-0.2) rectangle(0.8,0.4);
[dashed] (0,0.4)–(0,0.3);    [dashed] (0,-0.2)–(0,-0.1);

0        $B_0$;(-0.33,-0.33) rectangle(0.37, 0.37);

The matrix $\mathcal{B}(y^\ast=1, t, \partial_{y}, \partial_{t})$ when $q=1$, $q=3$, or $q=5$.

We say that the initial datum $g_0$ with the smoothness  $g_{0, l} \in C^{2k,\alpha}([0, 1])$, $\forall\, l$,
fulfills the compatibility conditions of order $0$ to the nonlinear parabolic system (<ref>)
and (<ref>)
(or (<ref>) and (<ref>))
at $y^{\ast} \in \{0,1\}$

\begin{align*}
\partial_{t}^{p}
\sum\limits_{\substack{1 \leq \ell \leq n \\1 \leq \nu \leq q}} \mathcal{B}_{i \, (\nu-1)n+\ell}
\left( y, t, \partial_{y}, \partial_{t} \right) g_{\nu}^{\ell}(t, y^{\ast}) \big{| }_{ t=0}=\partial_{t}^{p} \Phi_{ i }(t, y^{\ast}) \big{| }_{t=0}
% i \in \{1, \ldots, knq\}, p\in\mathbb{N}_0,
\end{align*}

hold for
$i \in \{1, \ldots, knq\}$ and $p\in\{0,1\}$
such that after converting $\partial_t$ into $\partial_y$ by (<ref>) the order of $\partial_y$ in
$\partial_t^p\mathcal{B}_{i(\nu-1)n+\ell}g_\nu^\ell$ is less than or equal to $2k$.

The map $g$ is derived from $\gamma$ by reversing the parametrization of $\gamma_l$ for a given $l$. The compatibility conditions in Definition <ref> have been adjusted to align with the geometric compatibility conditions of order $0$ in Definition <ref>. Consequently, we straightforwardly conclude that the initial datum $g_0$ obtained from $\gamma_0$ satisfies the compatibility conditions of order $0$ in Definition <ref> if $\gamma_0$ adheres to the geometric compatibility conditions of order $0$ in Definition <ref>.

§.§ The linear problem

\begin{equation*}
Y:= \left(C^{\frac{2k+\alpha}{2k}, 2k+\alpha} \left([0,T] \times [0,1]\right)\right)^{q}
\end{equation*}

be the Banach space associated with the norm
$\norm{g}_{Y}=\sum\limits^{q}_{l=1}\norm{g_{l}}_{C^{\frac{2k+\alpha}{2k},2k+\alpha}([0,T] \times [0,1])}$, $\forall \, g\in Y$.
For any $g_{0}\in C^{2k,\alpha}([0,1])$ and $T>0$, define $X^{T}_{g_{0}}\subset Y$ as

\begin{align*}
X^{T}_{g_{0}}=&\big{ \{ }g=(g_1, \ldots, g_{q}): g_{l} : [0,T] \times [0,1] \to \mathbb{R}^{n},\\
&~~~~~~~~~~~g_{l} \in C^{\frac{2k+\alpha}{2k}, 2k+\alpha} \left([0,T] \times [0,1]\right),
g_{l}(0, \cdot)=g_{0, l}(\cdot), \, l \in \{1,\ldots, q\} \big{ \} }.
\end{align*}

Let the subset $B_K\subset X^{T}_{g_{0}}$ be defined by

$$B_K=\{g \in X^{T}_{g_{0}}: \norm{g}_{Y} \leq K\}. $$

Note that $B_K$ is a bounded, closed, and convex subset of the Banach space $Y$.

Below we linearize the nonlinear parabolic system and verify the sufficient conditions in <cit.> so that we can apply Solonnikov's theory to obtain the existence of solutions to the linear parabolic system.

Observe that the equations in (<ref>) and (<ref>) are already linear.
We linearize the parabolic system in (<ref>) as the following:

\begin{equation}\label{AP-higher-order-linear}
\partial_t g_{l}+ (-1)^{k}\cdot \partial_y^{2k} g_{l}=
G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \partial_y \bar{g}_{l}, \bar{g}_{l})
~~~~ \text{ in } (0,T) \times (0,1), ~ l \in\{1,\cdots,q\},
\end{equation}

where $\bar{g} \in X^{T}_{g_0}$.
Hence, (<ref>) and (<ref>) are the linearized parabolic system of
(<ref>) and
while (<ref>) and (<ref>) are the linearized parabolic system of (<ref>) and (<ref>).

Observe that the left-hand side of the $2k$-th order linear parabolic system in (<ref>)
can be presented as $\mathcal{L}(y, t, \partial_{y}, \partial_{t})g^{t}$,

\begin{equation*}
\mathcal{L}(y, t, \partial_{y}, \partial_{t})=\text{diag }(\mathcal{L}_{\ell \ell})^{nq}_{\ell=1},
\end{equation*}

$$\mathcal{L}_{\ell \ell}(y,t,\partial_{y},\partial_{t})=\partial_{t}+(-1)^{k}\partial_{y}^{2k},
~~~~~~ \ell=(l-1)n+j, j \in \{1,\ldots, n\}, l \in \{1,\ldots, q\}. $$

In <cit.>, $\mathcal{L}_{0}$
denotes the principal part of $\mathcal{L}$.
As $\mathcal{L}$ coincides with its principal part, we simplify the notation to work only with $\mathcal{L}$, avoiding the use of $\mathcal{L}_{0}$.
We associate these differential operators with polynomials of coefficients that depend on $(t,y)$.
We make the replacement $\partial_{y}$ by $\mathrm{i} \xi$, where $\xi \in \mathbb{R}$, $\mathrm{i}=\sqrt{-1}$, and $\partial_{t}$ by $p\in \mathbb{C}$.

\begin{align*}
\mathcal{L}_{\ell \ell}(y, t, \mathrm{i} \xi, p)=p+(-1)^{k}(\mathrm{i} \xi)^{2k}=p+\xi^{2k},
~~~ \ell=(l-1)n+j, j \in \{1,\ldots, n\}, l \in \{1,\ldots, q\}.
\end{align*}

\begin{align*}
L(y,t,\mathrm{i}\xi,p):=\det \mathcal{L}(y,t,\mathrm{i} \xi,p)
=\left(p+ \xi^{2k}\right)^{nq}\text{,}
\end{align*}

\begin{align*}
\mathcal{\hat{L}}(y,t,\mathrm{i}\xi,p)
:=L(y, t, \mathrm{i}\xi, p)\mathcal{L}^{-1}(y, t, \mathrm{i}\xi, p)
=\text{diag}(\mathcal{\hat{L}}_{\ell \ell})^{nq}_{\ell=1}
\text{,}
\end{align*}

\begin{align*}
\mathcal{\hat{L}}_{\ell \ell}=\mathcal{\hat{L}}_{\ell \ell}(y,t,\mathrm{i} \xi, p)
=(p+ \xi^{2k})^{nq-1}.
\end{align*}

$\bullet$ The smoothness of coefficients.

Given the smoothness assumption on $\bar{g}$ and applying the technical lemmas,
Lemmas <ref> and <ref>, we can ascertain that the coefficients in the linear system (<ref>) satisfy

\begin{equation*}
G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \bar{g}_{l})
\in C^{\frac{1+\alpha}{2k}, 1+\alpha}([0,T] \times [0,1]) \subset C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1]),
\end{equation*}

for all $l \in \{1, \ldots, q\}$.

$\bullet$ The parabolicity condition.

From its definition,
the roots of the polynomial $L(x, t, \mathrm{i}\xi, p)$ with respect to the variable $p$ (considering multiples of $n$) satisfy the condition $p=-\xi^{2k}$, where $\xi$ is a real number. Consequently, uniform parabolicity is established (refer to <cit.>).

$\bullet$ The complementary conditions of the initial datum.

Following the definition in <cit.>, the matrix associated with the initial condition takes the form

\begin{align*}
\mathcal{C}_0(y,\partial_{y},\partial_{t})=Id_{nq \times nq}.
\end{align*}

From the definition of
$\mathcal{C}_0(y, 0, p)$ and
we obtain

\begin{align*}
\mathcal{C}_0(y, 0, p)\mathcal{\hat{L}}(y, 0, 0, p)=\text{diag}(p^{nq-1}),
\end{align*}

demonstrating that the rows of the matrix are linearly independent modulo $p^{nq}$.
Hence, the complementary conditions of the initial datum $g_0$ are satisfied.

$\bullet$ The compatibility conditions of the initial datum at the boundaries.

Given that the boundary operators for $g$ in either (<ref>) or (<ref>) are linear, and considering $g_{l}(0,\cdot)=\bar{g}_{l}(0,\cdot)=g_{0,l}$ for all $l$, we can deduce that $g_0$ satisfies the compatibility conditions of order $0$ for the linear equation (refer to <cit.>). This conclusion holds if $g_0$ also satisfies the compatibility conditions of order $0$ for the nonlinear parabolic system, as defined in Definition <ref>. Therefore, the compatibility conditions of order $0$ for the linear parabolic systems are confirmed.

$\bullet$ The polynomial $M^{+}$ defined in <cit.>.

As $\text{Re} \text{ } p \geq 0$ and $p \ne 0$
in the polynomial $L=L(y,t,\mathrm{i}\xi,p)$,
the polynomial $L$ of $\xi$ has $k$ roots with positive real parts, each with multiplicity $nq$, and $k$ roots with negative real parts and multiplicities $nq$ (see <cit.>).
Based on the assumption on $p$, we can write $p=|p|e^{\mathrm{i} \theta_{p}}$ with $ -\frac{1}{2}\pi \leq \theta_{p} \leq \frac{1}{2}\pi$, $|p| \ne 0$.

\begin{align*}
\xi_{\mu}(y^{\ast}, p)=\sqrt[2k]{|p|}e^{\mathrm{i} \left(\frac{\theta_{p}}{2k}+\frac{(2\mu-1)\pi}{2k}\right)},
~~~~~~~~~~~~~~ \mu \in \{1, \ldots, 2k\}.
\end{align*}

$\left\{\xi_{\mu}(y^{\ast}, p): \mu=1, \ldots, k \right\}$ is the set of roots of $p+\xi^{2k}=0$ with positive imaginary parts, and
$\left\{\xi_{\mu}(y^{\ast}, p): \mu=k+1, \ldots, 2k \right\}$ is the set of roots of $p+\xi^{2k}=0$ with negative imaginary parts.

\begin{align*}
p+\xi^{2k}=\prod\limits_{\mu=1}^{2k}\left(\xi-\xi_{\mu}(y^{\ast}, p)\right)
\text{.}
\end{align*}

Since each root has multiplicity $nq$, we define

\begin{align*}
M^{+}(y^{\ast}, \xi, p)=\prod\limits_{\mu=1}^{k} \left(\xi-\xi_{\mu}(y^{\ast}, p)\right) ^{nq}.
\end{align*}

$\bullet$ The complementary conditions along the boundaries $ [0,T]\times \partial{I}$.

As the matrices $\mathcal{B}(y^{\ast},t,\partial_{y},\partial_{t})$ are composed of block matrices, it suffices to verify the complementary conditions for each block matrix at the boundaries, whether $q$ is even or odd and $y^{\ast} \in \{0,1\}$.

From the definition of
and letting $\partial_y^k=(\mathrm{i} \xi)^k$,  $\partial_t^\ell=p^\ell$,

\begin{align*}
\mathcal{B}(y^\ast=0, t, \mathrm{i} \xi, p)=\text{diag}(B_0, \overbrace{B_1, \ldots, B_1}^{\frac{q}{2}-1}, B_0),
\mathcal{B}(y^\ast=1, t, \mathrm{i} \xi, p)=\text{diag}(\overbrace{B_1, \ldots, B_1}^{\frac{q}{2}}),
\end{align*}

are diagonal block matrices.

Notice that $B_0(y^\ast=0,t,\partial_y,\partial_t)=B_0(y^\ast=1,t,\partial_y,\partial_t)$ and
$B_1(y^\ast=0,t,\partial_y,\partial_t)= B_1(y^\ast=1,t,\partial_y,\partial_t)$.

To confirm the complementary conditions as outlined in Solonnikov's theory (see <cit.>), we must demonstrate that the rows of the matrix

\begin{align*}
\mathcal{A}(y^{\ast}, t, \mathrm{i}\xi, p)=\mathcal{B}(y^{\ast}, t, \mathrm{i} \xi, p)\hat{\mathcal{L}}(y^{\ast}, t, \mathrm{i} \xi, p)
\end{align*}

are linearly independent modulo $M^{+}(y^{\ast}, \xi, p)$, where $\text{Re}\{p\} \geq 0$, $p \ne 0$.

Note that $\mathcal{A}$ is a diagonal block matrix,
since $\mathcal{B}$ is a diagonal block matrix and $\hat{\mathcal{L}}$ is a diagonal matrix.
Thus, to obtain the linear independence for the rows of the matrix $\mathcal{A}$, it is sufficient to verify the linear independence on each block matrix separately.
Additionally, the block matrices $B_0(y^\ast,t,\mathrm{i}\xi,p), B_1(y^\ast,t,\mathrm{i}\xi,p)$,
which are types of block matrices in $\mathcal{B}(y^\ast,t,\mathrm{i}\xi,p)$,
are the same regardless of whether $y^\ast$ is 0 or 1.
Hence, we only need to separately verify the complementary conditions at $y^{\ast}=0$ for
$\mathcal{A}_0 :=B_0 \cdot \mathcal{\hat{L}}_{11} Id_{n \times n}$
and $\mathcal{A}_1 :=B_1 \cdot \mathcal{\hat{L}}_{11} Id_{2n \times 2n}$.

Observe that

\begin{align*}
\mathcal{A}_0
=B_{0}(y^\ast=0, t, \mathrm{i}\xi ,p) \cdot \mathcal{\hat{L}}_{11}(y^\ast=0,t,\mathrm{i} \xi , p)
Id_{n \times n}
\mathcal{\hat{L}}_{11} \cdot Id_{n\times n}\\
\mathcal{\hat{L}}_{11} \cdot \mathrm{i} \xi \cdot Id_{n \times n}\\
\mathcal{\hat{L}}_{11} \cdot (\mathrm{i} \xi)^{2} \cdot Id_{n \times n}\\
\vdots \\
\mathcal{\hat{L}}_{11} \cdot (\mathrm{i} \xi)^{k-1} \cdot Id_{n \times n}\\
\end{pmatrix}.
\end{align*}

To verify the linear independence of the rows of matrix $\mathcal{A}_0$ modulo $M^{+}(y^\ast=0, \xi, p)$, we need to show that

\begin{align}
\label{comple-xi-higher-order}
\omega \mathcal{A}_0(y^\ast=0, t, \mathrm{i}\xi, p)=0
~\text{  mod }  M^{+}(y^\ast=0, \xi, p)
~~ \Longrightarrow ~~ \omega=(\omega_{1}, \ldots, \omega_{kn})=0.
\end{align}

We can rewrite (<ref>) as,

\begin{align*}
\left(\sum\limits^{k-1}_{\mu=0}(\mathrm{i}  \xi)^{\mu} \cdot \omega_{\mu n+j}\right) (p+ \xi^{2k})^{nq-1} =0
~ \text{ mod }  M^{+}(y^\ast=0, \xi, p),
~~ \forall\, j \in \{1,\ldots, n\},
\end{align*}

which is also equivalent to

\begin{align}
\label{eq-first-m-column}
\omega_{j}+\mathrm{i}  \xi  \omega_{n+j}+(\mathrm{i}  \xi)^{2}  \omega_{2n+j}
+\cdots +(\mathrm{i}  \xi)^{k-1}  \omega_{(k-1)n+j}=0
~ \text{ mod } \prod\limits^{k}_{\mu=1}\left(\xi-\xi_{\mu}(y^\ast=0,p)\right), ~ \forall\, j
\text{,}
\end{align}

by using the same argument in <cit.> on simplification.
By letting $\xi=\xi_{\mu}$ in (<ref>),
we obtain

\begin{align}
\label{system-first-m-column}
C \cdot
\bar{\omega}_{j}^{t}
\end{align}

(\omega_{j}, \omega_{n+j}, \omega_{2n+j},\cdots, \omega_{(k-1)n+j})$

\begin{align*}
C= \begin{pmatrix}
1&\mathrm{i} \xi_{1}&(\mathrm{i} \xi_{1})^{2}&\cdots&(\mathrm{i} \xi_{1})^{k-1}\\
1&\mathrm{i} \xi_{2}&(\mathrm{i} \xi_{2})^{2}&\cdots&(\mathrm{i} \xi_{2})^{k-1}\\
1&\mathrm{i} \xi_{3}&(\mathrm{i} \xi_{3})^{2}&\cdots&(\mathrm{i} \xi_{3})^{k-1}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\mathrm{i} \xi_{k}&(\mathrm{i} \xi_{k})^{2}&\cdots&(\mathrm{i} \xi_{k})^{k-1}\\
\end{pmatrix}.
\end{align*}

Notice that the matrix $C$ is the so-called Vandermonde matrix with the property,

\begin{align*}
\det C
=\prod\limits_{1 \leq \mu< \nu \leq k}(\mathrm{i} \xi_{\nu}-\mathrm{i} \xi_{\mu}) \ne 0.
\end{align*}

This implies that any solution to
(<ref>) is trivial, i.e., $\bar{\omega}_{j}=0$, $\forall\, j \in \{1,\ldots,n\}$.
Thus, (<ref>) holds.

Next, observe that

\begin{align*}
\mathcal{A}_1
&=B_{1}(y^\ast=0, t, \mathrm{i}\xi ,p) \cdot \mathcal{\hat{L}}_{11}(y^\ast=0,t,\mathrm{i} \xi , p)
Id_{2n \times 2n} \\
\begin{pmatrix}
\mathcal{\hat{L}}_{11} Id_{n \times n}& 0 \\
0& \mathcal{\hat{L}}_{11} Id_{n \times n}\\
\mathrm{i} \xi \cdot \mathcal{\hat{L}}_{11}Id_{n \times n}& \mathrm{i} \xi \cdot \mathcal{\hat{L}}_{11}Id_{n \times n}\\
(\mathrm{i} \xi)^{2}\cdot \mathcal{\hat{L}}_{11} Id_{n \times n}& -(\mathrm{i} \xi)^{2} \cdot \mathcal{\hat{L}}_{11} Id_{n \times n}\\
\vdots  &\vdots \\
(\mathrm{i} \xi)^{2k-3}\cdot \mathcal{\hat{L}}_{11} Id_{n \times n}&  (\mathrm{i} \xi)^{2k-3} \cdot \mathcal{\hat{L}}_{11} Id_{n \times n}\\
(\mathrm{i} \xi)^{2k-2}\cdot \mathcal{\hat{L}}_{11} Id_{n \times n}& -(\mathrm{i} \xi)^{2k-2} \cdot \mathcal{\hat{L}}_{11} Id_{n \times n}\\
\end{pmatrix}.
\end{align*}

To verify the linear independence of the rows of matrix $\mathcal{A}_1$ modulo $M^{+}(y^\ast=0, \xi, p)$, we need to show that

\begin{align}\label{comple-xi-y-1}
\omega \cdot \mathcal{A}_1 =0
~ \text{ mod }  M^{+}(y^\ast=0, \xi, p)
~ \Longrightarrow ~ \omega=(\omega_{1}, \ldots,\omega_{2kn})=0.
\end{align}

From the same argument as in <cit.>, for simplicity,
(<ref>) is equivalent to

\begin{equation*}
\begin{cases}
\omega_{j}+\sum\limits_{\mu=1}^{2k-2} (\mathrm{i} \xi )^{\mu}\cdot \omega_{(\mu+1)n+j} =0
~~~~~~~~~~~~~~ \text{ mod} \prod\limits^{k}_{\mu=1}\left(\xi-\xi_{\mu}(y^\ast=0,p)\right),
\\
\omega_{n+j}+\sum\limits_{\mu=1}^{2k-2} (-1)^{\mu-1} (\mathrm{i} \xi )^{\mu}\cdot  \omega_{(\mu+1)n+j} =0
~ \text{ mod } \prod\limits^{k}_{\mu=1}\left(\xi-\xi_{\mu}(y^\ast=0,p)\right),
\end{cases}
\Rightarrow
\begin{cases}
\omega=0, \\
\forall\, j\in\{1, \ldots,n\}.
\end{cases}
\end{equation*}

Thus, by substituting $\xi=\xi_{r}$, $ r\in \{1, \ldots, k\}$,
we obtain

\begin{equation*}
\begin{cases}
\omega_{j}+\sum\limits_{\mu=1}^{2k-2} (\mathrm{i} \xi_{r})^{\mu} \cdot \omega_{(\mu+1)n+j} =0, \\
\omega_{n+j}+\sum\limits_{\mu=1}^{2k-2} (-1)^{\mu-1} (\mathrm{i} \xi_{r})^{\mu} \cdot  \omega_{(\mu+1)n+j} =0 ,
\end{cases}
~~ \forall\, j \in \{1,\ldots, n\}.
\end{equation*}

By subtraction and addition between the first and the second equations, we obtain

\begin{equation*}
\begin{cases}
\omega_{j}-\omega_{n+j}+2\sum\limits_{\mu=1}^{k-1} (\mathrm{i}\cdot  \xi_{r})^{2\mu} \cdot \omega_{(2\mu+1)n+j} =0, \\
\omega_{j}+\omega_{n+j}+2\sum\limits_{\mu=1}^{k-1} (\mathrm{i}\cdot  \xi_{r})^{2\mu-1} \cdot \omega_{2\mu n+j} =0,
\end{cases}
\forall\, r\in \{1, \ldots, k\}, \forall\, j \in \{1,\ldots, n\},
\end{equation*}

which can be rewritten as

\begin{equation}\label{D+E-2}
\begin{cases}
D \cdot \bar{\omega}^{t}=0 , \\
E \cdot \tilde{\omega}^{t}=0
\text{.}
\end{cases}
\end{equation}

Here, we let

\begin{equation*}
\begin{cases}
\bar{\omega} =\begin{pmatrix}
\frac{1}{2}\left(\omega_{j}-\omega_{n+j}\right), &\mathrm{i}^{2}\omega_{3n+j}, &\mathrm{i}^{4}\omega_{5n+j}, &\cdots, &\mathrm{i}^{2k-2}\omega_{(2k-1)n+j}\\
\end{pmatrix}, \\
\tilde{\omega} =\begin{pmatrix}
\frac{1}{2}\left(\omega_{j}+\omega_{n+j}\right), &\mathrm{i}\omega_{2n+j}, &\mathrm{i}^{3}\omega_{4n+j}, &\cdots,  &\mathrm{i}^{2k-3} \omega_{(2k-2)n+j}\\
\end{pmatrix},
\end{cases}
\end{equation*}

\begin{align*}
\vdots&\vdots&\vdots&\ddots&\vdots\\
\end{pmatrix},
\vdots&\vdots&\vdots&\ddots&\vdots\\
\end{pmatrix}.
\end{align*}

We want to show that the solutions to (<ref>) are trivial, i.e,
$\bar{\omega}=\tilde{\omega}=0$, or equivalently to say,
$\omega_j=0$, $\forall\, j \in \{1, \ldots, 2kn\}$.
To confirm it, we only need to show that $\det{D} \ne 0, \det{E} \ne 0$.
Notice that the matrix $D$ is the so-called Vandermonde matrix with the property,

\begin{align*}
\det{D}=\prod_{1 \leq \mu <\nu \leq k}\left(\xi_{\nu}^2-\xi_{\mu}^2\right)
\text{.}
\end{align*}

$\xi_{\mu} \ne \xi_{\nu}$, as $\mu\ne \nu$, we conclude that $\det D\ne 0$.

On the other hand, the matrix $E$ is a Vandermonde-like matrix.
We want to show that $\det E\ne 0$ below.

\begin{align*}
\xi_\mu=c(p) \cdot e^{\mathrm{i} \frac{\mu\pi}{k}},
c(p)=|p|^{\frac{1}{2k}} e^{\mathrm{i} \left(\frac{\theta_{p}}{2k}-\frac{\pi}{2k}\right)}.
\end{align*}

Let $F_{\mu}$ be the matrix obtained from $E$
by taking off the first column and the $\mu$-th row.
Then, for fixed $\mu$,

\begin{align*}
\det{F_{\mu}}=&\det\begin{pmatrix}
\xi_{1}&\xi_{1}^3&\cdots&\xi_{1}^{2k-5} &\xi_{1}^{2k-3}\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
\xi_{\mu-1}&\xi_{\mu-1}^3&\cdots &\xi_{\mu-1}^{2k-5}&\xi_{\mu-1}^{2k-3}\\
\xi_{\mu+1}&\xi_{\mu+1}^3&\cdots&\xi_{\mu+1}^{2k-5}& \xi_{\mu+1}^{2k-3}\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
\xi_{k}&\xi_{k}^3&\cdots&\xi_{k}^{2k-5}&\xi_{k}^{2k-3}\\
\end{pmatrix}
\underset{\overset{1 \leq \nu  \leq k}{\nu \ne \mu}}{\prod}   ~~\xi_{\nu} \cdot ~~
\underset{\overset{1 \leq r<s \leq k}{r, s \ne \mu}}{\prod} \left(\xi_{s}^2-\xi_{r}^2\right).
\end{align*}

\begin{align*}
=\sum_{\mu=1}^{k} (-1)^{\mu-1}  \det{F_{\mu}}
=\sum_{\mu=1}^{k} (-1)^{\mu-1}
\underset{\overset{1 \leq \nu  \leq k}{\nu \ne \mu}}{\prod} ~~\xi_{\nu}
\cdot
\underset{\overset{1 \leq r<s \leq k}{r, s \ne \mu}}{\prod} \left(\xi_{s}^2-\xi_{r}^2\right) \\
\nonumber
= &\left(\prod_{1 \leq \nu \leq k}\xi_{\nu} \right)
\left(\underset{1 \leq r <s \leq k}{\prod} (\xi_{s}^2-\xi_{r}^2) \right)
\left(\sum_{\mu=1}^{k} \frac{(-1)^{\mu-1} }{ \xi_{\mu} } \cdot
\frac{1}{\underset{1 \leq \nu< \mu }{\prod} (\xi_{\mu}^2-\xi_{\nu}^2)}  \cdot
\frac{1}{\underset{ \mu< \nu \leq k}{\prod} (\xi_{\nu}^2-\xi_{\mu}^2)}   \right).
\end{align*}

For any $r, s \in \{1, \ldots, k\}$,
we have

\begin{align}\label{eq:xi-nu^2-xi-mu^2}
\xi_{r}^2-\xi_{s}^2 =&|p|^{\frac{1}{k}} e^{\mathrm{i} \left(\frac{\theta_{p}}{k}-\frac{\pi}{k}\right)}
\left(e^{\frac{ \mathrm{i} 2r \pi}{k}}-e^{\frac{ \mathrm{i} 2 s \pi}{k}}\right)
=2 \mathrm{i}\cdot (c(p))^2 \sin{\frac{(r-s) \pi}{k}} \cdot e^{\frac{ \mathrm{i} (r+s) \pi}{k}}.
\end{align}

Hence, to prove $\det E\ne 0$, we only need to show

\begin{align*}
H:= \sum_{\mu=1}^{k} \frac{(-1)^{\mu-1} }{ \xi_{\mu} } \cdot
\frac{1}{\underset{1 \leq \nu< \mu }{\prod} (\xi_{\mu}^2-\xi_{\nu}^2)}
\cdot \frac{1}{\underset{ \mu< \nu \leq k}{\prod} (\xi_{\nu}^2-\xi_{\mu}^2)}
\ne 0
\text{.}
\end{align*}

From applying
(<ref>) and the definition of $\xi_{\mu}$,
we obtain

\begin{align*}
\nonumber
H=& \sum_{\mu=1}^{k} \frac{(-1)^{\mu-1} }{c(p)e^{\frac{ \mathrm{i}  \mu \pi}{k}}}
\cdot \frac{1}{ \underset{1 \leq \nu < \mu \leq k}{\prod}2\mathrm{i} c^2(p)  \sin{\frac{ (\mu-\nu) \pi}{k}} e^{\frac{ \mathrm{i} (\nu+\mu) \pi}{k}}
\cdot \underset{1 \leq \mu < \nu \leq k}{\prod}2\mathrm{i} \cdot (c(p))^2  \sin{\frac{ (\nu-\mu) \pi}{k}} e^{\frac{ \mathrm{i} (\nu+\mu) \pi}{k}}}
\\
=& \sum_{\mu=1}^{k} A_\mu \cdot C_{\mu} \cdot K_\mu,
\end{align*}

\begin{align*}
\frac{(-1)^{\mu-1} }{ c(p)\cdot [ 2\mathrm{i} \cdot (c(p))^2 ]^{k-1}},
\frac{1}{\underset{1 \leq \nu < \mu}{\prod} \sin{\frac{ (\mu-\nu)\cdot \pi}{k}}}
\frac{1}{\underset{\mu<\nu\leq k}{\prod} \sin{\frac{ (\nu-\mu)\cdot \pi}{k}}},
\end{align*}

\begin{align*}
e^{-\frac{ \mathrm{i}  \mu \pi}{k}} \cdot \underset{\overset{1 \leq \nu  \leq k}{\nu \ne \mu}}{\prod} e^{-\frac{\mathrm{i} (\nu+\mu) \pi}{k}}
= e^{-\frac{ \mathrm{i}  \mu \pi}{k}} \cdot \underset{\overset{1 \leq \nu  \leq k}{\nu \ne \mu}}{\prod} e^{-\frac{\mathrm{i} \nu \pi}{k}}
\cdot \underset{\overset{1 \leq \nu  \leq k}{\nu \ne \mu}}{\prod} e^{-\frac{\mathrm{i} \mu \pi}{k}}
=(-\mathrm{i})^{k+1} \cdot (-1)^{\mu} \cdot e^{\frac{\mathrm{i} \mu \pi}{k}}.
\end{align*}

A simple calculation shows that

\begin{align*}
H=\frac{(-1)^{k+1}}{2^{k-1} (c(p))^{2k-1} } \sum_{\mu=1}^{k}   C_{\mu} \cdot e^{\frac{\mathrm{i} \mu \pi}{k}}.
\end{align*}

Observe that the imaginary part of
$e^{\frac{ \mathrm{i} \cdot \mu \pi}{k} }$ is non-negative
for each $\mu\in\{1,\cdots,k\}$.
From its definition, it is straightforward to verify that $C_{\mu}>0$, which leads to

\begin{align*}
\sum_{\mu=1}^{k} C_{\mu}\cdot e^{\frac{ \mathrm{i} \cdot \mu \pi}{k} } \ne 0
\text{,}
\end{align*}

Thus, $H\ne 0$, and therefore $\det E\ne 0$.

Now, we conclude that any solution to (<ref>) is trivial and hence (<ref>) is verified.
The verification of the complementary conditions at the boundary is finished.

Let $\lambda\in (0,\infty)$, $\alpha\in(0,1)$, and $\{x\}=x-[x]$.
Suppose that  $g_0=(g_{0,1}, \ldots, g_{0,q})$, $g_{0, l} : [0,1] \to \mathbb{R}^n$,
is an initial datum to the linear parabolic system (<ref>) with the initial-boundary conditions (<ref>) as $q$ is even, or the initial-boundary conditions (<ref>) as $q$ is odd.
Assume that
$g_{0, l}\in C^{2k,\alpha}([0,1]), \forall\, l \in \{1,\ldots, q\}$
and fulfills the compatibility conditions of order $0$ to the linear parabolic system with the initial-boundary conditions.
Then, for any $T>0$ and $\bar{g} \in X^{T}_{g_{0}}$,
there is a unique solution $g=(g_{1}, \ldots, g_{q} )$ to the linear parabolic system
(<ref>) with the initial-boundary conditions
such that $g_{l} \in C^{\frac{2k+\alpha}{2k}, 2k+\alpha}([0,T] \times [0,1])$, $\forall\, l \in \{1,\ldots, q\}$, and the induced map $\gamma$,
defined in (<ref>),
fulfills $\gamma(t,\cdot) \in C^{2k-2}([x_0,x_q])$, $\forall\, t \in [0, T]$.

\begin{align*}
\nonumber
\norm{g}_{X^{T}_{g_0}}
\leq & C_{0}\cdot \bigg{(} \sum^{q}_{l=1}\norm{  G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \bar{g}_{l})
}_{C^{\frac{\alpha}{2k},\alpha}([0,T] \times [0,1])} \\
&~~~~~~+\sum^{q}_{l=1}(| p_{l-1}|+| p_{l}|)
+\sum\limits_{\mu=1}^{k-1} \left( |b^{\mu}_{1,0}|+ |b^{\mu}_{q, 2 \{ \frac{q}{2} \}}|   \right)
+\sum^{q}_{l=1}\Vert g_{0, l} \Vert_{C^{2k,\alpha}([0,1])}\bigg{ )},
\end{align*}

where $C_0=C_0(n)$.

The proof is a consequence of verifying all the sufficient conditions outlined in <cit.>.

§.§ The proof of Theorem <ref>

§.§.§ The local solutions

The proof is obtained from applying Theorem <ref> and Banach fixed point Theorem.

[The short-time existence]

Under the assumptions in Theorem <ref>, we further assume that
$g_0=(g_{0,1}, \ldots, g_{0,q}):[0,1]\rightarrow \mathbb{R}^{nq}$,
fulfilling $g_{0, l}\in C^{2k,\alpha} \left([0,1]\right)$,
satisfies the compatibility conditions of order $0$ in Definition <ref> and is an initial datum to the nonlinear parabolic system (<ref>) with the initial-boundary conditions
(<ref>) as $q$ is even, or with the initial-boundary conditions  (<ref>) as $q$ is odd.

Then, there exists

$$t_0=t_0 \left(m, n, \alpha, \lambda, \sum\limits^{q}_{l=1} \Vert g_{0, l}\Vert_{C^{2k,\alpha}([0,1])}, K_0, \|R_M\|_{C^{1}(M)} , \|\Pi_M\|_{C^{2k-1}(M)}\right)>0,

such that $g=(g_1,\ldots, g_{q})$ fulfills the regularity,

\begin{equation*}
g_l  \in C^{\frac{2k+\alpha}{2k},2k+\alpha}\left([0,t_0] \times [0,1]\right)\bigcap
C^{\infty}\left((0,t_0]\times [0,1]\right), ~~~~~~l\in\{1,\ldots, q\},
\end{equation*}

and $g$ is the unique solution  to the nonlinear parabolic system (<ref>) with the initial-boundary conditions
(<ref>) or (<ref>),
$\|R_M\|_{C^{1}(M)}$ is the $C^{1}$-norm of the Riemannian curvature tensor of manifold $M$,
$\|\Pi_M\|_{C^{2k-1}(M)}$ is the $C^{2k-1}$-norm of the second fundamental form $\Pi_M$ of the embedded $M$ in $\mathbb{R}^n$,

\begin{align*}
K_0&= 2C_0 \cdot \bigg{(} \sum^{q}_{l=1}
\Vert G_l(\partial_y^{2k-1}g_{0, l}, \ldots, g_{0, l})  \Vert_{C^{0,\alpha}([0,1])}
+\sum\limits_{l=1}^{q} (| p_{l-1}|+|p_{l}|)
\\
+\sum^{q}_{l=1}\Vert g_{0,l} \Vert_{C^{2k,\alpha}([0,1])} \bigg),
\end{align*}

for some constant
$C_0=C_0(m, n)>0$.

Moreover, the induced map $\gamma: [x_0, x_q]\rightarrow M \subset \mathbb{R}^n$, defined in (<ref>), fulfills
$\gamma(t,\cdot) \in C^{2k-2}([x_0, x_q]), ~ \forall\, t \in [0,t_0].$

Applying Theorem <ref> to (<ref>) with the initial-boundary conditions (<ref>) or (<ref>), we introduce the operator

\begin{align*}
\nonumber
\mathcal{G} :  X^{T}_{g_{0}} \cap B_K & \to  X^{T}_{g_{0}} \cap B_K\\
\bar{g}  &\mapsto  g
\text{,}
\end{align*}

where $g=(g_1, \ldots, g_{q})$ is a solution to (<ref>) with the initial-boundary conditions for any sufficiently small $T>0$.

$\bullet$ Self-maps. i.e., $\exists$ $t_2 >0$ so that
$\mathcal{G}\left( X^{T}_{g_{0}} \cap B_{K_0}\right) \subset X^{T}_{g_{0}} \cap B_{K_0}$,
$\forall\, T\in(0, t_2)$.

We will show the self-map property on $B_{K_0}$ by choosing a proper $K_0>0$ below.
Below we use the notation

\begin{align*}
=(R\circ\gamma) (\boldsymbol{\cdot},\boldsymbol{\cdot})\boldsymbol{\cdot}
\end{align*}

to represent the Riemannian curvature tensors of $M$.
By applying
the triangle inequalities in the
Hölder spaces, (<ref>), and
the technical lemmas,
Lemmas <ref>$\sim$<ref>,
with $v=\bar{g}_l-g_{0,l}$, we obtain

\begin{align}
\label{estimate-Riemmanian-curv-tensor}
\nonumber
\| R_{M}( \bar{g}_l) - R_{M}(g_{0,l}) \|_{C^{\frac{\alpha}{2k}, \alpha}([0,T]\times [0,1])} \leq &  \|R_M\|_{C^{1}(M)}  \| \bar{g}_l-g_{0,l}\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\
\leq & C  \left(m, n, K_0, \|R_M\|_{C^{1}(M)} \right) \cdot T^{\beta}
\text{,}
\end{align}

\begin{align}\label{estimate-W_i}
\nonumber
&\|W_i(\partial_y^{i}\bar{g}_l, \ldots,\bar{g}_l)- W_i(\partial_y^{i}g_{0,l}, \ldots,g_{0,l})\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\ \nonumber
&\leq C(m, n, K_0) \cdot
\|W_i \|_{C^{1}(\mathbb{R}^n \times \cdots \times \mathbb{R}^n \times M) } \cdot \sum_{j \leq i} \| \partial_y^{j}\bar{g}_l-\partial_y^{j}g_{0,l}\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])}
\\
& \leq  C(m, n, K_0, \|\Pi_M \|_{C^{2k-1}(M)}) T^{\beta}
\end{align}

for some $\beta \in (0,1)$, $\forall\, i \leq 2k-1$.
From (<ref>), (<ref>),  (<ref>), and (<ref>),
one can apply the triangle inequalities in the
Hölder spaces, the identity

\begin{align}\label{identity-a^n-b^n}
a^p-b^p=(a-b)(a^{p-1}+a^{p-2} b+\cdots+a b^{p-2}+b^{p-1}),
~~~~~~ p \in \{2,3,\cdots\},
\end{align}

and the technical lemmas,
Lemmas <ref>$\sim$<ref>,
to obtain

\begin{align*}
& \sum\limits^{q}_{l=1}\Vert G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \bar{g}_{l})
-G_l(\partial_y^{2k-1} g_{0, l}, \ldots, g_{0, l})
\Vert_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\  \nonumber
& ~\leq C \left(\sum\limits^{q}_{l=1}\Vert g_{0,l} \Vert_{C^{2k,\alpha}([0,1])}, K_0, \|\Pi_M \|_{C^{2k-1}(M)}  \right) \sum\limits^{q}_{l=1} \| R_{M}( \bar{g}_l) - R_{M}(g_{0,l}) \|_{C^{\frac{\alpha}{2k}, \alpha}([0,T]\times [0,1])} \\
&~ +C \left(\sum\limits^{q}_{l=1}\Vert g_{0,l} \Vert_{C^{2k,\alpha}([0,1])}, K_0, \|R_M\|_{C^{1}(M)}  \right) \times \\
&~~~~~ \times \sum\limits_{j \leq 2k-1}\sum\limits^{q}_{l=1}  \|W_i(\partial_y^{i}\bar{g}_l, \ldots,\bar{g}_l)- W_i(\partial_y^{i}g_{0,l}, \ldots,g_{0,l})\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\
& ~~~~~ +C \left(\lambda, \sum\limits^{q}_{l=1}\Vert g_{0,l} \Vert_{C^{2k,\alpha}([0,1])}, K_0, \|\Pi_M \|_{C^{2k-1}(M)}, \|R_M\|_{C^{0}(M)} \right) \times \\
&~~~~~~~~~~~~~~~~~~~~~~ \times \sum\limits_{i \leq 2k-2} \sum\limits_{l=1}^{q} \| \partial_y^{j}\bar{g}_l - \partial_y^{j}g_{0, l}\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\
\leq & C\left(m, n, \lambda, \sum\limits^{q}_{l=1} \|g_{0, l}\|_{C^{2k,\alpha}([0,1])}, K_0,
\|R_M\|_{C^{1}(M)}, \|\Pi_M \|_{C^{2k-1}(M)} \right) T^{\beta}
\end{align*}

where $\beta \in (0,1)$.
Building upon the inequality
in Theorem <ref> and leveraging the triangle inequalities within the Hölder spaces,
we obtain

\begin{align*}
&\norm{g}_{X^{T}_{g_0} }
\leq  C_{0}\cdot \sum\limits^{q}_{l=1}
\Vert G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \bar{g}_{l})
-G_l(\partial_y^{2k-1} g_{0, l}, \ldots, g_{0, l})
\Vert_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\  \nonumber
&~~~~~~~~~~~~~ +C_{0}\cdot \sum\limits^{q}_{l=1}\Vert G_l(\partial_y^{2k-1} g_{0, l}, \ldots, g_{0, l})  \Vert_{C^{0,\alpha}([0,1])} \\ \nonumber
&~~~~~~~~~~~~~ +C_{0}\cdot \bigg( \sum\limits_{l=1}^{q} (| p_{l-1}|+|p_{l}|)
+\sum\limits_{\mu=1}^{k-1} ( |b^{\mu}_{1,0}|+ |b^{\mu}_{q, 2 \{ \frac{q}{2} \}}|  )
+\sum\limits^{q}_{l=1}\Vert g_{0, l} \Vert_{C^{2k,\alpha}([0,1])} \bigg)
\\
\leq & C_1 T^{\beta} +C_0 \cdot
\bigg( \sum^{q}_{l=1} \Vert G_l(\partial_y^{2k-1}g_{0, l}, \ldots, g_{0, l})  \Vert_{C^{0,\alpha}([0,1])}
+\sum\limits_{l=1}^{q} (| p_{l-1}|+|p_{l}|)
\\
+\sum\limits_{\mu=1}^{k-1} ( |b^{\mu}_{1,0}|+ |b^{\mu}_{q, 2 \{ \frac{q}{2} \}}| )
+\sum^{q}_{l=1}\Vert g_{0,l} \Vert_{C^{2k,\alpha}([0,1])} \bigg),
\end{align*}

where $T>0$ is sufficiently small,
$\beta=\min \big{ \{ }\frac{1-\alpha}{2k}, \frac{\alpha}{2k} \big{ \} }$ and

=C_1  \left(m, n, \lambda, \sum\limits^{q}_{l=1} \|g_{0, l}\|_{C^{2k,\alpha}([0,1])}, K_0,
\|R_M\|_{C^{1}(M)}, \|\Pi_M\|_{C^{2k-1}(M)} \right).

Now, we obtain $\norm{g}_{X^{t_{2}}_{g_0}}\leq K_0$ by choosing $K_0 \in (0, \infty)$ fulfilling

\begin{align*}
\frac{K_0}{2}=& C_0 \cdot \bigg{ ( } \sum^{q}_{l=1}
\Vert G_l(\partial_y^{2k-1}g_{0, l}, \ldots, g_{0, l})  \Vert_{C^{0,\alpha}([0,1])}
+\sum\limits_{l=1}^{q} (| p_{l-1}|+|p_{l}|)
\\
+\sum^{q}_{l=1}\Vert g_{0,l}\Vert_{C^{2k,\alpha}([0,1])} \bigg),
\end{align*}

and by choosing $t_2 \in (0, t_1)$ fulfilling

\begin{align*}
C_1 t^{\beta}_{2} \le \frac{K_0}{2}.
\end{align*}

Therefore, we conclude that

$$\mathcal{G}\left(X^{T}_{g_{0}} \cap B_{K_0}\right) \subset X^{T}_{g_{0}} \cap B_{K_0},
~~~~~~\forall\, T \in (0,t_2]. $$

$\bullet$ Contraction-maps.
We will show that
$\mathcal{G} : X^{T}_{g_{0}} \cap B_{K_0} \to X^{T}_{g_{0}} \cap B_{K_0}$
is a contraction-map for some $T>0$.
More precisely, we prove that there exists
$t_1 \in (0,t_2]$ such that

\begin{align}
\label{est:contraction-map}
\norm{g-h}_{X^{T}_{g_0} }
\leq C_{2} \, T^{\beta} \,  \norm{\bar{g}-\bar{h}}_{X^{T}_{g_0}},
\end{align}

$\forall\, T \in (0,t_1]$,
$\forall$ $\bar{g},\bar{h}\in X^{T}_{g_{0}}\cap B_{K_0}$, and $g=\mathcal{G}(\bar{g})$, $h=\mathcal{G}(\bar{h})$,
where $\beta \in (0,1)$

C_{2}=C_{2}\left(m, n, \lambda, \sum\limits^{q}_{l=1} \|g_{0, l}\|_{C^{2k,\alpha}([0,1])}, K_0, \|R_M\|_{C^{1}(M)}, \|\Pi_M\|_{C^{2k-1}(M)} \right)>0.

Observe that $g-h$ fulfills the linear system

\begin{equation}\label{eq:linear-g-h}
\begin{cases}
\partial_t ( g_{l}-h_{l})+(-1)^{k} \partial_y^{2k}( g_{l}-h_{l})=
G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \bar{g}_{l})-G_l(\partial_y^{2k-1} \bar{h}_{l}, \ldots, \bar{h}_{l}) \\
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \text{ in } (0,T) \times (0,1), ~ l \in\{1,\ldots,q\},
\\
(g_{l}-h_{l})(0,y)= 0,  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \forall\, y \in [0,1], l \in \{1,\ldots, q\},\\
(g_{l}-h_{l})(t, y^{\ast})=0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (t, y^{\ast}) \in [0,T]\times \{0,1\}, ~~ l \in \{1,\ldots, q\},
\\
\partial_{y}^{\mu}(g_{l}-h_{l})(t, 0)=0, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  t \in [0,T] , ~~ 1 \leq \mu \leq k-1, ~~ l \in \{1, q\},
\\
\partial_{y}^{\mu}\left(g_{2l-y^{\ast}}-h_{2l-y^{\ast}}\right) (t,y^{\ast})
+(-1)^{\mu-1}\partial_{y}^{\mu} \left(g_{2l-y^{\ast}+1}-h_{2l-y^{\ast}+1}\right)(t,y^{\ast})=0,\\
~~~~~~~~~~~~ (t,y^{\ast}) \in [0,T] \times \{0,1\},~~1 \leq \mu \leq 2k-2, ~~ l \in  \{1, \ldots, \left[ \frac{q-1}{2}\right] +\left(1-2 \big{ \{ } \frac{q}{2} \big{ \}}\right) y^{\ast}  \}.
\end{cases}
\end{equation}

Given that $g = \mathcal{G}(\bar{g})$, $h = \mathcal{G}(\bar{h})$, and their initial data $\bar{g}_{l}(0,\cdot)=g_{0,l}(\cdot)$, $\bar{h}_{l}(0,\cdot)=g_{0,l}(\cdot)$ satisfy the compatibility conditions of order $0$ to (<ref>) with the initial-boundary conditions (<ref>) or (<ref>) (for $q$ odd), the initial datum of $g-h$ is null and adheres to the compatibility conditions of order $0$ for the linear parabolic system (<ref>). Consequently, we can apply the theory of linear parabolic systems from <cit.> to (<ref>), leading to the conclusion that $g-h$ is the unique solution to the linear system with the provided estimate.

\begin{align}
&\norm{g-h}_{X^{T}_{g_0}} \leq
C_0 \cdot \left(\sum\limits^{q}_{l=1}
\norm{ G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \bar{g}_{l})-G_l(\partial_y^{2k-1} \bar{h}_{l}, \ldots, \bar{h}_{l})}_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])}
\right) \text{.}
\label{est:contraction-map1}
\end{align}

By using (<ref>), the triangle inequality in the
Hölder spaces, (<ref>), (<ref>), and leveraging the technical lemmas,
Lemmas <ref>$\sim$<ref>,
with $v=\bar{g}_l-\bar{h}_l$, we can establish the following estimates:

\begin{align}\label{est:Riemmanian-curv-tensor-g}
\nonumber
&\| R_{M}( \bar{g}_l) - R_{M} (\bar{h}_l) \|_{C^{\frac{\alpha}{2k}, \alpha}([0,T]\times [0,1])} \\
& ~~~~~~~~~~~~~~~~~~~\leq
C \left(m, n, K_0, \|R_M\|_{C^{1}(M)} \right)
\cdot T^{\beta} \cdot \| \bar{g}_l-\bar{h}_l \|_{C^{\frac{2k+\alpha}{2k}, 2k+\alpha}([0,T]\times [0,1])},
\end{align}

\begin{align}\label{est:W_i-g-g}
\nonumber
&\|W_i(\partial_y^{i}\bar{g}_l, \ldots,\bar{g}_l)- W_i(\partial_y^{i}\bar{h}_l, \ldots,\bar{h}_l)\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\ \nonumber
C(m, n, K_0)
\cdot \|W_i\|_{C^1(\mathbb{R}^n \times \cdots \times \mathbb{R}^n \times M)} \cdot \sum\limits_{j \leq i} \| \partial_y^{i}\bar{g}_l-\partial_y^{i} \bar{h}_l\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])}
\\
& ~~~~~~~~ \leq
C\left((m, n, K_0, \|\Pi_M \|_{C^{2k-1}(M)} \right) \cdot T^{\beta} ,
\end{align}

where $\beta \in (0,1)$, $i \leq 2k-1$.
By applying (<ref>), (<ref>),  (<ref>), (<ref>), (<ref>), and
the technical lemmas,
Lemmas <ref>$\sim$<ref>, we obtain

\begin{align*}
\nonumber
&\norm{ G_l(\partial_y^{2k-1} \bar{g}_{l}, \ldots, \bar{g}_{l})-G_l(\partial_y^{2k-1} \bar{h}_{l}, \ldots, \bar{h}_{l})}_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\ \nonumber
&~~ \leq C \left( K_0, \|\Pi_M \|_{C^{2k-1}(M)}  \right) \sum\limits^{q}_{l=1} \| R_{M}( \bar{g}_l) - R_{M}(\bar{h}_{l}) \|_{C^{\frac{\alpha}{2k}, \alpha}([0,T]\times [0,1])} \\
&~~~~~~~~  +C \left( K_0, \|R_M\|_{C^{1}(M)}  \right) \sum\limits_{j \leq 2k-1}\sum\limits^{q}_{l=1}  \|W_i(\partial_y^{i}\bar{g}_l, \ldots,\bar{g}_l)- W_i(\partial_y^{i}\bar{h}_{l}, \ldots, \bar{h}_{l})\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\
& ~~~~~~~~ +C \left(\lambda, K_0, \|\Pi_M \|_{C^{2k-1}(M)}, \|R_M\|_{C^{0}(M)} \right) \sum\limits_{i \leq 2k-2} \sum\limits_{l=1}^{q} \| \partial_y^{j}\bar{g}_l - \partial_y^{j}\bar{h}_{l}\|_{C^{\frac{\alpha}{2k}, \alpha}([0,T] \times [0,1])} \\
& ~~ \leq  C\left(m, n, \lambda, K_0,
\|R_M\|_{C^{1}(M)}, \|\Pi_M \|_{C^{2k-1}(M)} \right)
\cdot T^{\beta} ,
\end{align*}

where $\beta \in (0,1)$.
Now, together with (<ref>),
we obtain  (<ref>).
Hence, the operator

$$\mathcal{G} :X^{T}_{g_{0}} \cap B_{K_0} \to X^{T}_{g_{0}} \cap B_{K_0}$$

not only represents a self-map but also serves as a strict contraction map for all $T\in(0,t_1]$.
By choosing any $t_0 \in (0,t_1)$, where

\begin{align*}
t_0=t_0\left(m, n, \alpha, \lambda, \sum\limits^{q}_{l=1} \|g_{0, l}\|_{C^{2k,\alpha}([0,1])}, K_0, \|R_M\|_{C^{1}(M)}, \|\Pi_M\|_{C^{2k-1}(M)}  \right),
\end{align*}

we then apply the Banach fixed-point theorem to the set $X^{t_{0}}_{g_{0}} \cap B_{K_0}$ and conclude that there exists a unique fixed point
$g \in X^{t_{0}}_{g_{0}} \cap B_{K_0}$.
Note that the fixed point $g$ is also a solution to
(<ref>) with the initial-boundary conditions
(<ref>) or (<ref>).

$\bullet$ Smoothing,
i.e., $g_{l} \in C^{\infty}\left( (0,t_0] \times [0,1] \right)$, $\forall\, l$.
The argument follows the similar approach in the proof of <cit.>, which is based on constructing smooth cut-off functions with null initial data, verifying the compatibility conditions of higher-order linear parabolic PDEs <cit.>, and applying bootstrapping on the regularity.
In essence, the regularity of the map $g$ can be improved by demonstrating that

\begin{align*}
g_{l} \in \bigcap_{\ell=1}^{\infty}C^{\frac{2k+\ell+\alpha}{2k},2k+\ell+\alpha}\left(\left[ \frac{2\ell-1}{2\ell}\varepsilon, t_0\right] \times [0,1]\right)
\end{align*}

for any $\varepsilon \in (0,t_0)$.
For further details, the reader is referred to [15, 22].

From (<ref>), we have

\begin{align*}
\gamma_{l}(t,x)=g_{l}\left(t, (-1)^{l+1}
\left(x-x_{2\left[\frac{l}{2}\right]} \right)\right) ,
~~~~~~~~~~x_{l-1} \leq x \leq x_{l}, ~~ l \in \{1,\ldots, q\}.
\end{align*}

The proof simply follows from applying Theorem <ref>.

§.§.§ The global solutions and the asymptotics

For tangent vector fields $v, v_{1},\ldots, v_{k}$ on $M$,
we denote by

\begin{equation*}
v_{1}\ast \cdot\cdot\cdot \ast v_{k}=\left\{
\begin{array}{l}
\langle v_{i_{1}}, v_{i_{2}}\rangle \cdot \cdot \cdot \langle v_{i_{k-1}}, v_{i_{k}}\rangle \text{ \ \ \ \ \ \ \ \ ,\ for }k\text{\ even,} \\
\langle v_{i_{1}}, v_{i_{2}}\rangle \cdot \cdot \cdot \langle v_{i_{k-2}}, v_{i_{k-1}}\rangle \cdot v_{i_{k}}
\text{ , for }k\text{\, odd,}
\end{array}
\right.
\end{equation*}

where $i_{1},\ldots,i_{k}$ is any permutation of $1,\ldots,k$, $\langle\cdot,\cdot\rangle=\langle\cdot,\cdot\rangle_g$, and $g$ is the Riemannian metric of $M$.
Slightly more generally, we allow some of the tangent vector fields $v_{i}$ to be scalar functions, in which the $\ast $-product reduces to multiplication.
The notation $P^{a,c}_b(v)$ represents the summation of finitely many terms
$D_{x}^{\mu_1} v \ast \cdots \ast D_{x}^{\mu_b} v $
with $\mu_1+\cdots+\mu_b=a$, $\max\{\mu_1, \ldots, \mu_b\}\le c$, and with universal coefficients,
depending on the Riemannian curvature tensor $R_M$ and Riemannian metric $g_M$ and their derivatives.

The universal coefficients and the numbers of terms in $P_b^{a,c}(v)$
are assumed to be bounded from above and below by constants depending only on
$a$, $b$, $K_0$, where $K_0$ is determined by the $C^{\infty}$-norm of the Riemannian curvature tensor $R_M$.
Note that it is easy to verify the following:

\begin{equation*}
\left\{
\begin{array}{l}
D_{x}\left( P_{b_{1}}^{a_{1},c_{1}}\left( v \right) \ast P_{b_{2}}^{a_{2},c_{2}}\left( v \right) \right)
=D_{x} P_{b_{1}}^{a_{1},c_{1}}\left( v \right) \ast P_{b_{2}}^{a_{2},c_{2}}\left( v \right)
+P_{b_{1}}^{a_{1},c_{1}}\left( v \right) \ast D_{x} P_{b_{2}}^{a_{2},c_{2}}\left( v \right) \text{,}
\\ \\
P_{b_{1}}^{a_{1},c_{1}}\left( v \right) \ast P_{b_{2}}^{a_{2},c_{2}}\left( v \right)
=P_{b_{1}+b_{2}}^{a_{1}+a_{2},max\{c_1,c_2\}}\left( v \right) \text{, }D_{x}P_{b_{2}}^{a_{2},c_{2}}\left( v \right)
=P_{b_{2}}^{a_{2}+1,c_{2}+1}\left( v \right)
\text{.}
\end{array}
\right.
\end{equation*}

The notation $Q_{b}^{a, c}(\gamma_{x})$ is introduced to simplify the summation of polynomials. Namely,

\begin{equation*}
Q_{b}^{a, c}(\gamma_{x})
:= \sum_{ \substack{  [\![\mu, \nu]\!]\leq [\![a,b]\!] \\\xi \leq c} }   P_{\nu}^{\mu, \xi}(\gamma_{x})
=\sum_{\mu=0}^{a}\sum_{\nu=1}^{2a+b-2\mu}\sum_{\xi=0}^{c}\text{ } P_{\nu}^{\mu, \xi}(\gamma_{x})
\text{,}
\end{equation*}

$[\![\mu,\nu]\!]:=2\mu+\nu$, and $a, c \in \mathbb{N}_{0}$, $b \in \mathbb{N}$.
For simplicity, we also use the notation $P^{a}_{b}\left(\gamma_{x} \right)$ and
to represent $P^{a,a}_{b}\left(\gamma_{x} \right)$ and
$Q_{b}^{a, a}(\gamma_{x})$

The notation $D^{0}_{x}\gamma=\gamma$ and $D^{0}_{x}\gamma_{l}=\gamma_{l}$ represent an extension for simplicity, allowing the use of the covariant differentiation notation $D_x$ in a way that $D_{x}\gamma=\gamma_{x}$ and $D_{x}\gamma_{l}=\gamma_{l,x}$.
A similar extension applies to covariant differentiation with respect to $t$, denoted as $D_t$.

Let $\gamma : (t_0,t_1)\times (0,1)\rightarrow M$ be a smooth map,
and $V: (t_0,t_1)\times (0,1)\rightarrow T_{\gamma}M$ be a smooth vector field associated with $\gamma$.

\begin{align}
D_{t} D_{x} \gamma=& D_{x} D_{t} \gamma
\text{,}
\label{eq:D_t-D_x^j=D_x^j-D_t}
\\
D_{t}D_{x}^{\mu}V=& D_{x}^{\mu}D_{t}V
+\sum_{\nu=1}^{\mu} D_{x}^{\mu-\nu}\left[R(\gamma_{t}, \gamma_{x}) D_{x}^{\nu-1} V \right],
~~\forall\, \mu \geq 1
\text{,}
\label{eq:D_t-D_x^j->D_x^j-D_t}
\\
D_{t}^{\mu} D_{x}V=& D_{x}D_{t}^{\mu}V +\sum_{\nu=1}^{\mu} D_{t}^{\mu-\nu}\left[R(\gamma_{t}, \gamma_{x}) D_{t}^{\nu-1} V \right],
~~\forall\, \mu \geq 1
\text{.}
\label{eq:D_t^j-D_x->D_x-D_t^j}
\end{align}

By applying Nash's isometric embedding theorem, we have

\begin{align*}
=\partial_{t}\partial_{x}\gamma-\sum_{\alpha=m+1}^{n}\langle\partial_{t}\partial_{x}\gamma, e_\alpha\rangle e_\alpha,
\\
=\partial_{x}\partial_{t}\gamma-\sum_{\alpha=m+1}^{n}\langle\partial_{x}\partial_{t}\gamma, e_\alpha\rangle e_\alpha.
\end{align*}

Since $\partial_{t}\partial_{x}\gamma=\partial_{x}\partial_{t}\gamma$, we conclude (<ref>).

We prove (<ref>) inductively.
As $\mu=1$, the equality

\begin{align}
\label{eq:D_t-D_x->D_x-D_t}
D_{t}D_{x}V=&D_{x}D_{t}V+R( \gamma_{t}, \gamma_{x})V
\end{align}

comes from using the definition of Riemannian curvature tensors $R(\boldsymbol{\cdot},\boldsymbol{\cdot})\boldsymbol{\cdot}$.
By an induction argument and applying
(<ref>), we obtain

\begin{align*}
D_{t}D_{x}^{\mu+1}V=&D_{t}D_{x} \left(D_{x}^{\mu}V\right)
=D_{x}D_t \left(D_{x}^{\mu}V\right)+R(\gamma_{t}, \gamma_{x})D_{x}^{\mu}V \\
=&D_{x} \left(D_{x}^{\mu}D_{t}V
+\sum_{\nu=1}^{\mu} D_{x}^{\mu-\nu}\left[R( \gamma_{t}, \gamma_{x}) D_{x}^{\nu-1} V \right] \right)
+R(\gamma_{t}, \gamma_{x})D_{x}^{\mu}V \\
+\sum_{\nu=1}^{\mu+1} D_{x}^{\mu+1-\nu}\left[R(\gamma_{t}, \gamma_{x}) D_{x}^{\nu-1} V \right],
\end{align*}

which finishes the proof of

Below, we prove (<ref>).
The case $\mu=1$ is proven in (<ref>).
By an induction argument and (<ref>),
we obtain

\begin{align*}
D_{t}^{\mu+1}D_{x}V=&D_{t} \left(D_{x}D_{t}^{\mu}V
+\sum_{\nu=1}^{\mu} D_{t}^{\mu-\nu}\left[R(\gamma_{t}, \gamma_{x}) D_{t}^{\nu-1} V \right]\right)\\
=&D_{x}D_t \left(D_{t}^{\mu}V\right)+R(\gamma_{t}, \gamma_{x})D_{t}^{\mu}V
+\sum_{\nu=1}^{\mu} D_{t}^{\mu-\nu+1}\left[R(\gamma_{t}, \gamma_{x}) D_{t}^{\nu-1} V \right] \\
=&D_{x} D_{t}^{\mu+1}V
+\sum_{\nu=1}^{\mu+1} D_{t}^{\mu+1-\nu}\left[R(\gamma_{t}, \gamma_{x}) D_{t}^{\nu-1} V \right],
\end{align*}

which finishes the proof of (<ref>).

$\gamma_{l}:(t_0,t_1)\times I_{l}\rightarrow M$ is a smooth solution of
Then, for any $\mu \in\mathbb{N}$, $\ell \in\mathbb{N}_{0}$,
we obtain the following formula,

\begin{align*}
D_{t}^{\mu} D^{\ell}_{x} \gamma_{l}=&(-1)^{(k+1) \mu} D_{x}^{2k\mu+\ell-1} \gamma_{l,x}+Q^{2k\mu+\ell-3}_{3}(\gamma_{l,x}).
\end{align*}

Below, we prove the lemma inductively on $\mu$.
For $\mu=1$ and $\ell=0$, it follows directly from applying (<ref>).
Now, consider the case when $\mu=1$ and $\ell \geq 2$. We apply (<ref>) and (<ref>), resulting in

\begin{align}\label{D_t-D^s_x}
\nonumber
D_{t} D^{\ell}_{x} \gamma_{l}=& D_{t} D^{\ell-1}_{x} \gamma_{l, x}=
D^{\ell-1}_{x}D_{t} \gamma_{l,x}+\sum^{\ell-1}_{r=1}D_{x}^{\ell-1-r}\left[R(\gamma_{l, t}, \gamma_{l, x}) D^{r-1}_{x}\gamma_{l, x}\right] \\  \nonumber
=&D^{\ell}_{x} \gamma_{l, t}+\sum^{\ell-1}_{r=1}D_{x}^{\ell-1-r}\left[R(\gamma_{l, t}, \gamma_{l, x}) D^{r-1}_{x}\gamma_{l, x}\right] \\
=& (-1)^{k+1}D_{x}^{2k+\ell-1}\gamma_{l, x}+Q^{2k+\ell-3}_{3}(\gamma_{l,x}).
\end{align}

Hence, the lemma holds for $\mu=1$.

By applying (<ref>) and the definition of $P^{a, c}_{b}(\gamma_{l,x})$ along with an induction argument, we obtain

\begin{align}
\label{general-formula-Q^a-c_b}
\end{align}

Using (<ref>) and (<ref>), we obtain

\begin{align*}
D_{t}^{\mu} D^{\ell}_{x} \gamma_{l} =& D_{t} (D_{t}^{\mu-1} D^{\ell}_{x} \gamma_{l} )
=D_{t}\left((-1)^{(k+1) (\mu-1)} D_{x}^{2k(\mu-1)+\ell-1} \gamma_{l,x}+Q^{2k(\mu-1)+\ell-3}_{3}(\gamma_{l,x})\right)\\
=&(-1)^{(k+1)\mu} D_{x}^{2k\mu+\ell-1} \gamma_{l,x}+Q^{2k\mu+\ell-3}_{3}(\gamma_{l,x}).
\end{align*}

This completes the proof.

Let $\gamma_{l}: (t_0,t_1)\times I_{l}\rightarrow M$ be a smooth solution to
(<ref>) for any $l\in\{1,\cdots,q\}$.
Then, we have the followings.

\begin{align}
\label{BC-1}
D^{\mu}_{x}D_{t}^{\ell} \gamma_{1}(t,x_0)=0,
D^{\mu}_{x}D_{t}^{\ell} \gamma_{q}(t,x_q)=0,
~ \forall \, \mu \in \{0, \ldots, k-1\}, \, \ell \in \mathbb{N}, \,  t \in (t_0,t_1)
\text{.}
\end{align}

$l \in \{1, \ldots, q-1\}$,

\begin{align}
\label{BC-2}
~ \forall\, \mu \in \{0,\ldots, 2k-2\}, \, \ell \in \mathbb{N}_{0}, \, t \in (t_0,t_1)
\text{.}
\end{align}

We prove (<ref>) and (<ref>) inductively on $\ell$.

As $\ell=1$, we apply (<ref>) and (<ref>) to obtain

\begin{align} \label{eq:identity-1}
D_{x}^{\mu} D_{t} \gamma_{l}=D^{\mu-1}_{x} D_{t} D_x \gamma_{l}
= & D_t D^{\mu-1}_{x} \gamma_{l,x}
- \sum_{\nu=1}^{\mu-1} D_{x}^{\mu-1-\nu}\left[R(\gamma_{l,t}, \gamma_{l,x}) D_{x}^{\nu-1} \gamma_{l,x}  \right]
\text{.}
\end{align}

From applying the boundary conditions
(<ref>), (<ref>) at the two end-points,
the right hand side of (<ref>) equals zero.
Thus, (<ref>) holds as $\ell=1$.

As $\ell\ge 2$, we apply (<ref>) to obtain

\begin{align}
\label{eq:identity-2}
D^{\mu}_{x}D_{t}^{\ell} \gamma_{l}
=&D^{\mu}_{x}D_t (D_{t}^{\ell-1} \gamma_{l} )
=D_t D^{\mu}_{x} (D_{t}^{\ell-1} \gamma_{l} )
- \sum_{\nu=1}^{\mu} D_{x}^{\mu-\nu}\left[R(\gamma_{l,t}, \gamma_{l,x}) D_{x}^{\nu-1} D_{t}^{\ell-1} \gamma_{l} \right]
\text{.}
\end{align}

By applying the boundary conditions (<ref>), (<ref>) at the two end-points, and by an induction argument on $\ell$,
the right hand side of (<ref>) vanishes.
Therefore, (<ref>) holds.

The proof of (<ref>) is also an induction argument. We proceed as above by applying the boundary conditions (<ref>) and (<ref>) to (<ref>)
to conclude that the right hand side of (<ref>) is continuous at the boundary points $\{x_1,\cdots,x_{q-1}\}$, which offers the required continuity. Thus, (<ref>) is proven.

Let $\gamma: [t_0,t_1) \times I \rightarrow M$ be a smooth map and $\ell \in \mathbb{N}_{0}$.
Assume that the Riemannian curvature tensor $R_M$ is smooth and uniformly bounded in $C^{\infty}$-norm.
If $A, B\in \mathbb{N}$ with $B \geq 2$ and $A+\frac12 B<2\ell+5$, then we have

\begin{align*}
\sum_{\substack{[\![a,b]\!] \leq [\![A,B]\!]\\c \leq \ell+2, ~ b\geq 2}} \int_{I}\big\vert P^{a,c}_{b}( \gamma_{x}  ) \big\vert\text{ }dx \leq
C \cdot \| \gamma_{x}  \|_{\ell+2, 2}^{\alpha} \, \| \gamma_{x}  \|_{L^2(I)}^{B-\frac{A}{\ell+2}}
\end{align*}

where  $\alpha= (A + \frac12 B-1)/(\ell+2)$ and $C=C(n, \ell, A, B, N_{L})$
for some constant $L\in\mathbb{N}_0$ such that $\|R_M\|_{C^{L}} \leq N_L$.

Assume that $i_{1}, \ldots, i_{b} \in \{0, 1, \ldots, c\}$, $i_{1}+\cdots+i_{b}=a$.
From the definition of $P^{a,c}_{b}( \gamma_{x} )$
and by applying the Hölder's inequality, and Lemma <ref> (vi.),
we obtain

\begin{align*}
\nonumber
&\int_{I}| P^{a,c}_{b}( \gamma_{x}  )| \text{ }dx
\leq C \cdot \int_{I} \prod^{b}_{\mu=1}|D_{x}^{i_{\mu}} \gamma_{x} |
\text{ }dx
\leq C \cdot \prod^{b}_{\mu=1} \| D_{x}^{i_{\mu}} \gamma_{x}  \|_{L^{b}(I)}
\\
& ~~~~~~~~~
\leq C \cdot \prod^{b}_{\mu=1} \|  \gamma_{x}  \|_{\ell+2, 2}^{\frac{1}{\ell+2}\left[i_{\mu}+\frac{1}{2}-\frac{1}{b}\right]}
\cdot \| \gamma_{x}  \|_{L^2(I)}^{\frac{\ell+2-i_{\mu}}{\ell+2}}
\leq C \cdot \|   \gamma_{x}  \|_{\ell+2, 2}^{\frac{1}{\ell+2}\left(a+\frac{b}{2}-1\right)}
\cdot \|  \gamma_{x}  \|_{L^2(I)}^{b-\frac{a}{\ell+2}}
\end{align*}

where $C=C\left(m,n, \ell, a, b, N_L\right)$.
The proof of the lemma follows is obtained from applying the inequalities.

From Theorem <ref>,
we only need to prove the long-time existence below.
The proof is an argument by contradiction.
From the short-time existence, we assume on the contrary that any smooth solution to
(<ref>)$\sim$(<ref>) only exists up to $t_1\in [t_0, \infty)$, i.e.,
the smoothness of $\gamma_{l}(t_1,\cdot)$ fails for some $t_1$ and $l$.
Then, we show that
$\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}(t,\cdot)\|^2_{L^2(I_{l})}$
remains uniformly bounded up to $t_1$ for any $\mu \in \mathbb{N}$.
This will give the contradiction.

By applying Lemmas <ref>, <ref>, and (<ref>), we obtain

\begin{align}
\nonumber
\frac{d}{dt}& \frac{1}{2}\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}\|^2_{L^2(I_{i})}
=\sum_{l=1}^{q}\int_{I_{l}}\langle D_t^{\mu} \partial_{t} \gamma_{l}, D_t^{\mu-1} \partial_{t} \gamma_{l} \rangle \text{ }dx \\ \nonumber
=& (-1)^{k+1} \sum_{l=1}^{q}\int_{I_{l}}\big{\langle} D_t^{\mu} D_{x}^{2k-1}\gamma_{l, x}, D_t^{\mu-1} \partial_{t} \gamma_{l} \big{\rangle}\text{ }dx
\\ \nonumber
&+  \sum\limits_{l=1}^{q} \sum\limits^{k-1}_{r=2}  (-1)^{r+k+1} \int_{I_{l}}\bigg{\langle}
D_t^{\mu} \left[ R\left(D^{2k-r-1}_{x}\gamma_{l, x}, D^{r-2}_{x} \gamma_{l, x} \right) \gamma_{l, x}
+\lambda D_x\gamma_{l,x} \right], D_t^{\mu-1} \partial_{t} \gamma_{l} \bigg{\rangle}\text{ }dx
\\
=& I_1+I_2+\sum_{l=1}^{q}\int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu+k-1}_{4}(\gamma_{l,x}) \text{ }dx,
\label{identity-D-t-j-1}
\end{align}

\begin{equation*}
\begin{cases}
I_1=(-1)^{k+1} \sum\limits_{l=1}^{q}\int\limits_{I_{l}}\big{\langle} D_t^{\mu}  D_{x}^{2k-1}\gamma_{l, x}, D_t^{\mu-1} \partial_{t} \gamma_{l}\big{\rangle}\text{ }dx,\\
I_2= \sum\limits_{l=1}^{q} \sum\limits^{k-1}_{r=2}  (-1)^{r+k+1} \int\limits_{I_{l}}\big{\langle}
R\left(D_t^{\mu}D^{2k-r-1}_{x}\gamma_{l, x} , D^{r-2}_{x}\gamma_{l, x}\right) \gamma_{l, x}, D_t^{\mu-1}\partial_{t} \gamma_{l} \big{\rangle}\text{ }dx,
\end{cases}
\end{equation*}

are terms involving high-order derivatives.
We need to apply integration by parts to these terms so that the Gronwall  inequality can be obtained from using interpolation inequalities.

By letting $V=D_{x}^{2k-2}\gamma_{l, x}$ in  (<ref>) of Lemma <ref>, by applying Lemma <ref>, and using integration by parts,
we obtain

\begin{align*}
=& (-1)^{k+1}\left(\sum_{l=1}^{q}\int_{I_{l}} \big{\langle} D_{x}D_t^{\mu}  D_{x}^{2k-2}\gamma_{l, x} , D_t^{\mu-1} \partial_{t} \gamma_{l} \big{\rangle}\text{ }dx
+\sum_{l=1}^{q} \int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu-1}_{4}(\gamma_{l,x}) \text{ }dx \right)
\\
=& (-1)^{k+1}\sum_{l=1}^{q}\big{\langle} D_t^{\mu}  D_{x}^{2k-2}\gamma_{l, x}, D_t^{\mu-1} \partial_{t} \gamma_{l} \big{\rangle} \bigg{|}^{x_{l}}_{x_{l-1}} \\
-& (-1)^{k+1}\left(\sum_{l=1}^{q}\int_{I_{l}} \big{\langle} D_t^{\mu}  D_{x}^{2k-2}\gamma_{l, x} , D_{x} D_t^{\mu-1} \partial_{t} \gamma_{l} \big{\rangle}\text{ }dx
-\sum_{l=1}^{q} \int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu-1 }_{4}(\gamma_{l,x}) \text{ }dx \right)
\\
=&(-1)^{k} \left( \sum_{l=1}^{q}\int_{I_{l}} \big{\langle} D_t^{\mu}  D_{x}^{2k-2}\gamma_{l, x}, D_{x}D_t^{\mu-1} \partial_{t} \gamma_{l} \big{\rangle}\text{ }dx
-\sum_{l=1}^{q} \int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu-1 }_{4}(\gamma_{l,x}) \text{ }dx \right)
\text{,}
\end{align*}

where the boundary terms vanish as applying
Lemma <ref>, i.e.,

\begin{align*}
&\big{\langle} D_t^{\mu}  D_{x}^{2k-2}\gamma_{q, x}(t,x_{q}), D_t^{\mu-1} \partial_{t}\gamma_{q}(t,x_{q}) \big{\rangle} =0,
\\
&\big{\langle} D_t^{\mu}  D_{x}^{2k-2}\gamma_{1, x}(t,x_{0}), D_t^{\mu-1}\partial_{t}\gamma_{1}(t,x_{0}) \big{\rangle} =0,
\\
&\big{\langle} D_t^{\mu}( \Delta_{l} D_{x}^{2k-2} \partial_{x} \gamma )(t), D_t^{\mu-1}\partial_{t} \gamma_{l}(t,x_{l}^{-}) \big{\rangle} =0,
~~~~~~ \forall\, l=1,\cdots, q-1.
\end{align*}

Note that the terms written as
$\sum\, \int \, Q^{a,c}_{b}(\gamma_{l,x})$
represent lower-order terms, which will be taken care of by applying the interpolation inequalities in Lemma <ref> later.
In other words, we obtain

I_1=(-1)^{k} \sum_{l=1}^{q}\int_{I_{l}} \big{\langle} D_t^{\mu}  D_{x}^{2k-2}\gamma_{l, x}, D_{x}D_t^{\mu-1} \partial_{t} \gamma_{l} \big{\rangle}\text{ }dx
+\text{ lower-order terms}.

By the same procedure as above,
i.e., applying Lemmas <ref>, <ref>, <ref>,
and integration by parts to the main terms in $I_1$
inductively, we finally derive

\begin{align*}
I_1=& -\sum_{l=1}^{q}\int_{I_{l}} \big{\langle} D_t^{\mu}  D_{x}^{k-1}\gamma_{l, x}, D_{x}^{k}D_t^{\mu-1} \partial_{t} \gamma_{l} \big{\rangle}\text{ }dx
+\sum_{l=1}^{q} \int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu+k-1}_{4}(\gamma_{l,x}) \text{ }dx \\
=&-\sum_{l=1}^{q}\int_{I_{l}} | D_{x}^{2k\mu+k-1}\gamma_{l, x} |^2 \text{ }dx
+\sum_{l=1}^{q} \int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu+k-1}_{4}(\gamma_{l,x}) \text{ }dx.
\end{align*}

Utilizing the symmetry property inherent in Riemannian curvature tensors, the terms within the Riemannian curvature tensors of $I_2$ can be interchanged up to sign changes.
This enables the application of integration by parts and Lemma <ref>, leading to the following result,

\begin{align*}
=\sum_{l=1}^{q} \int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu+k-1}_{4}(\gamma_{l,x}) \text{ }dx,
\end{align*}

which are lower-order terms.

Thus, by applying
Lemmas <ref>, <ref>, <ref> to (<ref>) and the uniform-bound estimate of the Dirichlet integral of $\gamma$ in
we obtain

\begin{align*}
&\frac{d}{dt}\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}\|^2_{L^2(I_{l})}
+\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}\|^2_{L^2(I_{l})}
+2 \sum_{l=1}^{q}\int_{I_{l}} | D_{x}^{2k\mu+k-1}\gamma_{l, x} |^2
\text{ }dx
\\
=&\sum_{l=1}^{q} \int_{I_{l}} Q^{4k\mu+2k-4, 2k\mu+k-1}_{4}(\gamma_{l,x}) \text{ }dx
+\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}\|^2_{L^2(I_{l})}
\\
\leq & \varepsilon \sum_{l=1}^{q}\int_{I_{l}} | D_{x}^{2k\mu+k-1}\gamma_{l, x} |^2 \text{ }dx
+C\left(m, n, \varepsilon, \lambda, k, \mu, N_{2k\mu+k-1},
\mathcal{E}_{k}^\lambda[\gamma(t_0,\cdot)] \right),
\end{align*}

where $\varepsilon>0$ is chosen arbitrarily small and
$\{N_0, N_1,\cdots,N_L,\cdots\}$ is a sequence of positive numbers such that
$\|R_M\|_{C^L}\le N_L$, $\forall\, L\in\mathbb{N}_0$.
By letting $\varepsilon\in(0,1]$, we obtain

\begin{align}\label{inequality-D^{j-1}_t-partial_{t}gamma_i}
\frac{d}{dt}\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}\|^2_{L^2(I_{l})}
&+\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}\|^2_{L^2(I_{l})}
\leq
C(m, n, \lambda, k, \mu, N_{2k\mu+k-1}, \mathcal{E}_{k}^\lambda[\gamma(t_0,\cdot)]).
\end{align}

By applying the Gronwall inequality, we obtain

\begin{align}\label{ineq:D-4k-higher-order}
\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}\|^2_{L^2(I_{l})}
\leq
C\left(m, n, k, \lambda, N_{2k\mu+k-1}, \mathcal{E}_{k}^{\lambda}[\gamma(t_0,\cdot)],
\sum\limits_{l=1}^{q}\| D_t^{\mu-1}\partial_{t} \gamma_{l}(t_0,\cdot)\|^2_{L^2(I_{l})} \right).
\end{align}

By applying Lemmas <ref>, <ref>, and <ref>,
the estimate (<ref>) is equivalent to

\begin{align}\label{unif-bounded:D-2kj+k-3-derivative}
\sum\limits_{l=1}^{q} \| D_{x}^{2k\mu-1} \gamma_{l,x} \|^2_{L^2(I_{l})}
\leq
C\left(m, n, k, \lambda, N_{2k\mu+k-1}, \mathcal{E}_{k}^{\lambda}[\gamma(t_0,\cdot)],
\sum\limits_{l=1}^{q}\| D_t^{\mu-1} \partial_{t} \gamma_{l}(t_0,\cdot)\|^2_{L^2(I_{l})}
\right).
\end{align}

The asymptotic behavior of the global solutions is established using a similar argument to that found in [15, 14]. For the sake of completeness, we restate the proof below, albeit concisely.
From the uniform bounds in (<ref>),
we may choose a subsequence of
$\gamma \left(t_{j},\cdot \right)$,
which converges smoothly to $\gamma_{\infty }(\cdot)$ as $t_{j}\rightarrow \infty$.

\begin{align*}
u\left( t\right)
:=\overset{q}{\underset{l=1}\sum}{\underset{I_l}{\int}}\,|\partial _{t}\gamma_l|^{2}\, dx
\text{.}
\end{align*}

From (<ref>), we obtain

\begin{equation*}
u\left( t\right) \leq
C\left( m,n,k, \lambda, N_{2k\mu+k-1}, \mathcal{E}_{k}^{\lambda}[\gamma(t_0,\cdot)] \right)
\text{, }
\end{equation*}

which implies $u \in L^{1}\left([t_0,\infty )\right)$.
Together with (<ref>),
we obtain

\begin{align*}
\left| u^{\prime }\left( t\right) \right|
\leq
C(m,n,k, \lambda, N_{2k\mu+k-1}, \mathcal{E}_{k}^\lambda[\gamma(t_0,\cdot)]).
\end{align*}

Hence $u\left( t\right) \rightarrow 0$ as $t\rightarrow \infty$.
Again, together with the smoothness of $\gamma_{\infty,l}$, $\forall\, l$,
we conclude that
$\gamma_{\infty }(\cdot) \in
\Phi_{\mathcal{P}}$
and fulfills
$\mathcal{L}^{2k}_{x}(\gamma_{\infty,l})=0$, $\forall\, l$.

The proof for showing $\gamma_{\infty} \in C^{2k-2,1}([x_0,x_q])$ simply follows from
the smoothness of $\gamma_{\infty,l} \in C^{\infty}([x_{l-1}, x_{l})]$,
and the uniformly bounds of $\|\gamma_{\infty,l}\|_{W^{\mu,2}([x_{l-1}, x_{l}])}$, $\forall\, \mu\in\mathbb{N}$.

§ APPENDIX

Here we introduce the (parabolic) Hölder spaces, by
following <cit.>. Let $g: [0,T] \times [0,1] \to \mathbb{R}$, and $\rho\in (0,1)$. We introduce the notation

\begin{align*}
[ g ]_{\rho,y}=&\underset{(t,y),(t,y^{\prime}) \in  [0,T] \times [0,1]}{\sup} \frac{| g(t,y)-g(t,y^{\prime}) |}{|y-y^{\prime}|^{\rho}},\\
[ g ]_{\rho,t}=& \underset{(t,y),(t^\prime,y) \in  [0,T] \times [0,1]}{\sup} \frac{| g(t,y)-g(t^\prime,y)|}{|t-t^{\prime}|^{\rho}}.
\end{align*}

\begin{align*}
C^{\frac{\ell+\alpha}{2k}, \ell+\alpha}\left([0,T] \times [0,1]\right),
~~~~\text{ for } \alpha \in (0,1) \text{ and } \ell \in \mathbb{N}_{0},
\end{align*}

represents the set of of all functions $g:  [0,T] \times [0,1] \to \mathbb{R}$,
which have continuous derivatives
$\partial_{y}^{\mu}\partial_{t}^{\nu}g$ for all $\nu, \mu \in \mathbb{N}_{0}$, $2\nu+\mu \leq \ell$,
and with the finite norm

\begin{align*}
\norm{g}_{C^{\frac{\ell+\alpha}{2k}, \ell+\alpha}\left( [0,T] \times [0,1]\right)}:
=&\sum\limits^{\ell}_{2k\nu+\mu=0} \underset{(t,y) \in  [0,T] \times [0,1]}{\sup} | \partial_{y}^{\mu}\partial_{t}^{\nu}g(t,y) |
+\sum\limits_{2k\nu+\mu=\ell}[\partial_{y}^{\mu}\partial_{t}^{\nu}g]_{\alpha, y}\\
\end{align*}

From the definition it follows that for $m\leq k$, $m,k \in \mathbb{N}_0$

$$  C^{ \frac{m+\alpha}{4},m+\alpha}([0,T]\times[0,1]) \subset  C^{ \frac{k+\alpha}{4}, k+\alpha}([0,T]\times[0,1]) \, ,$$

and if $v \in C^{ \frac{k+\alpha}{4}, k+\alpha}([0,T]\times[0,1])$, then $\pa_x^l v \in C^{ \frac{k-l+\alpha}{4}, k-l+\alpha}([0,T]\times[0,1])$ for all $0 \leq l\leq k$ so that

$$\|\pa_x^l v \|_{C^{ \frac{k-l+\alpha}{4}, k-l+\alpha}([0,T]\times[0,1])} \leq \| v\|_{C^{ \frac{k+\alpha}{4}, k+\alpha}([0,T]\times[0,1])} \, .$$

In particular at each fixed $x \in [0,1]$, we have $\pa_x^l v (\cdot,x )\in C^{s,\beta}([0,T])$
with $s =[\frac{k-l+\alpha}{4}]$ and $\beta = \frac{k-l+\alpha}{4}-s$.

For $k \in \mathbb{N}_0$, $\alpha,\beta \in (0,1)$ and $T>0$, we have

* if $v,w \in C^{\frac{k+\alpha}{4},k+\alpha}([0,T]\times [0,1])$, then

$$ \| v w \|_{C^{\frac{k+\alpha}{4},k+\alpha}} \leq C \| v \|_{C^{\frac{k+\alpha}{4},k+\alpha}} \| w \|_{C^{\frac{k+\alpha}{4},k+\alpha}} \, ,$$

with $C=C(k)>0$;

* if $v\in C^{\frac{\alpha}{4},\alpha}([0,T]\times [0,1])$, $v(t,x)\ne 0$ for all $(t,x)$, then

$$ \Big\| \frac{1}{v} \Big\|_{C^{\frac{\alpha}{4},\alpha}} \leq \Big\| \frac{1}{v} \Big\|^2_{C^{0}([0,T]\times [0,1])} \| v \|_{C^{\frac{\alpha}{4},\alpha}} \, .$$

Similar statements are true for functions in $C^{k,\beta}([0,T])$ and $C^{k,\beta}([0,1])$.

Let $T<1$ and $v \in C^{ \frac{4+\alpha}{4}, 4+\alpha}([0,T]\times[0,1])$ such that $v(0,x)= 0$, for any $x \in [0,1]$ then

$$ \| \pa_x^l v \|_{C^{\frac{m+\alpha}{4}, m+\alpha}} \leq C(m) T^{\beta} \|  v \|_{C^{\frac{4+\alpha}{4}, 4+\alpha}} $$

for all $l,m \in \mathbb{N}_0$ such that $l+m<4$. Here $\beta=\max \{ \frac{1-\alpha}{4}, \frac{\alpha}{4} \} \in (0,1)$; more precisely for $l\geq 1$ then $\beta=\frac{\alpha}{4}$.
In particular, for each $x \in [0,1]$ fixed

$$ \| \pa_x^l v (\cdot, x) \|_{C^{0,\frac{m+\alpha}{4}}([0,T])} \leq C(m) T^{\beta} \|  v \|_{C^{\frac{4+\alpha}{4}, 4+\alpha}} $$

for all $l,m \in \mathbb{N}_0$ such that $l+m<4$.

Let $0=x_0<x_1 <\cdots<x_{q}=1$, $I_{l}=(x_{l-1},x_{l})$, and $\Omega$ denote the set of all $C^{2k-2}$ paths $\gamma : [0,q] \to M$ such that
$\gamma_{l}:=\gamma_{|_{\bar{I}_{l}}}$ is smooth $C^{\infty}$. We define the tangent space of $\Omega$ at a path $\gamma$, $T_{\gamma}\Omega$,
as the set of all $C^{k-1}$ vector fields $W: [0,q] \to TM$ such that $W_{|_{[x_{l-1},x_{l}]}}$ is smooth.
If $\alpha : (-a,a) \times [0,q] \to M$ is a one-parameter variation of $\gamma \in \Omega$ is defined by

\begin{align*}
\alpha(\varepsilon,x)=exp_{\gamma}(\varepsilon W(x))
\end{align*}

and $W \in T_{\gamma}\Omega$ is
the variational vector field associated to $\alpha$, then

\begin{align*}
\nonumber
&\frac{d}{d\varepsilon} \big\vert_{\varepsilon=0} \mathcal{E}_{k}[\alpha(\varepsilon,x)]
\\ \nonumber
(-1)^{k} \sum^{q}_{l=1} \int_{I_{l}} \bigg{\langle}
D_{x}^{2k-1} \partial_{x} \gamma_{l}
%\\ \nonumber
+\sum^{k}_{\mu=2} (-1)^{\mu}
R\left(D^{2k-\mu-1}_{x} \partial_{x} \gamma_{l},
D^{\mu-2}_{x} \partial_{x} \gamma_{l}\right) \partial_{x} \gamma_{l} ,
W \bigg{\rangle} \text{ }dx
\\ \nonumber
&~~~+\sum^{k}_{\ell=1} \sum^{q}_{l=1}(-1)^{\ell-1}
\bigg{\langle} D^{k-\ell}_{x}W,
D^{k+\ell-2}_{x}\partial_{x} \gamma_{l} \bigg{\rangle} \bigg{|}^{x^{-}_{l}}_{x^{+}_{l-1}}
\\
&~~~+\sum^{k-1}_{\mu=2}\sum^{k-j}_{\ell=1}\sum^{q}_{l=1} (-1)^{\ell-1}  \bigg{\langle}
D^{k-\mu-\ell}_{x}R(W, \partial_{x} \gamma_{l}) D^{\mu-2}_{x}\partial_{x} \gamma_{l},
D^{k+\ell-2}_{x}\partial_{x} \gamma_{l}
\bigg{\rangle} \bigg{|}^{x^{-}_{l}}_{x^{+}_{l-1}},
\end{align*}

\begin{align*}
\frac{d}{d\varepsilon} \big\vert_{\varepsilon=0} \mathcal{T}[\alpha(\varepsilon,x)]=&
-\sum_{l=1}^{q}\int_{I_{l}} \langle D_{x}\partial_{x} \gamma_{l} , W\rangle \text{ }dx
+\sum_{l=1}^{q} \langle W,  \partial_{x} \gamma_{l} \rangle \bigg{\vert}^{x_{l}}_{x_{l-1}}.
\end{align*}

Let $a_1,\ldots, a_n$ be positive numbers and assume $0<i_1+\cdots +i_n<2$,
where $i_\mu>0$ for all $\mu=1, \ldots, n$.
Then, for any $\varepsilon>0$, we have

\begin{align*}
a^{i_1}_{1} \cdot a^{i_2}_{2}\cdot \ldots \cdot a^{i_n}_{n} \leq \varepsilon (a^2_1+\ldots +a^2_{n})+C_{\varepsilon}\text{.}
\end{align*}

For simplicity, we follow [7] and our previous papers to define the Sobolev norms,

\begin{equation*}
\| g\|_{k,p} := \sum_{\ell=0}^{k} \| D_x^{\ell} g \|_{p}
\quad \mbox{ with } \quad
\| D_x^{\ell} g\|_{p} :=  \Big( \int_I |D_x^{\ell} g |^p dx \Big)^{1/p} \,
\text{,}
\end{equation*}

which are not the same as the scale-invariant Sobolev norms in [7].

Let $\gamma: I \rightarrow M \subset \mathbb{R}^n$ be a smooth map, and $\phi$ be a smooth vector field on $M$.
Then, we have the followings.

\begin{align*}
(i.)~ & | \partial_{x}( |\phi|) | \leq  |D_{x}\phi| ~~~~~ \text{ almost everywhere},
\\
(ii.)~ &
\| D_{x}\phi \|_{2} \leq c(n) \left(\varepsilon  \| \phi \|_{2, 2}+\frac{1}{\varepsilon} \| \phi \|_{2} \right)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\forall\,  \varepsilon \in (0,1),
\\
(iii.)~ & \| D_{x}^{i}\phi \|_{2} \leq  c(i, n, k) \cdot (\varepsilon  \| \phi \|_{k, 2}
+\varepsilon^{\frac{i}{i-k}} \| \phi \|_{2}),
~~\forall\, \varepsilon \in (0,1), k\geq 2, 1\le i\le k-1,
\\
(iv.)~ & \| D_{x}^{i}\phi \|_{2} \leq  c(i, n, k)  \| \phi \|_{k, 2}^{\frac{i}{k}} \| \phi \|_{2}^{\frac{k-i}{k}},
~~~~~~~~~~~~~~~~~~~~~~~~~~\forall\, k\geq 2, 1\le i\le k-1,
\\
\nonumber
(v.)~ &  \| D_{x}^{i}\phi \|_{p} \leq
c(i, p, n, k) \cdot \| D_{x}^{i}\phi \|_{k-i, 2}^{\frac{1}{k-i}\left(\frac{1}{2}-\frac{1}{p}\right)} \,  \,  \| D_{x}^{i}\phi \|_{2}^{1-\frac{1}{k-i}\left(\frac{1}{2}-\frac{1}{p}\right)},
\\
& ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\forall\, p \geq 2, k \geq 1, 1\le i\le k-1,
\\
\nonumber
(vi.)~ & \| D_x^{i}\phi \|_{p} \leq  c(i, p, n, k) \cdot \|\phi \|_{k, 2}^{\alpha}  \cdot \|\phi \|_{2}^{1-\alpha},
\\
& ~~~~~~~~~~~~~~~~~~~~~~~~~~~ \forall\,  0 \leq i \leq k, k \in \mathbb{N}, p \geq 2,
~\text{ where }~ \alpha=\frac{1}{k}\left(i+\frac{1}{2}-\frac{1}{p}\right).
\end{align*}

The proof follows exactly the same argument in Lemmas C.3 $\sim$ C.6 and 4.1 in [7],
except that the covariant differentiation $D$ is
the Levi-Civita connection of a Riemannian manifold $M$
defined in (<ref>), i.e.,
$D_x \phi (\gamma)=\partial_x \phi (\gamma)-\Pi_\gamma (\phi, \partial_x\gamma)$.
Note that the covariant differentiation $D_x$ in an extrinsic setting is also a projection to the tangent subspace $T_{\gamma}M\subset\mathbb{R}^n$ after taking partial differentiation of a tangent vector field $\phi$ along $\gamma(x)$ in $M$, while the covariant differentiation $\nabla_s$ in [7] is defined as
a projection to the normal subspace after taking partial differentiation of a normal vector field $\phi$ w.r.t. the arclength parameter of $\gamma$.
Since the proof is the same as the one in [7], we leave the proof to the interested reader.

This project is supported by the following grants.
C.-C. L received support from the research grant
(MoST 111-2115-M-003-015) of the
National Science and Technology Council of Taiwan,
the iCAG program of Higher Education Sprout Project of National Taiwan Normal University and the Ministry of Education (MoE) of Taiwan,
and the National Center for Theoretical Sciences at Taipei;
T.D.T. received support from the research grant
(MoST 111-2115-M-003-014) of National Science and Technology Council of Taiwan.

Data availability.
Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

The authors have no conflict of interests related to this publication.

[3] M. Camarinha, F. Silva Leite, and P. Crouch,
Splines of class $C^k$ on non-Euclidean spaces.
IMA J. Math. Control Inform. 12 (1995), no. 4, 399–410.

[4] M. Camarinha, F. Silva Leite, and P. Crouch,
Riemannian cubics close to geodesics at the boundaries.
J. Geom. Mech. 14 (2022), no. 4, 545–558.

[5]
A. Dall'Acqua, C.-C. Lin, and P. Pozzi, Flow of elastic networks: long time existence results,
Geom. Flows. 4 (2019), no. 1, 83–136.

[6]
A. Dall'Acqua, C.-C. Lin, and P. Pozzi, Elastic flow of networks: short-time existence result.
J. Evolution Eq.
21 (2021), no. 2, 1299–1344.

[7] A. Dall'Acqua, and P. Pozzi,
A Willmore-Helfrich $L^2$-flow of curves with natural boundary conditions,
Comm. Anal. Geom.  22 (2014), no. 4,  617–669.

[8] S. A. Gabriel, and J. T. Kajiya,
Spline interpolation in curved space,
State of the Art in Image Synthesis, SIGGRAPH’85 Course Notes,
vol. 11, ACM Press, New York, 1985, pp. 1–14.

[9] F. Gay-Balmaz, D. D. Holm, D. M. Meier, T. S. Ratiu, and F.-X. Vialard, Invariant higher-order variational problems,
Commun. Math. Phys. 309 (2012), 413–458.

[10]
F. Gay-Balmaz, D. D. Holm, D. M. Meier, T. S. Ratiu, and F.-X. Vialard,
Invariant higher-order variational problems II,
J. Nonlinear Sci. 22 (2012), 553–597.

[11] R. Giambò, F. Giannoni, and P. Piccione,
An analytical theory for Riemannian cubic polynomials,
IMA J. Math. Control Inform. 19 (2002), 445–460.

[12] R. Giambò, F. Giannoni, and P. Piccione,
Optimal Control on Riemannian Manifolds by Interpolation,
Math. Control Signals Syst. 16 (2003), no. 4,  278–296.

[13]
B. Heeren, M. Rumpf, and B. Wirth,
Variational time discretization of Riemannian splines.
IMA J. Numer. Anal. 39 (2019), no. 1, 61–104.

[14] C.-C. Lin, Y.-K. Lue, and D. T. Tran,
A second-order elastic flow for path planning in $\mathbb{R}^2$.
Commun. Contemp. Math., (2023) Published online.

[15] C.-C. Lin, H. R. Schwetlick, and D. T. Tran,
An elastic flow for nonlinear spline interpolations in $\mathbb{R}^n$.
Trans. Amer. Math. Soc. 375 (2022), no. 7, 4893–4942.

[16] L. Machado, F. Silva Leite, and K. Krakowski,
Higher-order smoothing splines versus least squares problems on Riemannian manifolds.
J. Dyn. Control Syst. 16 (2010), no. 1, 121–148.

[17] L. Machado, T. Teresa, and M. Monteiro,
A numerical optimization approach to generate smoothing spherical splines.
J. Geom. Phys. 111 (2017), 71–81.

[19] L. Noakes,
Null cubics and Lie quadratics.
J. Math. Phys. 44 (2003), no. 3, 1436–1448.

[20]
L. Noakes, G. Heinzinger, and B. Paden,
Cubic splines on curved spaces.
IMA J. Math. Control Inform. 6 (1989), no. 4, 465–473.

[21] C. Samir, P.-A. Absil, A. Srivastava, and E. Klassen,
A gradient-descent method for curve fitting on Riemannian manifolds.
Found. Comput. Math., 12 (2012), no.1, 49–73.

[22] V. A. Solonnikov, Boundary Value Problems of Mathematical Physics. III,
Proceedings of the Steklov institute of Mathematics (1965), Amer. Math. Soc., Providence, R. I., No. 83, 1967.