text:
\section{Introduction}\label{sec:introduction}

Physiological flows are unsteady in nature and are characterized by complex geometries and fluid-structure interaction.
In healthy individuals, arterial flow is generally assumed laminar, but complex (\emph{disturbed}) flow patterns are acknowledged to play an important mechanistic role in the development of vascular diseases {\citep{Ku1997,Chiu2011}}. Even for the simple case of pulsatile flow in a straight pipe, the mechanisms of instability and transition to turbulence are poorly understood and particularly the dependence on the pulsation amplitude ($\A=U_o/U_s$, where $U_o$ and $U_s$ are the magnitude of the oscillatory and steady components of the velocity) is largely unknown. This makes it difficult to assess whether disturbed flow patterns in arterial and respiratory flow are solely due to geometric and structural effects (e.g.\ vessel curvature and flexibility, bifurcations), or are also related to the stability of pulsatile pipe flow.

Pulsatile flow in a straight pipe is governed by the pulsation amplitude $\A$, the {Womersley} number $\Wo=D/2\sqrt{\omega/\nu}$ and the Reynolds number $\Rey_s=U_s D/\nu$. Here $D$ the pipe diameter, $\omega$ the angular frequency of the pulsation and  $\nu$ the kinematic viscosity of the fluid. The limit of small pulsation amplitude $\A\rightarrow 0$ (steady flow), is relevant to laminar blood flow in capillaries, whereas the opposite limit $\A\rightarrow\infty$ (oscillatory flow) is relevant to respiratory flow.  In humans, the airflow may be laminar, transitional or turbulent depending on the airway segment \citep{kleinstreuer2010}. The intermediate regime, in which the pulsatile flow component is similar to the steady one ($\A\gtrsim 1$), is typical of blood flow in the large arteries.

Steady laminar pipe flow is linearly stable and  transition can only be triggered with finite-amplitude disturbances \citep{reynolds83}.  Following transition, turbulence persists provided that $\Rey_s\gtrsim2040$ \citep{Avila11}. Despite the non-linear nature of the transition, the key underlying mechanism is linear \citep{SchmidHenningson_springer2001}. The Navier--Stokes equations linearized about the laminar flow are non-normal and disturbances can be transiently amplified before asymptotically decaying. \citet{schmid1994} showed that in pipe flow the optimal (non-modal) disturbance consists of a pair of stream-wise rolls, which generate a pair of stream-wise velocity streaks (lift-up mechanism). The perturbation's energy gain in this process scales as $G\propto Re_s^2$, and keeping the flow laminar as $Re_s$ increases becomes an arduous task.

Oscillatory pipe flow is linearly (Floquet) unstable when the Reynolds number based on the Stokes-layer thickness, $\Rey_\delta=U_o\delta/\nu\gtrsim 10^3$, where  $\delta = \sqrt{2\nu/\omega} $ is the thickness of the Stokes layer \citep{Thomas12}. In experiments turbulence was observed already for $280 \lesssim \Rey_\delta \lesssim 550$ \citep{Sergeev66,Merkli75,Hino76,Eckmann91,zhao1996}, indicating that oscillatory pipe flow also undergoes transition via finite-amplitude disturbances. {\citet{Feldmann12,Feldmann16} performed direct numerical simulations of oscillatory pipe flow initialized with fully turbulent fields and confirmed the existence of sustained turbulence in the sub-critical regime.} \citet{Thomas11} showed that the linear instability of oscillatory pipe flow persists also for pulsatile flow, however in experiments transition occurs much earlier \citep[see e.g.][]{Sarpkaya66, Stettler86,Xu17,xu2020}. Taken together these results suggest that pulsatile pipe flow undergoes a sub-critical transition to turbulence in all regimes.

\citet{xu2020} recently reported on a nonlinear instability of pulsatile pipe flow, which occurs at pulsation amplitudes relevant for arterial flow. In their experiments, geometric imperfections triggered a helical wave pattern which emerged cyclically during the deceleration phase and broke down into turbulence, before decaying.  For $\A=1$ they observed transition at Reynolds numbers as low as $\Rey_s\approx 800$. In addition, \citet{xu2020} performed also linear non-modal transient growth computations at a selected parameter set and showed that the most amplified disturbance is a helical wave. Direct numerical simulations initialized with this helical wave reproduced the flow patterns and the time of turbulence breakdown observed experimentally, which suggests an important role of transient growth. {\citet{Nebauer2019} anticipated that in certain regimes of pulsatile pipe flow, the energy growth scales exponentially with the Reynolds number. This was also reported for the oscillatory Stokes layer over a flat plate \citep{Biau16} and for the flow following an axisymmetric stenosis \citep{Blackburn08JFM}.  Large non-modal transient amplification of disturbances has also been found in pulsatile channel flow \citep{Tsigklifis17}.} In what follows, we present a comprehensive non-modal linear analysis of the sub-critical regimes of pulsatile and oscillatory pipe flows.

\section{Methods}\label{sec:methods}

We consider an incompressible viscous fluid driven at a pulsatile flow rate in a straight pipe of circular cross-section. Lengths, velocities and time are rendered dimensionless with $D$,  $U_s$ and $D/U_s$, respectively. The dimensionless fluid velocity averaged over the circular cross-section reads
\begin{equation}
U(t) = U_s\cdot[1 + \A\cdot \textrm{sin}(2\pi \cdot t/T)],
\label{eq:Re}
\end{equation}
where $T=\pi \Rey_s/(2\Wo^2)$ is the dimensionless pulsation period. In oscillatory flow, we rendered lengths, velocities and time dimensionless with $D$, $U_o$ and $D/U_o$, respectively. {The  Navier--Stokes equations linearized about the Sexl-Womersely solution \citep{Womersley55}, and the corresponding adjoint equations,} were discretized in cylindrical coordinates $(r, \theta, z)$ with a Chebyshev collocation method for each Fourier mode $(k, m)$, where $k$ and $m$ are the axial and azimuthal wavenumbers of the perturbation, respectively. A second-order projection scheme was used to integrate the equations in time \citep[see][for details of the method]{xu2020}. {The code was validated against the Floquet analysis of \citet{Thomas12}}. {In this study, we used $N=96$ radial points and time step $\Delta t=0.002$} (convergence was checked at selected parameter values with $\Delta t=0.0005$ and  $N=128$). The optimal transient energy growth of a disturbance $\mathbf{u'}_{km}$ with wavenumbers $(k, m)$ was computed as
	\begin{equation}
	G_{km}(t_0,\tau) = \max\limits_{||\mathbf{u'}_{km}(t_0)||_2\neq0} \; \frac{E_{km}(t_0+\tau)}{E_{km}(t_0)},
	\label{eq:G}
	\end{equation}
with the adjoint method of \citet{Barkley08}. Here $E_{km}(t)$ is the kinetic energy of the disturbance at time $t$,  $t_0$ the time (phase) at which the perturbation is applied and $t_f$ the point at which the growth is evaluated ($\tau=t_f-t_0$ is the perturbation evolution time).
We varied $\Rey_s$, $\A$, and $\Wo$ independently and for each set of parameter values the maximum transient growth $G_\text{max}$ was optimized  over $t_0$, $\tau$, $k$ and $m$. We found that in most regimes the optimal azimuthal wavenumber is $m=1$ (except for some regimes of oscillatory flow). In all the results shown below  $m=1$ unless otherwise specified.

\section{Dynamics of the optimal disturbance}

\begin{figure}
	\centering
	\includegraphics[width=1\textwidth, trim={0cm 0cm 0cm 0cm}, clip]{figure1abc.eps}\\
	\includegraphics[width=1\textwidth, trim={0cm 0cm 0cm 0cm}, clip]{figure1defg.eps}\\
	\caption{\label{fig:Re2000Wo15A1_OptimalPerturbation} (a) Time series of the kinetic energy $E(t)/E(t_0)$ of the optimal helical, $(k,m)=(3.24,1)$, and classic, $(k,m)=(0,1)$, disturbances at $(Re_s,\A,\Wo)=(2000,1,15)$. (b)--(c) Time series of kinetic energy contribution of each velocity component  for the classic (b) and helical (c) disturbances.  (d)--(e) Contours of stream-wise vorticity (on a $r$-$\theta$ cross-section) and of stream-wise velocity (on a $z$-$r$ cross-section) of the classic disturbance, and the corresponding base flow profile ${\bm U}_{b,z}$ at $t_0/T=0.25$ (d) and $t_f/T=1.75$; see supplementary movie 1 for an animation of the disturbance dynamics  (e). (f)--(g) The same as (d)--(e), but for the helical disturbance at $t_0/T=0.5$ and $t_f/T=1.2$; see also supplementary movie 2. The dashed line denotes the Stokes-layer thickness.}
\end{figure}

The temporal evolution of the optimal disturbance's energy at $\Rey_s=2000$, $\A=1$ and $\Wo=15$  is shown as a dashed line in figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(a). At these parameter values the optimal disturbance has a helical structure with $(k,m)=(3.24,1)$ and is localized at the outer half of the pipe (exceeding the Stokes layer thickness), see figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(f)--(g). The optimal point to disturb is during the deceleration phase, at $t_0/T=0.5$, whereas the maximum amplification is reached during the acceleration phase, at $t_f/T=1.2$. The classic $k=0$, $m=1$  optimal disturbance of steady pipe flow is also amplified significantly in this case, albeit an order of magnitude less than the helical one. The classic disturbance initially consists of stream-wise vortices, as shown in figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(d)--(e), and the energy is subsequently transferred to the stream-wise velocity components, while the cross-stream components decay monotonically, see figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(b) and supplementary movie 1. Overall the classic perturbation's behavior appears to be rather insensitive to the change in flow profile throughout the cycle and the decay is very slow.

The kinetic energy of the optimal helical perturbation is mostly distributed in the stream-wise and azimuthal components, which self-amplify rapidly during the deceleration phase and a bit slower during the subsequent acceleration phase, see figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(c). Initially the disturbance spirals clock-wisely towards the pipe center while leaning against the background shear profile, see figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(f) and supplementary movie 2. As the energy grows, the perturbation switches the spiraling direction and is tilted by the shear until it aligns with it and the disturbance finally decays, see figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(g). This is reminiscent of the Orr mechanism \citep{Orr1907}. However, approximately $96\%$ of the kinetic energy is shared in equal parts between the azimuthal and stream-wise components, which indicates a strong three-dimensional effect, distinct from the two-dimensional Orr mechanism reported for many flows \citep[see e.g.][]{boyd1983,Farrell88,Maretzke2014, Biau16}.

{\citet{vonkerczek1982,Tsigklifis17} studied the stability of pulsatile channel flows and reported large \emph{modal} growth of disturbances at low $\Wo$, which also sets in during the deceleration phase and reaches the maximum at the end of the period of pulsation. These authors associated the modal energy growth with the inflectional velocity profiles occurring in the deceleration phase. Here, we found that at the optimal perturbation point $t_0$, the velocity profile has also inflection points, see figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(f).  We froze this flow profile and did a modal (eigenvalue) analysis to obtain the most dangerous mode. This mode was subsequently used as initial condition for the linearized Navier--Stokes equations and exhibited large energy growth (albeit an order of magnitude lower than the non-modal optimal). We found similar results in additional runs at $\Wo=5$ and $8$, suggesting that inflection points and modal mechanisms may also play an important role in pulsatile pipe flow.}

\begin{figure}
	\centering
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 3cm 0cm}, clip]{figure2a.eps}
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 3cm 0cm}, clip]{figure2b.eps}
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 3cm 0cm}, clip]{figure2c.eps}
	\caption{\label{fig:Re2000Wo15A1} Transient growth at $(Re_s,\A,\Wo)=(2000,1,15)$: (a) The red dashed lines denote the temporal evolution of the optimal classic disturbances $(k,m)=(0,1)$ for four different initial times $t_0$, whereas the blue solid lines correspond to the optimal helical disturbances $(k,m)=(3.24,1)$ initialized at the same $t_0$. The thick black line is the maximum gain $G(t_f)$ that can be achieved at a given time $t_f$ (optimized over $k$, $m$ and $t_0$ disturbances). (b) Dependence of the optimal axial wavenumber $k$ (associated to the thick line of a) on $t_f$. (c) Colormap of the maximum gain $G(t_0,\tau)$ (optimized over $k$ and $m$) in the $t_0-\tau$ plane. The black cross marks the maximum of $G$.}
\end{figure}

The black thick line in figure~\ref{fig:Re2000Wo15A1}(a) depicts the maximum energy amplification $G(t_f)$ over all $(k,m)$ and initial disturbances time $t=t_0$. The maximum amplification is reached during the acceleration phase via helical disturbances, as shown in figure~\ref{fig:Re2000Wo15A1}(b), whereas the classic disturbance achieves larger growth only during the second half of the deceleration phase. The colormap of figure~\ref{fig:Re2000Wo15A1}(c) shows that the optimal time to perturb the flow is during middle of the deceleration phase ($t_0/T\approx0.5$); perturbing during the acceleration phase leads to much lower growth (yielded by the classic disturbance during the deceleration phase). Clearly the helical mechanism is efficient only in the deceleration phase. The time needed to reach the maximum growth is about $\tau\approx 3/4T$, which explains why the maximum growth occurs during the acceleration phase.

\section{Parametric study of transient growth}

\begin{figure}
	\centering
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 1.5cm 0cm}, clip]{figure3a.eps}
	\includegraphics[width=0.32\textwidth, trim={-0.5cm 0cm 2cm 0cm}, clip]{figure3b.eps}
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 1.5cm 0cm}, clip]{figure3c.eps}
	\caption{\label{fig:Re_dependence} Optimal transient growth as a function of $\Rey_s$ for $(\A,\Wo)=(1,15)$. (a) The symbols show the computed maximum energy gain. The black solid line $G_\text{max}=\exp(\Rey_s/232-2.12)$ is a fit to the data for $\Rey_s\geq4000$. The orange dotted line $G_\text{max}=2.72\,\Rey_s^2$ is a fit to data for $(\A,\Wo)=(1,25)$. (b)--(c) Optimal initial and final disturbance times and axial wavenumber.}
\end{figure}

We show in figure~\ref{fig:Re_dependence}(a) that keeping $(\A,\Wo)=(1,15)$ fixed and increasing the Reynolds number leads to an approximately exponential increase of the energy gain $G_\text{max}$. Independently of the Reynolds number, the maximum energy amplification occurs always during flow acceleration with helical disturbances introduced during the deceleration phase, see figure~\ref{fig:Re_dependence}(b). As shown in figure~\ref{fig:Re_dependence}(c) an asymptotic behavior is approached for $\Rey_s\gtrsim 5000$, with $k\approx3.4$, $t_f/T\approx1.15$ and $t_0/T\approx0.35$.

\begin{figure}
	\centering
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 1.5cm 0cm}, clip]{figure4a.eps}
	\includegraphics[width=0.32\textwidth, trim={-0.5cm 0cm 2cm 0cm}, clip]{figure4b.eps}
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 1.5cm 0cm}, clip]{figure4c.eps}
	\caption{\label{fig:amplitude_dependence} Optimal transient growth as a function of $\A$ for $(\Rey_s,\Wo)=(2000,15)$. (a) The black circles show the computed maximum energy gain, the green squares and the orange diamonds are for $(\Rey_s,\Wo)=(2000,10)$ and $(\Rey_s,\Wo)=(2000,20)$, respectively. (b)--(c) Optimal initial and final disturbance times and axial wavenumber. The maximum gain as a function of $t_f/T$ is shown in the inset in (b)  for $\A=2.4$ and $2.6$ and illustrates the competition between two distinct different disturbances. }
\end{figure}

The effect of the pulsation amplitude $\A$, whilst keeping  $(\Rey_s,\Wo)=(2000,15)$ fixed is illustrated in figure~\ref{fig:amplitude_dependence}.  For steady pipe flow ($\A=0$), \citet{Meseguer03} obtained $G_\text{max}\approx288.7$, $\tau\approx 24.5$ and $(k,m)=(0,1)$, which is in line with our results were for pulsatile flow of low amplitude $A\lesssim0.55$. Helical disturbances dominate thereafter {(see the black line in figure~\ref{fig:amplitude_dependence}a)}, but their energy gain does not grow monotonically (exponentially) with the amplitude until $\A\gtrsim4$. For intermediate amplitudes there is competition between two distinct types of helical disturbances, as shown in figures~\ref{fig:amplitude_dependence}(b). For $0.55\lesssim\A\lesssim2.5$, the dominant helical disturbance is similar to that examined in detail in figures~\ref{fig:Re2000Wo15A1}--\ref{fig:Re_dependence}, whereas for $\A\gtrsim 2.5$ another helical disturbance, with much shorter axial length and shorter evolution time $\tau=t_f-t_0$, takes over. Both types of helical disturbances are triggered in the deceleration phase. Similar results were obtained for lower frequency, at $(\Rey_s,\Wo)=(2000,10)$, see the green lines in figures~\ref{fig:amplitude_dependence}(a) and (c). Here the maximum growth $G$ is over two orders of magnitude larger and helical disturbances dominate earlier (already for $A\gtrsim 0.41$). From the data in figures~\ref{fig:Re_dependence}(a) and \ref{fig:amplitude_dependence}(a) we conclude that the maximum energy gain of helical disturbances scales exponentially with the Reynolds number of the oscillatory flow component $\Rey_o = U_o D /\nu=A \Rey_s$.

\begin{figure}
	\centering
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 1.5cm 0cm}, clip]{figure5a.eps}
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 1.5cm 0cm}, clip]{figure5b.eps}
	\includegraphics[width=0.32\textwidth, trim={0cm 0cm 1.5cm 0cm}, clip]{figure5c.eps}
	\caption{\label{fig:frequency_dependence} Effect of the {Womersley} number $\Wo$ on  the transient growth. (a) Maximum energy gain for $(\Rey_s,\A)=(2000,1)$ (black circles), $(\Rey_s,\A)=(2000,4)$ (green squares) and oscillatory flow at $\Rey_o=8000$ $(\Rey_s,\A)=(0,\infty)$ (orange diamonds). The hollow symbols show the optimal growth of classic disturbances $(k,m)=(0,1)$ in the regime where the helical disturbance domiates. (b) Optimal initial and final disturbance times at $(\Rey_s,\A)=(2000,1)$. The pink dashed line shows $t_f/T = (\tau_\text{steady} + t_0)/T$, where $\tau_\text{steady} $ is the optimal evolution time of the classic disturbance in steady pipe flow at $\Rey_s=2000$. {The inset shows a zoom for $\Wo\leqslant18$.} (c) Optimal axial wavenumber.}
\end{figure}

In what follows we  examine the influence of the pulsation frequency on the dominant mechanism of transient growth in more detail. We begin by focusing on $\Rey_s=2000$ and $\A=1$. The circles in figure~\ref{fig:frequency_dependence}(a) show that  for sufficiently large $\Wo\gtrsim 20$ the classic disturbance dominates and the optimal gain $G$ of steady pipe flow is recovered (exactly as for the data shown as a orange lines in figures~\ref{fig:Re_dependence} and \ref{fig:amplitude_dependence}). As shown in figure~\ref{fig:frequency_dependence}(b), the perturbation growth time $\tau\approx25$ is also indistinguishable from that of steady pipe flow.  It can be concluded that in the limit $\Wo\rightarrow \infty$ the disturbance response is solely governed by $\Rey_s$. The dynamics of the steady pipe flow is also recovered in the quasi-steady limit $\Wo\rightarrow 0$, where the classic perturbation dominates  as well. In this limit, the maximum transient growth is governed by the maximum Reynolds number $\Rey_\text{max} = (1+\A)\Rey_s=\Rey_s+\Rey_o$ and occurs for $t_0,t_f\rightarrow T/4$, see the inset in figure~\ref{fig:frequency_dependence}(b).  For instance, for $\Rey_s=2000$, $\A=1$ and $\Wo=3$ the maximum gain is $G\approx1100$, which is slightly less than the gain for steady pipe flow $G\approx1155$ at $\Rey_s=4000$. As shown in figure~\ref{fig:frequency_dependence}(c), helical disturbances dominate for intermediate $4\lesssim\Wo\lesssim18$ and exhibit a sharp maximum in gain at $\Wo\approx 7$, figure~\ref{fig:frequency_dependence}(a). The results for $\A=4$ (shown as green squares) are qualitatively similar to those for $\A=1$, but the energy gain is much larger.

\section{Transient growth in oscillatory pipe flow} \label{sec:oscillation}

\begin{figure}
	\centering
	\includegraphics[width=0.32\textwidth, trim={0.5cm 0cm 0cm 0cm}, clip]{figure6a.eps}
	\includegraphics[width=0.32\textwidth, trim={0.5cm 0cm 0cm 0cm}, clip]{figure6b.eps}
	\includegraphics[width=0.32\textwidth, trim={0.5cm 0cm 0cm 0cm}, clip]{figure6c.eps}
	\includegraphics[width=1\textwidth, trim={-0.2cm 0cm 0.2cm 0cm}, clip]{figure6d2k.eps}
	\caption{\label{fig:OPF} {Transient growth of disturbances in oscillatory pipe flow for $\Wo=10$, $15$ and $20$ as a function of $\Rey_{\delta}$: (a) energy amplification, (b)--(c) optimal initial and final disturbance times and optimal wavenumbers, where the axial wavenumber is scaled with $\delta^{-1}$. (d)--(g) Contours of span-wise vorticity (on a $z-r$ cross-section) of the optimal helical disturbance ($m=1$) for $\Wo=10$ and $\Rey_\delta=530$ at $t_0/T=0.35$ (d), $t/T= 0.4$ (e), $t/T=0.5$ (f) and $t_f/T=0.75$ (g); see supplementary movie 3 for an animation. (h)--(k) The same as (d)--(g), but for the optimal axisymmetric disturbance ($m=0$) for $\Wo=15$ and $\Rey_\delta=589$ at $t_0/T=0.35$ (h), $t/T= 0.4$ (i), $t/T=0.5$ (j) and $t_f/T=0.75$ (k); see supplementary movie 4.}
	}
\end{figure}

The dependence of the energy gain $G$ on the frequency for the specific case of oscillatory flow  at $\Rey_o=8000$ is shown as orange diamonds in figure~\ref{fig:frequency_dependence}(a) and follows the trend of pulsatile flow. In figure~\ref{fig:OPF}(a) these data are shown as a function of the Reynolds number of the Stokes layer $\Rey_{\delta}=U_o\delta/\nu=\Rey_o/(\sqrt{2}\Wo)$, together with three additional sets for $\Wo=10$, $15$ and $20$ covering wide ranges of $\Rey_o$. At  low $\Rey_{\delta}$, the maximum gain $G=1$ is reached for $t_0=t_f$, implying that all perturbations decay monotonously, and the disturbance with ($k=0, m=0$) is the least damped, see figure~\ref{fig:OPF}(b)--(c). This can also be seen in figure~\ref{fig:frequency_dependence}(a), where the  limit $G\rightarrow 1$ is reached by the orange data points as {$\Wo\rightarrow\infty$} (and hence $\Rey_{\delta}\rightarrow0$ for constant $\Rey_o=8000$). For $\Rey_{\delta} \gtrsim 100$ the energy gain increases exponentially with $\Rey_{\delta}$ and depends only weakly on the {Womersley} number. All data sets shown in figure~\ref{fig:OPF}(a) exhibit excellent collapse, except for the data set at $\Rey_o=8000$, which deviates from the exponential scaling for $\Rey_\delta \gtrsim 630$ (corresponding to $\Wo\lesssim9$, i.e.\ near and below the peak in figure~\ref{fig:frequency_dependence}a). The optimal point to disturb the flow is during the deceleration phase at $t_0/T\approx0.35$ (as in figure~\ref{fig:Re_dependence}b at sufficiently high $\Rey_o$). The maximum energy gain is attained toward the end of the deceleration phase at low $\Rey_{\delta}$, and moves progressively into the acceleration phase as $\Rey_{\delta}$ increases.

Snapshots of the span-wise vorticity of the optimal perturbation are shown in figure~\ref{fig:OPF}(d)--(g) for $\Wo=10$ and $\Rey_\delta=530$. Initially, the perturbation leans against the background shear and then tilts to align with the stream-wise direction as it propagates radially inwards (see supplementary movie 3). The helical structure and dynamics of the disturbance are very similar to that shown in figure~\ref{fig:Re2000Wo15A1_OptimalPerturbation}(f), but here the disturbance is initially more confined to the Stokes layer. In figure~\ref{fig:OPF}(h)--(k) and supplementary movie 4 we show that at higher $\Wo=15$, the optimal disturbance is very similar in dynamics, but is axisymmetric (see also figure~\ref{fig:OPF}(c)). The energy growth occurring here is as in the classic Orr mechanism, with the only difference that the perturbation travels radially inwards as it changes the tilt direction. We stress that despite this difference in the disturbance geometry (helical versus axisymmetric), the energy gain, the point of disturbance and the point of maximum gain are indistinguishable (also for larger $\Wo=20$, see figure~\ref{fig:OPF}a--c).

\section{Conclusion}\label{sec:conclusion}

We showed that the classic lift-up mechanism produces the largest transient growth in pulsatile flow of low amplitude, $A\lesssim0.4$. This is in agreement with recent experiments of \citet{Xu17} and direct numerical simulations of \citet{Xu18}, exhibiting turbulent puff and slugs as in steady pipe flow. At higher amplitudes, helical disturbances begin to dominate in a band of intermediate {Womersley} number $4\lesssim\Wo\lesssim18$ which progressively widens toward larger $\Wo$ as the oscillatory Reynolds number $\Rey_o$ is increased. For the specific case of oscillatory flow ($\A\rightarrow\infty$), our results are in qualitative agreement with the transient growth analysis of \citet{Biau16} for oscillatory Stokes flows over a flat plate. A key difference between oscillatory pipe flow and oscillatory flow over a flat plate is that the former is not solely governed by $\Rey_\delta$, but also by $\Wo$. In particular, the thickness of the Stokes layer in oscillatory pipe flow scales as $\Wo^{-1}$ and at sufficiently low $\Wo$ it fills the pipe. Hence, oscillatory pipe flow is only exactly similar to oscillatory flow over a flat plate in the limit $\Wo\rightarrow \infty$, where curvature effects become negligible \citep{Thomas12}. This convergence can be observed in figure~\ref{fig:OPF}(c); axisymmetric disturbances dominate for $\Rey_\delta\gtrsim200$ at $\Wo=20$, $\Rey_\delta\gtrsim400$ at $\Wo=15$, whereas for $\Wo=10$ helical disturbances dominate even up to  $\Rey_\delta\approx 800$. Note also that the axial wavenumber of the optimal disturbance scaled with $\delta^{-1}$ is close to $0.4$ in figure~\ref{fig:OPF}, which is in excellent agreement with figure~2 of \citet{Biau16}.

We here considered linear transient growth of disturbances, but transition to turbulence can only be completed with nonlinear effects. \citet{xu2020} showed that initializing direct numerical simulations with the linear optimal helical disturbance can trigger turbulent flow patterns as those observed in their experiments. In experiments of oscillatory pipe flow, the transition threshold is independent of $\Wo$ \citep[with different critical numbers $280 \lesssim \Rey_\delta \lesssim 550$ depending on the setup, see][]{Sergeev66,Merkli75,Hino76,Eckmann91,zhao1996}, which is consistent with our observation that the maximum energy gain depends only on $\Rey_\delta$. {We found large energy growth already at moderately low $\Rey_o = \mathcal{O}(1000)$, which suggests that non-modal disturbance energy growth may be at the root of disturbed flow patterns observed in physiology. Testing this hypothesis requires detailed direct numerical simulation of flows at physiologically relevant parameter regimes and the consideration of non-harmonic waveforms of the flow rate. We anticipate that the steep deceleration phase in aortic flow may be even more prone to instability than in the case of harmonic pulsation.}

\section*{Acknowledgment}
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the research
unit FOR 2688 `Instabilities, Bifurcations and Migration in Pulsatile Flows' under grant AV 120/6-1. Baofang Song acknowledges the financial support from the National Natural Science Foundation of China under grant number 91852105.

\bibliographystyle{jfm}


text:
\section{Introduction}\label{intro}

	In this paper we will establish a local null controllability result for the degenerate system, with nonlocal terms, given by
	\begin{equation}
	\begin{cases}
	u_t-\left(\mu_1 \left(x,\int_0^1u \   \right)u_x \right)_x+f_1(t,x,u,v) =h\chi_\omega, & (t,x)\in (0,T)\times (0,1),\\
	v_t-\left(\mu_2 \left(x,\int_0^1v \   \right)v_x \right)_x+f_2(t,x,u,v)=0, & (t,x) \in (0,T)\times (0,1),\\
	u(t,0)=u(t,1)=v(t,0)=v(t,1)=0, &  t\in (0,T),\\
	u(0,x)=u_0(x) \ \ \mbox{ and }\ \ v(0,x)=v_0(x), & x\in  (0,1),
	\end{cases}
	\label{pb1}
	\end{equation}
	where $T>0$ is given and  $(u_0,v_0)$ is the initial data. Moreover, $h$ is the control function, $(u,v)$ is the associated state and $\chi_{\omega}$ represent the characteristic function of     $\omega=(\alpha,\beta)\subset\subset (0,1)$. Regarding the functions $\mu_1 ,\mu_2, f_1$ and $f_2$, we make the following assumptions:

	\begin{assumption}{1} \label{hyp_a}  Let $\ell_1,\ell_2:\R\to \R$ be  $C^1$ functions with bounded derivative and suppose that $\ell_i(0)>0$, for each $i=1,2$. We also consider $a \in C([0,1])\cap C^1((0,1])$ satisfying $a(0)=0$, $a>0$ in $(0,1]$, $a'\geq 0$ and  \begin{equation}\label{prop_a}
		xa'(x)\leq Ka(x),\ \ \forall x\in [0,1] \mbox{ and some } K\in [0,1).
		\end{equation}
		In other words, the function $a$ behaves $x^{\alpha}$, with $\alpha \in (0,1)$.

		Throughout this article, we will consider the functions $\mu_1,\mu_2:[0,1]\times \R\to \R$, given by
		\[	\mu_1 (x,r)=\ell _{1} (r)a(x)\ \ \ \mbox{ and } \ \ \ \mu_2 (x,r)=\ell_2(r)a(x).\]
		And, for the sake of simplicity, we will also consider $\ell_i(0)=1$,   for each $i=1,2$.
	\end{assumption}

	\begin{assumption}{2}\label{hyp_f} For $i\in \{ 1,2\}$, we suppose that $f_i:[0,T]\times [0,1]\times \R^2\to \R$ is a $C^1$ function, with bounded derivatives, such that $f_i(t,x,0,0)=0$.  We also consider $b_{ij}(t,x):=\partial_{j+2}f_i(t,x,0,0)\in  L^\infty((0,T)\times (0,1))$, for any $i,j\in \{1,2\}$. And, let us assume that there exists $\omega_1\subset\subset \omega$ such that
		\[\inf \{b_{21}(t,x);\ (t,x)\in [0,T]\times \omega_1 \}>0.\]
	\end{assumption}

	The main purpose of this work is to prove the local null controllability of ($\ref{pb1}$) by means of one control. Precisely,
	we will obtain $h\in L^2((0,T)\times (0,1))$ such that the associated state $(u,v)=(u(t,x),v(t,x))$ of \eqref{pb1} satisfies
	\begin{equation*}
	u(T,x)\equiv v(T,x)\equiv 0 \mbox{ for any } x\in [0,1],
	\end{equation*}
	at least if $\|(u_0,v_0)\|_{H_a^1 \times H_a^1}$ is sufficiently small, where $H_a^1$ is an appropriate weighted Hilbert space which will be defined later in Section \ref{sec-pre} .

	The main difficulty comes from  the fact that the diffusion coefficients degenerate at $x=0$ and have nonlocal terms, namely
	\begin{equation*}
	\left(\mu_1 \left(\cdot ,\int_0^1u \   \right)u_x \right)_x \text{  and } \left(\mu_2 \left(\cdot,\int_0^1u \   \right)u_x \right)_x
	\end{equation*}
	satisfy assumption \eqref{hyp_a}.

	It is important to remark that semilinear nondegenerate problems have been extensively studied over the last decades, see \cite{fattorini1971,fernandez2012,fernandez2000cost,fursikov1996,lebeau1995} for example.

	However, it seems to us that there is also a large interest in degenerate operators when the degeneracy occurs at the boundary of the space domain. For instance, in \cite{oleinik1999mathematical}, it was developed a study about the Prandtl system for stationary flows, in which the related boundary layer system was reduced to a quasilinear degenerate parabolic equation. Degenerate operators also appear in probabilistic models, see \cite{feller1952parabolic,feller1954diffusion}, and in climate science, see \cite{floridia2014approximate}.

	In the context of degenerated  systems, controllability was studied in the case of two coupled equations in \cite{cannarsa2009controllability,hassi2011null,hassi2011carleman}. Recently, Ait Benhassi et al., in  \cite{benhassi2017algebraic},
	generalize the Kalman rank condition for the null controllability to $n$-coupled linear degenerate parabolic systems with  $m$-controls.

	On the other hand, as it was pointed out in \cite{fernandez2012}, nonlocal terms type can be found in several natural phenomena, such as in the reaction-diffusion systems, see \cite{chang2003some}, and in nonlinear vibration theory, see \cite{medeiros2002vibrations} for example.

	In \cite{alabau2006carleman},  it was obtained the null controllability for the semilinear equation
	\begin{align}\label{cf}
	u_t-(a(x)u_x)_x+f(t,x,u)=h(t,x)\chi_{\omega }, \text{ where }(t,x)\in (0,T)\times (0,1).
	\end{align}

	Based on this work, in \cite{jrl2016}, we have considered (\ref{cf}), replacing the second-order term $(au_x)_x$ by a specific degenerate nonlocal operator. In that new context, we have achieved a local null controllability result. For systems of parabolic equations, the
	main issue is often to reduce the number of control functions acting on the system (see \cite{benhassi2017algebraic,carreno2013local,coron2009null,fadili2017null}, for example), besides that,  as it was pointed out in \cite{ammar2011recent}, the problem of controlling coupled parabolic equations has a very different behavior with respect to the scalar case, for instance, boundary controllability is not equivalent to distributed controllability, approximate controllability is not equivalent to null controllability, and
	``\textit{the list of open
		problems is long and there is a lot of work to be done in order to fully understand
		this challenging subject}''\cite{ammar2011recent}.  In this direction, the current work may be seen as a natural continuation of \cite{jrl2016} and a first step in order to understand parabolic system with nonlocal and degenerate diffusion coefficients of the type $\left(\mu \left(\cdot ,\int_0^1u \   \right)u_x \right)_x $.

	Our main result is the following:

	\begin{thm}\label{main}
		Under the assumptions on $\mu_1,\mu_2,f_1$ and $f_2$, the nonlinear system \eqref{pb1} is locally null-controllable at any time $T>0$, i.e., there exists $r>0$ such that, whenever $u_0,v_0\in H_a^1$ and $\n{(u_0,v_0)}{H_a^1}\leq r$, there exists a control $h\in L^2((0,T)\times \omega)$ associated to a state $(u,v)$ satisfying
		\begin{equation*}
		u(T,x)=v(x,T)=0, \text{ for every } x\in [0,1].
		\end{equation*}
	\end{thm}

	The proof of Theorem \ref{main} will follow standard arguments (see for instance \cite{fernandez2012}, \cite{clark2013theoretical}, \cite{imanuvilov2001remarks}, \cite{fernandez2004local} and \cite{fernandez2006some}), based on Lyusternik's Inverse Mapping Theorem, which can be found in \cite{fursikov1996} and \cite{imanuvilov2003carleman}. To be more specific, we will see that the desired result is equivalent to find a solution to the equation
	\begin{equation}\label{eq}
	H(u,v,h)=(0,0,u_0,v_0),
	\end{equation}
	where $H:E\to F$ is a $C^1$ mapping between two appropriate Hilbert spaces, defined by
	\begin{equation*}
	H (u,v,h) = (H_{1} (u,v,h), H_{2} (u,v,h), u(\cdot ,0), v(\cdot ,0) ),
	\end{equation*}
	where
	\[\displaystyle H_{1} (u,v,h) = u_t-\left(\mu_1 \left(x,\int_0^1u \   \right)u_x \right)_x+f_1(t,x,u,v) -h\chi_\omega\]
	and
	\[\displaystyle H_{2} (u,v,h) = v_t-\left(\mu_2 \left(x,\int_0^1v \   \right)v_x \right)_x+f_2(t,x,u,v) . \]

	In order to use Lysternik's Theorem, we need to prove that $H'(0)$ is onto. It is equivalent to prove  a global null controllability result  to  the  linearization of \eqref{pb1} (see the system (\ref{pb-lin}) below). This approach relies on a suitable Carleman estimate for the solutions of the adjoint problem associated to \eqref{pb-lin} (see Proposition \ref{prop3.2}).

	This paper is organized as follows: Section \ref{sec-pre} contains notations we use and a preliminary result. In Section \ref{sec-carl}, we prove a Carleman type inequality to solutions of \eqref{adj}, which also allows us to obtain an Observality inequality. Section \ref{sec-NCLS} is concerned with the global null controllability of \eqref{pb-lin} as well as two crucial additional estimates. Section \ref{sec-main} is devoted to the proof of Theorem \ref{main}. In Section \ref{sec6}, we present some comments and remarks. At the end, we  provide three appendices where we sketch the proof of wellposedness of problem {\eqref{pb1}} and two minors results we use throughout the paper.
	\section{Preliminaries}\label{sec-pre}

	\noindent

	In this section we will  state some notations and results which are necessary to prove Theorem \ref{main}. At first, we need to introduce some weighted spaces related to the function $a$,  namely
			\begin{align*}
		& \begin{multlined}[t][0.9\textwidth]
		H_a^1:= \{  u\in L^2(0,1);\ u\mbox{ is absolutely continuous in } (0,1],\\
		\sqrt{a}u_x\in L^2(0,1) \mbox{ and } u(1)=u(0)=0\},
		\end{multlined}\\
		& H_a^2:= \{  u\in H_a^1;\ au_x\in H^1(0,1) \},
		\end{align*}
		with  the norms defined by\\
		$\|u\|_{H_a^1}^2:=\|u\|_{L^2(0,1)}^2+\|\sqrt{a}u_x\|_{L^2(0,1)}^2$
		and $\|u\|_{H_a^2}^2:=\|u\|_{H_a^1}^2+\|(au_x)_x\|_{L^2(0,1)}^2$.

		Alabau-Boussouira at al. in \cite{alabau2006carleman} introduced and studied some of the main properties of these spaces.

		As we pointed out in the Introduction, Lyusternik's Theorem requires the proof of a global null controllability result to a linear problem. Hence, we will present  a wellposedness result to this kind of problem that will be used later on the definition of our spaces in Section \ref{sec-main}. That result can be seen  in \cite{fadili2017null}.  Let us consider the linear system
		\begin{equation}
		\begin{cases}
		u_t-\left(a(x)u_x \right)_x+b_{11}(t,x)u+b_{12}(t,x)v =G_1, & (t,x) \mbox{ in }(0,T)\times (0,1),\\
		v_t-\left(a(x)v_x \right)_x+b_{21}(t,x)u+b_{22}(t,x)v =G_2, & (t,x) \mbox{ in }(0,T)\times (0,1),\\
		u(t,1)=u(t,0)=v(t,0)=v(t,1)=0, &  t \mbox{ in } (0,T),\\
		u(0,x)=u_0(x) \ \ \mbox{ and }\ \ v(0,x)=v_0(x), & x\mbox{ in } (0,1).
		\end{cases}
		\label{pb-lin1}
		\end{equation}

		\begin{prop}\label{prop-WP-lin}
			For all $G_1,G_2\in L^2((0,T)\times (0,1))$ and $u_0,v_0\in L^2(0,1)$, there exists a unique weak solution $u,v\in C^0([0,T];L^2(0,1))\cap L^2(0,T;H_a^1)$ of \eqref{pb-lin1} and there exists a positive constant $C_T$ such that
			\begin{multline*}
			\sup_{t\in[0,T]}\left(\n{u(t)}{L^2(0,1)}^2+\n{v(t)}{L^2(0,1)}^2\right) +\int_0^T\left(\n{\sqrt{a}u_x}{L^2(0,1)}^2+ \n{\sqrt{a}v_x}{L^2(0,1)}^2\right)\\
			\leq C_T\left(\n{u_0}{L^2(0,1)}^2+\n{v_0}{L^2(0,1)}^2 +\n{G_1}{L^2((0,T)\times (0,1))}+\n{G_2}{L^2((0,T)\times (0,1))} \right).
			\end{multline*}
			Moreover, if $u_0,v_0\in H_a^1$, then
			\[u,v\in H^1(0,T;L^2(0,1))\cap L^2(0,T;H_a^2)\cap C^0([0,T];H_a^1), \]
			and there exists a positive constant $C_T$ such that
			\begin{multline}\label{ineq1}
			\sup_{t\in[0,T]}\left(\n{u(t)}{H_a^1}^2+\n{v(t)}{H_a^1}^2\right)\\
			+\int_0^T\left(\n{u_t}{L^2(0,1)}^2+\n{v_t}{L^2(0,1)}^2+\n{(au_x)_x}{L^2(0,1)}^2+ \n{(av_x)_x}{L^2(0,1)}^2\right)\\
			\leq C_T\left(\n{u_0}{H_a^1}^2+\n{v_0}{H_a^1}^2 +\n{G_1}{L^2((0,T)\times (0,1))}+\n{G_2}{L^2((0,T)\times (0,1))} \right)
			\end{multline}
		\end{prop}

	Now, let us present a Carleman estimate that will play an important role in the next section. First, consider $\omega'=(\alpha',\beta')\subset\subset \omega$ and let $\psi:[0,1]\to \R$ be a $C^2 $  function such that
	\begin{align}\label{functions1}
	\psi(x):=\begin{cases}
	\phantom{-}\int_0^x \frac{y}{a(y)}dy,\ x\in [0,\alpha')\\
	-\int_{\beta'}^x \frac{y}{a(y)}dy,\ x\in [\beta',1].
	\end{cases}
	\end{align}
	For $\lambda$ sufficiently large, define
	\begin{multline}\label{functions}
	\theta(t):=\frac{1}{[t(T-t)]^4}, \ \eta(x):=e^{\lambda(|\psi|_\infty+\psi)},\ \sigma(x,t):=\theta(t)\eta(x) \mbox{ and }\\
	\varphi(x,t):=\theta(t)(e^{\lambda(|\psi|_\infty+\psi)}-e^{3\lambda|\psi|_\infty}).
	\end{multline}

	Let us also  consider the linear system
	\begin{equation}\label{adj-jlr}
	\left\{\begin{array}{ll}
	-\xi_t-\left(a\left(x\right)\xi_x \right)_x+c(t,x)\xi=F(t,x), & (t,x)\in (0,T)\times (0,1), \\
	\xi(t,1)=\xi(t,0)=0, & t\in (0,T),\\
	\xi(T,x)=\xi_T(x), & x\in(0,1),
	\end{array}\right.   \end{equation}
	where $a$ satisfies assumption \ref{hyp_a}, $c\in L^\infty((0,T)\times (0,1))$,  $F\in L^2((0,T)\times (0,1))$ and $\xi_T\in L^2(0,1)$.

	The following Carleman estimate, proved in \cite{jrl2016}, holds for the solution to \eqref{adj-jlr}:

	\begin{prop}\label{car-jlr}
		There exist $C>0$ and $\lambda_0,s_0>0$ such that every solution  $\xi$ of (\ref{adj-jlr}) satisfies, for all $s\geq s_0$ and $\lambda\geq \lambda_0$,
		\begin{equation}
		\int_0^T\int_0^1 e^{2s\varphi}\left((s\lambda)\sigma a\xi_x^2+(s\lambda)^2\sigma^2\xi^2 \right)
		\leq C\left(\int_0^T\int_0^1 e^{2s\varphi}|F|^2\   +(\lambda s)^3\int_0^T\int_{\omega} e^{2s\varphi}\sigma^3\xi^2\   \right),\label{carl_jlr}
		\end{equation}
		where the constants $C, \lambda_0,s_0$  only depends on $\omega$, $a$, $\n{c}{L^\infty((0,T)\times (0,1))}$ and $T$.
	\end{prop}

	For the sake of simplicity, we will  introduce the operator
	\begin{equation*}
	I(s,\lambda,\xi):=\I{\xi}.
	\end{equation*}

	Finally, let us present the following version of \textit{Lyusternik's Inverse Mapping Theorem} that can be found for instance in \cite{alekseev1987optimal,fursikov1996}.
\begin{thm}[Lyusternik]\label{lyusternik}
		Let $E$ and $F$ be two Banach spaces, $H: E\to F$ a $C^1$ mapping and put $\eta_0=H(0)$. If $H'(0):E\to F$ is onto, then there exist $r>0$ and $\tilde{H}:B_r(\eta_0)\subset F\to E$ such that
\[		H(\tilde{H}(\xi))=\xi,\ \forall \xi\in B_r(\eta_0),\]
		that is, $\tilde{H}$ is a right inverse of $H$. In addition, there exists $K>0$ such that
		\[\n{\tilde{H}(\xi)}{E}\leq K\n{\xi-\eta_0}{F}, \forall \xi\in B_r(\eta_0).  \]
	\end{thm}

	\section{Carleman and Observability inequalities} \label{sec-carl}

	In order to prove that map $H$ is onto, we have   to prove a global null controllability result  to  the  linearization of \eqref{pb1}, given by
	\begin{equation}
	\begin{cases}
	u_t-\left(a(x)u_x \right)_x+b_{11}(t,x)u+b_{12}(t,x)v =h\chi_\omega+g_1, & (t,x) \mbox{ in }(0,T)\times (0,1),\\
	v_t-\left(a(x)v_x \right)_x+b_{21}(t,x)u+b_{22}(t,x)v =g_2, & (t,x) \mbox{ in }(0,T)\times (0,1),\\
	u(t,1)=u(t,0)=v(t,0)=v(t,1)=0, &  t \mbox{ in } (0,T),\\
	u(0,x)=u_0(x) \ \ \mbox{ and }\ \ v(0,x)=v_0(x), & x\mbox{ in } (0,1),
	\end{cases}
	\label{pb-lin}
	\end{equation}
	where $g_1$, $g_2$ and $h$ belong to appropriate $L^2$-weighted spaces which we will specify later on.  To this purpose it is crucial to obtain an appropriate  Carleman estimate for solutions to
	\begin{equation}
	\begin{cases}
	-y_t-\left(a(x)y_x \right)_x+b_{11}(t,x)y+b_{21}(t,x)z =F_1, & (t,x) \mbox{ in }(0,T)\times (0,1),\\
	-z_t-\left(a(x)z_x \right)_x+b_{12}(t,x)y+b_{22}(t,x)z =F_2, & (t,x) \mbox{ in }(0,T)\times (0,1),\\
	y(t,0)=y(t,1)=z(t,0)=z(t,1)=0, &  t \mbox{ in } (0,T),\\
	y(T,x)=y_T(x),\ z(T,x)=z_T(x), & x\in (0,1),
	\end{cases}
	\label{adj}
	\end{equation}
	which is the adjoint problem of \eqref{pb-lin}.

	\begin{prop}\label{carl1}There exist positive constants $C,\lambda_0$ and $s_0$ such that, for any $s\geq s_0$, $\lambda\geq \lambda_0$ and any $y_{_T},z_{_T}\in L^2((0,T)\times (0,1))$, the corresponding solution $(y,z)$ to \eqref{adj} satisfies
		\begin{multline}
		\int_0^T\int_0^1 e^{2s\varphi}\left((s\lambda)\sigma a(y_x^2+z_x^2)+(s\lambda)^2\sigma^2(y^2+z^2) \right) \\
		\leq C\left(\int_0^T\int_0^1 e^{2s\varphi}s^4\lambda^4\sigma^4(|F_1|^2+|F_2|^2)\   +\int_0^T\int_{\omega} e^{2s\varphi}s^8\lambda^8\sigma^8y^2\   \right).\label{carl}
		\end{multline}
	\end{prop}

	\begin{proof}
		Firstly, we rewrite the first equation in \eqref{adj} as
		\[-y_t-\left(a(x)y_x \right)_x+b_{11}(t,x)y=F_1-b_{21}z,  (t,x) \mbox{ in }(0,T)\times (0,1).\]
		So, we apply Proposition \ref{car-jlr} for the case in which $\xi=y$, $F=F_1-b_{21}z$, $c=b_{11}$ and $\omega=\omega_1$ to obtain
		\begin{equation*}
		I(s,\lambda,y)
		\leq  C\bigg(\int_0^T\int_0^1 e^{2s\varphi}|F_1|^2+\n{b_{21}}{_{L^\infty}}\int_0^T\int_0^1 e^{2s\varphi}|z|^2 +\int_0^T\int_{\omega_1} e^{2s\varphi}\lambda^3s^3\sigma^3|y|^2 \bigg).
		\end{equation*}

		Proceeding in the same way for the second equation, we get an analogous inequality
		\begin{equation*}
		I(s,\lambda,z)
		\leq  C\bigg(\int_0^T\int_0^1 e^{2s\varphi}|F_2|^2+\n{b_{12}}{_{L^\infty}}\int_0^T\int_0^1 e^{2s\varphi}|y|^2 +\int_0^T\int_{\omega_1} e^{2s\varphi}\lambda^3s^3\sigma^3|z|^2 \bigg).
		\end{equation*}
		Now, we add this two inequalities and take  $s$ and $\lambda$ sufficiently large such that $I(s,\lambda,y)$ will  absorb the integral  depend on $|y|^2$, and $I(s,\lambda,z)$  the integral depend on $|z|^2$. This will give us the following inequality
		\begin{equation*}
		I(s,\lambda,y)+I(s,\lambda,z)\leq C\left(\int_0^T\int_0^1 e^{2s\varphi}(|F_1|^2+|F_2|^2)+\int_0^1\int_{\omega_1} e^{2s\varphi}\lambda^3s^3\sigma^3(|y|^2+|z|^2) \right).
		\end{equation*}

		Thus, in order to obtain \eqref{carl},  it is sufficient to show that there exists a small  $\varepsilon>0$  such that
		\begin{multline*}
		\int_0^1\int_{\omega_1}e^{2s\varphi}\lambda^3s^3\sigma^3|z|^2
		\leq \varepsilon I(s,\lambda,z)\\ + C\left(\int_0^T\int_0^1 e^{2s\varphi}s^4\lambda^4\sigma^4(|F_1|^2+|F_2|^2)+\int_0^1\int_{\omega} e^{2s\varphi}\lambda^8s^8\sigma^8|y|^2 \right).
		\end{multline*}

		Let us take $\chi\in C_0^\infty(\omega)$ satisfying $0\leq \chi\leq 1$ and $\chi\equiv 1$ in $\omega_1$. Since $\dps\inf b_{21}>0$, we can easily see that
		\[\int_0^1\int_{\omega_1}e^{2s\varphi}\lambda^3s^3\sigma^3|z|^2\leq C\int_0^T\int_{\omega} \chi b_{21}e^{2s\varphi}\lambda^3s^3\sigma^3|z|^2.\]

		Now, multiplying the first equation in \eqref{adj} by $e^{2s\varphi}s^3\lambda^3\sigma^3\chi z$ and integrating over $(0,T)\times (0,1)$, we get
		\begin{align*}
		\int_0^T\int_0^1 \chi b_{21} e^{2s\varphi}s^3\lambda^3\sigma^3 |z|^2  = &\int_0^T\int_0^1  \chi e^{2s\varphi}s^3\lambda^3\sigma^3 zF_1+\int_0^T\int_0^1 \chi e^{2s\varphi}s^3\lambda^3\sigma^3zy_t
		\\
		& +\int_0^T\int_0^1 \chi e^{2s\varphi}s^3\lambda^3\sigma^3(ay_x)_xz-\int_0^T\int_0^1 \chi b_{11} e^{2s\varphi}s^3\lambda^3\sigma^3 yz\\
		=: & I_1+I_2+I_3+I_4.
		\end{align*}

		Next, we need to estimate $I_1,I_2,I_3$ and $I_4$. Firstly, from Young's inequality, we have
		\begin{align}\label{I1}
		I_1\leq &  C\int_0^T\int_0^1 \chi  e^{2s\varphi}s^3\lambda^3\sigma^3 |z||F_1|
		\leq  \varepsilon \int_0^T\int_0^1 e^{2s\varphi}s^2\lambda^2\sigma^2 |z|^2+C_\varepsilon\int_0^T\int_0^1   \chi^2e^{2s\varphi}s^4\lambda^4\sigma^4 |F_1|^2\nonumber\\
		\leq & \varepsilon I(s,\lambda,z)+C_\varepsilon\int_0^T\int_0^1 e^{2s\varphi}s^4\lambda^4\sigma^4 |F_1|^2.
		\end{align}

		In the same way, since $b_{11}$ is bounded, it is immediate that
		\begin{align}\label{I4}
		I_4\leq &  C\int_0^T\int_0^1 \chi  e^{2s\varphi}s^3\lambda^3\sigma^3 |y||z|
		\leq  \varepsilon \int_0^T\int_0^1 e^{2s\varphi}s^2\lambda^2\sigma^2 |z|^2+C_\varepsilon\int_0^T\int_0^1   \chi^2e^{2s\varphi}s^4\lambda^4\sigma^4 |y|^2\nonumber\\
		\leq & \varepsilon I(s,\lambda,z)+C_\varepsilon\int_0^T\int_{\omega} e^{2s\varphi}s^8\lambda^8\sigma^8 |y|^2.
		\end{align}

		Using integration by parts,  we will split up $I_2$ and $I_3$ in several integrals. In fact,
		\begin{multline*}
		I_3= -\int_0^T\int_0^1 \chi s^3\lambda^3\sigma^3 e^{2s\varphi} az_xy_x+\int_0^T\int_0^1\left(a(\chi  e^{2s\varphi}\sigma^3)_x\right)_xs^3\lambda^3yz\\+\int_0^T\int_0^1 s^3\lambda^3(\chi  e^{2s\varphi}\sigma^3)_xaz_xy.
		\end{multline*}
		and, recalling  that $e^{2s\varphi} $ vanishes at $0$ and $T$ and  using the second equation of \eqref{adj}, we have
		\begin{align}
		I_2= & -\int_0^T\int_0^1 \left[\chi s^3\lambda^3(e^{2s\varphi}\sigma^3)_{_t}+\chi s^3\lambda^3\sigma^3e^{2s\varphi}b_{22}\right]yz +\int_0^T\int_0^1\chi s^3\lambda^3\sigma^3e^{2s\varphi}(az_x)_xy\nonumber\\
		&\ \ \ \ -\int_0^T\int_0^1\chi s^3\sigma^3e^{2s\varphi}b_{12}y^2 +\int_0^T\int_0^1\chi s^3\lambda^3\sigma^3e^{2s\varphi}yF_2\nonumber\\

		=& -\int_0^T\int_0^1 \left[\chi s^3\lambda^3(e^{2s\varphi}\sigma^3)_{_t}+\chi s^3\lambda^3\sigma^3e^{2s\varphi}b_{22}\right]yz  -\int_0^T\int_0^1 s^3\lambda^3(\chi e^{2s\varphi}\sigma^3)_xaz_xy\nonumber\\
		& -\int_0^T\int_0^1 \chi s^3\lambda^3\sigma^3e^{2s\varphi}az_xy_x -\int_0^T\int_0^1\chi s^3\lambda^3\sigma^3e^{2s\varphi}b_{12}y^2 +\int_0^T\int_0^1\chi s^3\lambda^3\sigma^3e^{2s\varphi}yF_2. \label{I2}
		\end{align}

		Thus,
		\begin{align}\label{I2I3}
		I_2+I_3= &  \int_0^T\int_0^1 \left[-\chi s^3\lambda^3(e^{2s\varphi}\sigma^3)_{_t}-\chi s^3\lambda^3\sigma^3e^{2s\varphi}b_{22}+s^3\lambda^3\left(a(\chi  e^{2s\varphi}\sigma^3)_x\right)_x\right]yz \nonumber\\
		& -2\int_0^T\int_0^1 \chi s^3\lambda^3\sigma^3e^{2s\varphi}az_xy_x -\int_0^T\int_0^1\chi s^3\lambda^3\sigma^3e^{2s\varphi}b_{12}y^2\nonumber\\
		& +\int_0^T\int_0^1\chi s^3\lambda^3\sigma^3e^{2s\varphi}yF_2 \nonumber \\
		=& J_1+J_2+J_3+J_4.
		\end{align}

		Now, it remains estimates these four integrals. It is immediate that
		\begin{equation}\label{J3}
		J_3\leq C\int_0^T\int_{\omega} e^{2s\varphi} s^8\lambda^8\sigma^8y^2.
		\end{equation}
		and, from Young's inequality, that
		\[J_4\leq \int_0^T\int_0^1 e^{2s\varphi}s^4\lambda^4\sigma^4|F_2|^2+\int_0^T\int_{\omega} e^{2s\varphi}s^8\lambda^8\sigma^8 y^2.\]

		In order to estimate $J_1$, we will analyze each term between brackets. Firstly, we observe that all the terms are  multiplied by $\chi$, which vanishes outside of $\omega$. Clearly,
		\[|\chi s^3\lambda^3\sigma^3e^{2s\varphi}b_{22}|\leq C\chi s^5\lambda^5 \sigma^5 e^{2s\varphi}. \]
		Since $|\sigma_x|\leq C\lambda \sigma$, $|\sigma_{xx}|\leq C\lambda^2 \sigma$  and $a\in C^1(\omega)$, after distributing the  derivatives with respect to $x$, we can see that
		\[|s^3\lambda^3\left(a(\chi  e^{2s\varphi}\sigma^3)_x\right)_x|\leq  C\chi s^5\lambda^5 \sigma^5 e^{2s\varphi}.\]
		Likewise, the relations $|\varphi_t|\leq C\sigma^2$ and  $|\sigma_t|\leq C\sigma^2$ yield
		\[|-\chi s^3\lambda^3(e^{2s\varphi}\sigma^3)_{_t}|\leq  C\chi s^5\lambda^5 \sigma^5 e^{2s\varphi}.\]

		As a conclusion,
		\begin{equation*}
		J_1\leq C\int_0^T\int_0^1 \chi s^5\lambda^5 \sigma^5 e^{2s\varphi}|yz|\leq \varepsilon I(s,\lambda,z)+C_\varepsilon \int_0^T\int_{\omega} s^8\lambda^8\sigma^8e^{2s\varphi}y^2.
		\end{equation*}

		The last step is to deal with $J_2$. To do this, we notice that
		\begin{align*}
		J_2\leq & \varepsilon\int_0^T\int_0^1 e^{2s\varphi}s\lambda\sigma az_x^2+C_\varepsilon\int_0^T\int_0^1 \chi^5e^{2s\varphi}s^5\lambda^5\sigma^5y_x^2\\
		\leq & \varepsilon I(s,\lambda,z)+C_\varepsilon\int_0^T\int_0^1 \chi^5e^{2s\varphi}s^5\lambda^5\sigma^5y_x^2.
		\end{align*}
		Hence, we only need to estimate the last integral. Multiplying the first equation in \eqref{adj} by $\chi^2e^{2s\varphi}s^5\lambda^5\sigma^5y$, integrate over $(0,T)\times (0,1)$ and integrating by parts we get that
		\begin{align*}
		&    \int_0^T\int_0^1 \chi^2e^{2s\varphi}s^5\lambda^5\sigma^5y_x^2 \\
		&\phantom{\leq } \leq  \int_0^T\int_0^1 \left[-\frac{1}{2}\chi^2s^5\lambda^5(e^{2s\varphi}\sigma^5)_{_t}+ \frac{1}{2}s^5\lambda^5\left[a(e^{2s\varphi}\sigma^5\chi^2)_{_x}\right]_x-\chi^2b_{11}e^{2s\varphi}s^5\lambda^5\sigma^5\right]y^2\\
		&\phantom{\leq }\phantom{\leq } -\int_0^T\int_0^1 \chi^2b_{21}e^{2s\varphi}s^5\lambda^5\sigma^5yz+\int_0^T\int_0^1 \chi^2e^{2s\varphi}s^5\lambda^5\sigma^5yF_1.
		\end{align*}

		We can see that all the integrals here are of the same type of those in \eqref{I2I3}. Following the same arguments developed there, we have the result.
	\end{proof}

	Now we need to prove a Carleman inequality  for solutions of problem \eqref{adj}  with weights which do not vanish at $t=0$. It is necessary in order to guarantee  the null controllability results in Theorem {\ref{linearcontrol}} and Theorem {\ref{main}}. We will give more details in Remark {\ref{rem1}}.

	Consider a function $m\in C^\infty([0,T])$ satisfying

	\[\left\{\begin{array}{ll}
	m(t)\geq t^4(T-t)^4, & t\in (0, T/2];\\
	m(t)= t^4(T-t)^4, & t\in \left[T/2,T\right];\\
	m(0)>  0, &
	\end{array}
	\right.\]
	and define

	\[\tau(t):=\frac{1}{m(t)},\ \ \zeta(x,t):=\tau(t)\eta(x)\ \ \mbox{ and }\ \ A(t,x):=\tau(t)(e^{\lambda(|\psi|_\infty+\psi)}-e^{3\lambda|\psi|_\infty}),\]
	where $(t,x)\in [0,T)\times [0,1]$. As usual, we introduce the operators
	\begin{align*}
	&\Gamma(s,\xi,\vartheta):=\G{\xi}{\vartheta}, \\
	&\Gamma_1(s,\xi,\vartheta):=\Ga{\xi}{\vartheta}\\
	\text{and }& \\
	&  \Gamma_2(s,\xi,\vartheta):=\Gb{\xi}{\vartheta}.
	\end{align*}

	\begin{prop}\label{prop3.2}
		There exist positive constants $C,\lambda_0$ and $s_0$ such that, for any $s\geq s_0$, $\lambda\geq \lambda_0$ and any $y_{_T},z_{_T}\in L^2((0,T)\times (0,1))$ the corresponding solution $(y,z)$ to \eqref{adj} satisfies
		\begin{multline}
		\G{y}{z}\\
		\leq C\left(\int_0^T\int_0^1 e^{2sA}s^4\lambda^4\zeta^4(|F_1|^2+|F_2|^2)\   +\int_0^T\int_{\omega} e^{2sA}s^8\lambda^8\zeta^8|y|^2\   \right).\label{carlA}
		\end{multline}
	\end{prop}

	\begin{proof} In order to estimate $\Gamma_2(s,y,z)$, let us   observe that   $e^{2s\varphi}\sigma^n\leq Ce^{2sA}\zeta^n$ for all $(t,x)\in [0,T]\times[0,1]$ and $n\geq 0$. Since $\tau=\theta$ and $A=\varphi$ in $[T/2,T]$, Carleman inequality \eqref{carl} implies
		\[\Gamma_2(s,y,z)\leq  C\left(\int_0^T\int_0^1 e^{2sA}s^4\lambda^4\zeta^4(|F_1|^2+|F_2|^2)\   +\int_0^T\int_{\omega} e^{2sA}s^8\lambda^8\zeta^8|y|^2\   \right).\]

		Now, we will prove an analogous estimate for $\Gamma_1(s,y,z)$, arguing as in \cite{clark2013theoretical}. Multiplying the first and the second equations of \eqref{adj} by $-y$ and $-z$, respectively, and integrating over $[0,1]$, we obtain
		\begin{multline}\label{eq3.5}
		-\frac{d}{dt}\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right)-C \left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right)\\+2\left(\n{\sqrt{a}y_x}{L^2(0,1)}^2+\n{\sqrt{a}z_x}{L^2(0,1)}^2\right)
		\leq \n{F_1}{L^2(0,1)}^2+\n{F_2}{L^2(0,1)}^2,
		\end{multline}
		which implies
		\[-\frac{d}{dt}\left[e^{Ct}\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right)\right]\leq e^{Ct}\left(\n{F_1}{L^2(0,1)}^2+\n{F_2}{L^2(0,1)}^2\right).\]

		Integrating from a $t\in [0,T/2]$ to $t+T/4$, we get
		\begin{multline*}
		\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2
		\leq   e^{CT}\int_0^{3T/4}\left( \n{F_1}{L^2(0,1)}^2+\n{F_2}{L^2(0,1)}^2\right)\\
		+e^{3CT/4}\left( \n{y(t+T/4)}{L^2(0,1)}^2+\n{z(t+T/4)}{L^2(0,1)}^2\right).
		\end{multline*}

		Hence, we conclude that
		\begin{multline}\label{eq3.6}
		\int_{0}^{T/2}	\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right)\leq  e^{CT}\frac{ T}{2}\int_0^{3T/4}\left( \n{F_1}{L^2(0,1)}^2+\n{F_2}{L^2(0,1)}^2\right)\\
		+e^{3CT/4}\int_{T/4}^{3T/4}\left( \n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right).
		\end{multline}

		Now, integrating inequality \eqref{eq3.5}  over $[0,t]$, where $t\in [0,T]$, we take
		\begin{multline}\label{eq3.6-2}
		\int_0^t\left(\n{\sqrt{a}y_x}{L^2(0,1)}^2+\n{\sqrt{a}z_x}{L^2(0,1)}^2\right)\leq \frac{1}{2}\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right)
		\\+C\left[ \int_0^t\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right) +\int_0^t\left(\n{F_1}{L^2(0,1)}^2+\n{F_2}{L^2(0,1)}^2\right)\right].
		\end{multline}

		In order to establish our next inequality, first we recall that if a function $f$ is non-negative, then the function $t\mapsto \int_{t_0}^t f(t) $ is non-decreasing. As a consequence, for all $  t\in [T/2,3T/4]$, we have that
		\begin{align*}
		\int_0^{T/2}\left(\n{\sqrt{a}y_x}{L^2(0,1)}^2+\n{\sqrt{a}z_x}{L^2(0,1)}^2\right) & \leq \int_{0}^{t}\left(\n{\sqrt{a}y_x}{L^2(0,1)}^2+\n{\sqrt{a}z_x}{L^2(0,1)}^2\right)\\
		& 	\leq \int_{0}^{3T/4}\left(\n{\sqrt{a}y_x}{L^2(0,1)}^2+\n{\sqrt{a}z_x}{L^2(0,1)}^2\right).	\end{align*}
		Thus, integrating inequality \eqref{eq3.6-2} from $T/2$ to $3T/4$ and  using \eqref{eq3.6} we have
		\begin{multline}\label{eq3.7}
		\int_0^{3T/4}\left(\n{\sqrt{a}y_x}{L^2(0,1)}^2+\n{\sqrt{a}z_x}{L^2(0,1)}^2\right)\\ \leq C\left[\int_{T/2}^{3T/4}\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right)+\int_0^{3T/4}\left(\n{F_1}{L^2(0,1)}^2+\n{F_2}{L^2(0,1)}^2\right)\right].
		\end{multline}

		Finally, we observe that $e^{2sA}(s\lambda \zeta)^n$ and $e^{2s\varphi}(s\lambda \sigma)^n$ are bounded in $[0,T/2]$ and $[T/4,3T/4]$ \\ respectively, for all $n\in \mathbb{Z}$. Hence, \eqref{eq3.6}, \eqref{eq3.7} and Carleman Inequality \eqref{carl} imply
		\begin{align*}
		\Gamma_1(s,y,z)& =\Ga{y}{z} \\
		& \leq C\left(\int_{T/4}^{3T/4}\int_{0}^{1}\left(|y|^2+|z|^2\right)+\int_{0}^{3T/4} \int_{0}^{1}\left(|F_1|^2+|F_2|^2\right)\right)\\
		& \leq C\left(\int_{T/4}^{3T/4}\int_{0}^{1}e^{2s\varphi}(s\lambda)^2\sigma^2\left(|y|^2+|z|^2\right) +\int_{0}^{3T/4}\int_{0}^{1}e^{2sA}(s\lambda)^4\sigma^4\left(|F_1|^2+|F_2|^2\right)\right)\\
		& 	\leq C\bigg(\int_0^T\int_0^1 e^{2sA}s^4\lambda^4\zeta^4(|F_1|^2+|F_2|^2)\   +\int_0^T\int_{\omega} e^{2sA}s^8\lambda^8\zeta^8|y|^2\bigg),
		\end{align*}
		which concludes the proof.
	\end{proof}

	\begin{cor}\label{DO-2}
		There exist positive constants $C,\lambda_0$ and $s_0$ such that, for any $s\geq s_0$, $\lambda\geq \lambda_0$ and any  $y_{_T},z_{_T}\in L^2(0,1)$ the corresponding solution $(y,z)$ to \eqref{adj},  with $F_1\equiv F_2\equiv0$,  satisfies
		\begin{equation}\label{obs.ineq}
		\n{y(0)}{L^2(0,1)}^2+\n{z(0)}{L^2(0,1)}^2\leq C\int_0^T\int_{\omega} e^{2sA}s^8\lambda^8\zeta^8|y|^2
		\end{equation}
	\end{cor}

	\begin{proof}

		Using  standard energy inequalities for each equations in \eqref{adj}, we obtain
		\[-\frac{d}{2dt}\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right)\leq 2C \left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right).\]
		Hence, we get that
		\begin{equation}\label{eq1}
		\n{y(0)}{L^2(0,1)}^2+\n{z(0)}{L^2(0,1)}^2\leq e^{4CT}\left(\n{y}{L^2(0,1)}^2+\n{z}{L^2(0,1)}^2\right).
		\end{equation}
		Finally, integrating the last inequality in $[0,3T/4]$,  recalling that $e^{2sA}(s\lambda)^2\zeta^2$ is bounded from below in $[0,3T/4]$ and using \eqref{carlA} with $F_1\equiv F_2\equiv 0$,  we obtain
		\begin{equation*}
		\n{y(0)}{L^2(0,1)}^2+\n{z(0)}{L^2(0,1)}^2 \leq C\int_{0}^{\frac{3T}{4}}\int_{0}^{1}e^{2sA}(s\lambda)^2 \zeta^2(|y|^2+|z|^2) \leq  C\int_0^T\int_{\omega} e^{2sA} s^8\lambda^8\zeta^8|y|^2.
		\end{equation*}

	\end{proof}

	\section{Global null controllability for the linear system} \label{sec-NCLS}

	\noindent

	The goal of this section is to prove a  null controllability result for the linear system \eqref{pb-lin} and establish some important additional estimates.  In order to state this result, we need to define the weights functions:
	\[\rho=e^{-sA},\ \rho_0=e^{-sA}\zeta^{-1},\ \hat{\rho}=e^{-sA}\zeta^{-5/2} \text{ and } \rho_\ast=e^{-sA}\zeta^{-4}, \]
	which satisfy $\rho_\ast\leq C\hat{\rho}\leq C\rho_0\leq C\rho$ and $\hat{\rho}^2=\rho_\ast\rho_0$.

	\begin{thm} \label{linearcontrol}
		If $u_0,v_0\in H_a^1(0,1)$ and the functions $g_1$ and $g_2$ fulfill
		\[\int_0^T\int_0^1 \rho_0^2\left(|g_1|^2+|g_2|^2\right)<+\infty,\]
		then the system \eqref{pb-lin} is null-controllable. More precisely, there exists a control $h\in L^2((0,T)\times \omega)$ with associated state $(u,v)$ satisfying
		\begin{equation}\label{eq25}
		\int_0^T\int_{\omega}\rho_\ast^2|h|^2<+\infty \text{ and } \int_0^T\int_0^1 \rho_0^2(|u|^2+|v|^2)<+\infty.
		\end{equation}
		In particular, $u(T,x)=v(T,x)=0,$ for all $x\in [0,1]$.
	\end{thm}

	\begin{rem}\label{rem1}
		\noindent

		\begin{itemize}
			\item[(a)]
				Recalling that  $\rho_0(t)\to +\infty$, as $t\to T^{-}$, and $\rho_{0} (0) >0$ (since $m(0)>0$), the second relation in ($\ref{eq25}$) garantees $u(T,x)=v(T,x)=0$.
			\item[(b)]  If we had chosen $m\in C^{\infty} ([0,T])$ satisfying $m(0)=0$, we would verify $\rho_0(t)\to +\infty$, as $t\to 0^{+}$. As a consequence, the second relation in ($\ref{eq25}$) would imply $u(0,x)=v(0,x)\equiv 0$. However, in general, this fact is not true, because $u_{0} \in H_{a}^{1} (0,1)$ and $v_{0} \in H_{a}^{1} (0,1)$ must be taken arbitrarily.
		\end{itemize}

	\end{rem}

	\begin{proof}
	For  each $n\in \mathbb{N}^\ast$, we define
		\[A_n(t,x)=\frac{A(T-t)^4}{(T-t)^4+\frac{1}{n}},\ (t,x) \in  [0,T]\times [0,1].

		\]
		We also consider
		\[\rho_n=e^{-sA_n}, \ \rho_{0,n}=\rho_{n}\zeta^{-2} \text{ and } \rho_{\ast,n}=\rho_n\zeta^{-4} m_n, \text{ where }m_n=\begin{cases}
		1,&  x\in \omega,\\
		n,&  x\notin \omega.
		\end{cases}\]

		These weight functions are built in such a way that $\rho_{0,n}$ and $\rho_{\ast, n}$ are  bounded from below by a constant depending only on $T$ and from above by another one depending on $n$ and $T$, see Lemma \ref{lemma-const}. It will allows us to obtain a sequence $(u_n,v_n,h_n)_{n\in \mathbb{N}^\ast}$ of solutions to \eqref{pb-lin} which will converge to a solution of \eqref{pb-lin} satisfying \eqref{eq25}.

		To do that,  for any functions $u,v,h\in L^2((0,T)\times (0,1))$, let us  define the functional
		\[J_n(u,v,h)=\frac{1}{2}\int_0^T\int_0^1 \rho_{0,n}^2\left(|u|^2+|v|^2\right)+\frac{1}{2}\int_0^T\int_0^1 \rho_{\ast,n}^2|h|^2.\]

		Since each $J_n$ is lower semi-continuous, strictly convex and coercive (see Appendix \ref{Jn}), Proposition 1.2 in \cite{ekeland1999convex} yields  a unique $(u_n,v_n,h_n)$ that minimizes $J_n(u,v,h)$ subject to the condition $\mathcal{C}=\{(u,v,h)\in \mathbf{ \left[L^2((0,T)\times (0,1))\right]^3}; (u,v,h)\text{ solves } \eqref{pb-lin}\}$. In this case, $(u_n,v_n,h_n)$ satisfies \eqref{pb-lin} and, by virtue of \textit{Lagrange's Principle}, there exist  functions $p_n,q_n$ solving the following system
		\begin{equation}\label{S2}
		\begin{cases}
		-{p_{n}}_t-({ap_{n}}_x)_x+b_{11}p_{n}+b_{21}q_{n}=-\rho_{0,n}^2u_n, & (t,x) \in (0,T)\times (0,1),\\
		-{q_{n}}_t-(a{q_n}_x)_x+b_{12}p_{n}+b_{22}q_{{n}}=-\rho_{0,n}^2v_n, & (t,x)\in (0,T)\times (0,1),\\
		p_n(t,0)=p_n(t,1)=q_n(t,0)=q_n(t,1)=0,& t\in (0,T),\\
		p_n(T,x)=q_n(T,x)=0, & x\in (0,1),\\
		p_n=\rho_{_{\ast,n}}^2h_n \chi_{\omega }, & \in (0,T)\times (0,1).
		\end{cases}
		\end{equation}

		We wan to prove that $(J_n(u_n,v_n,h_n))_{n=1}^\infty$ is a numerical bounded sequence. To do that, we will use Carleman  and   Observability inequalities to prove that  $J_{n} (u_{n},v_n ,h_{n} ) \leq C \sqrt{J_{n} (u_{n},v_n ,h_{n} )}$.

		In fact, multiplying the PDEs in (\ref{S2}) by $u_{n}$ and $v_n$, integrating over $(0,T)\times (0,1)$ and  using integration by parts, we have
		\begin{align*}
		0&=
		\begin{multlined}[t][0.9\textwidth]
		\int_0^T\int_0^1 [-{p_{n}}_t - (a {p_{n}}_x )_{x} +b_{11} p_{n} +b_{21}q_n +\rho_{0,n}^{2} u_{n}] u_{n} \\
		+	\int_0^T\int_0^1 [-{q_{n}}_t - (a {q_{n}}_x )_{x} +b_{12} p_{n} +b_{22}q_n  +\rho_{0,n}^{2} v_{n}] v_{n}
		\end{multlined}\\
		&=\begin{multlined}[t][0.9\textwidth]
		\int_0^T\int_0^1 [{u_{n}}_t - (a {u_{n}}_x )_{x} +b_{11}u_{n}+b_{12}v_n ]p_{n} + \int_0^T\int_0^1 \rho_{0,n}^{2} |u_{n} |^{2} + \int_{0}^{1} p_{n} (0,x) u_{0} \\
		+\int_0^T\int_0^1 [{v_{n}}_t - (a {v_{n}}_x )_{x} +b_{21}u_{n}+b_{22}v_n ]q_{n} + \int_0^T\int_0^1 \rho_{0,n}^{2} |v_{n} |^{2} + \int_{0}^{1} q_{n} (0,x) v_{0}
		\end{multlined}
		\end{align*}

		\begin{align*}
		&=\begin{multlined}[t][0.9\textwidth]
		\int_0^T\int_0^1 (h_{n} \chi _{\omega} +g_1) p_{n} + \int_0^T\int_0^1 \rho_{0,n}^{2} |u_{n} |^{2} + \int_{0}^{1} p_{n} (0,x) u_{0}\\
		+\int_0^T\int_0^1 g_2 q_{n} + \int_0^T\int_0^1 \rho_{0,n}^{2} |v_{n} |^{2} + \int_{0}^{1} q_{n} (0,x) v_{0}
		\end{multlined}
		\end{align*}
		Hence, since $p_n=p_{\ast,n}^2h_n\chi_{\omega }$, we obtain
		\begin{align*}
		\displaystyle J_{n} (u_{n}, v_n ,h_{n}) & =\frac{1}{2}\int_0^T\int_0^1 \rho_{0,n}^2\left(|u_n|^2+|v_n|^2\right)+\frac{1}{2}\int_0^T\int_0^1 \rho_{\ast,n}^2|h_n|^2\chi_{\omega } \\
		&\begin{multlined}[t][0.8\textwidth]
		=-\frac{1}{2} \int_0^T\int_0^1 (h_{n} \chi _{\omega} +g_1) p_{n} -\frac{1}{2} \int_{0}^{1} p_{n} (0,x) u_{0} \\
		-\frac{1}{2} \int_0^T\int_0^1 g_2q_{n} - \frac{1}{2} \int_{0}^{1} q_{n} (0,x)v_{0} + \frac{1}{2} \int_0^T\int_0^1 p_{n} h_{n}\chi_\omega
		\end{multlined}\\
		& =-\frac{1}{2} \int_0^T\int_0^1 (g_1 p_{n} + g_2q_{n}) -\frac{1}{2} \int_{0}^{1} p_{n} (0,x) u_{0} - \frac{1}{2} \int_{0}^{1} q_{n} (0,x)v_{0}.
		\end{align*}
		Using  H\"older inequality, we have
		\begin{multline*}
		\displaystyle J_{n} (u_{n},v_n ,h_{n})		 \leq \frac{1}{2} \n{g_1\rho _{0}}{ L^2}\n{\rho _{0}^{-1} p_{n}}{L^2}+\frac{1}{2} \n{g_2\rho _{0}}{ L^2}\n{\rho _{0}^{-1} q_{n}}{L^2}\\
		+\frac{1}{2}\n{p_{n} (0,\cdot)}{L^2(0,1)}\n{u_{0}}{L^2(0,1)} +\frac{1}{2}\n{q_{n} (0,\cdot)}{L^2(0,1)}\n{v_{0}}{L^2(0,1)}.
		\end{multline*}

		Applying the classical Cauchy-Schwartz inequality $\Sigma_{i=1}^4(a_ib_i)\leq (\Sigma_{i=1}^4a_i^2)^{1/2}(\Sigma_{i=1}^4b_i^2)^{1/2}$ and then using hypothesis $\displaystyle \int_0^T\int_0^1 \rho _{0}^{2}( |g_1|^{2}+|g_2|^{2}) <+\infty$ , we get
		\begin{align*}
		\displaystyle J_{n} (u_{n},v_n ,h_{n})	 &\begin{multlined}[t][0.9\textwidth]
		\leq \frac{1}{2} \left(\n{g_1\rho _{0}}{ L^2}^2+\n{g_2\rho _{0}}{ L^2}^2 +\|u_{0}\|^2 +\|v_{0}\|^2\right)^{1/2}\\
		\cdot\left(\n{\rho _{0}^{-1} p_{n}}{L^2}^2 +\n{\rho _{0}^{-1} q_{n}}{L^2}^2 + \n{p_{n} (0,\cdot)}{L^2(0,1)}^2+ \n{q_{n} (0,\cdot)}{L^2(0,1)}^2 \right)^{1/2}
		\end{multlined} \\
		&\leq C \left( \int_0^T\int_0^1 \rho _{0}^{-2} (p_{n}^{2}+ q_{n}^{2})  +\n{p_{n} (0,\cdot)}{L^2(0,1)}^2+ \n{q_{n} (0,\cdot)}{L^2(0,1)}^2 \right) ^{1/2}.
		\end{align*}

		Now, it is enough to prove that each term in the last inequality is  bounded by $J(u_n,v_n,h_n)$. In order to estimate the first term, we will apply Carleman inequality \eqref{carlA}  to the solution $(p_n,q_n)$ to \eqref{S2} and then use that  $\rho^{-2}\rho_n^4 \leq \rho_n^{2}$  and  $p_n=\rho_{\ast,n}^2h_n\chi_{\omega }$ in $(0,T)\times (0,1)$.

		Indeed,
		\begin{align*}
		\displaystyle \int_0^T\int_0^1 \rho _{0}^{-2} (p_{n}^{2}+q_n^2)
		&=\frac{1}{(s\lambda)^2} \int_0^T\int_0^1 e^{2s A}(s\lambda)^2\zeta^2(p_{n}^{2}+q_n^2)\\
		&\leq C\left( \int_0^T\int_0^1 e^{2s A}\zeta^4(s\lambda)^4\rho_{0,n}^{4}(  |u_{n} |^{2}+|v_n|^2) + \int_0^T\int_{\omega}  e^{2s  A}\zeta^8(s\lambda)^8 |p_{n} |^{2} \right) \\
		&\leq C \left( \int_0^T\int_0^1 \rho^{-2} \rho_{n}^{4}\zeta^{-4}( |u_{n} |^{2}+|v_{n} |^{2}) + \int_0^T\int_{\omega} \rho^{-2}\zeta^8 \rho _{\ast ,n}^{4} h_{n}^{2} \right) \\
		&\leq C \left( \int_0^T\int_0^1 \rho_{0,n}^{2} ( |u_{n} |^{2}+|v_{n} |^{2}) + \int_0^T\int_{\omega} \rho _{\ast ,n}^{2} |h_{n} |^{2} \right) \\
		&=C J_{n} (u_{n}, v_n ,h_{n}).
		\end{align*}

		The remain terms are readily estimated  from the Observability Inequality \eqref{obs.ineq}. In fact,
		\begin{align*}
		\n{p_{n} (0,\cdot)}{L^2(0,1)}^2+ \n{q_{n} (0,\cdot)}{L^2(0,1)}^2 & \leq C\left( \int_0^T\int_{\omega}  e^{2s  A}\zeta^8(s\lambda)^8 |p_{n} |^{2} \right)\leq C J_{n} (u_{n}, v_n ,h_{n}).
		\end{align*}

		Hence, we have proven that $J_{n} (u_{n},v_n ,h_{n} ) \leq C \sqrt{J_{n} (u_{n},v_n ,h_{n} )}$.  As a consequence, $(J_n(u_n,v_n,h_n))_{n\in \mathbb{N}}$ is a bounded sequence. Since $\rho_{0,n}^2\geq C_T$ and $\rho_{\ast,n}^2\geq C_Tm_n$, we deduce that
		\[\n{u_n}{L^2}^2+\n{v_n}{L^2}^2+\int_0^T\int_{\omega} |h_n|^2 +n\int_0^T\int_{[0,1]\setminus\omega}|h_n|^2\leq C J_{n} (u_{n},v_n ,h_{n} )\leq  C.
		\]

		It implies that   there exist $u,v\in L^2((0,T)\times (0,1))$ and $h\in L^2((0,T)\times \omega)$  such that
		\[u_n\rightharpoonup u, \ \ v_n\rightharpoonup v \ \mbox{ in } L^2((0,T)\times (0,1)) \ \mbox{ and } \ \ h_n\chi_{\omega }\rightharpoonup h\ \ \mbox{ in } L^2((0,T)\times \omega),\]
		up to subsequences.
		From this, we take
		\begin{equation}\label{wc}
		\rho_{0,n}u_n\rightharpoonup \rho_0 u \ \ \rho_{0,n}v_n\rightharpoonup \rho_0 v \ \ \mbox{ and } \ \ \rho_{\ast,n}h_n\chi_{\omega }\rightharpoonup \rho_{\ast}h\chi_{\omega } \ \ \mbox{ in } L^2((0,T)\times (0,1)).
		\end{equation}

		Consequently, passing to limits as $n\to +\infty$, we conclude that $(u,v,h)$ solves \eqref{pb-lin}.		Furthermore, (\ref{eq25}) follows from  (\ref{wc}) and this establishes the result.
			\end{proof}

	%
	The next step is to prove two crucial estimates which will needed later.

	\begin{prop}\label{ad1}
		Assume the same hypothesis of Theorem $\ref{linearcontrol}$. Then

		\begin{multline}
		\int_0^T\int_0^1 \hat{\rho} ^2 a(|u_x |^2 + |v_x |^2 )
		\leq C \left( \int_0^T\int_0^1 \rho_{0}^{2} (|u|^2 + |v|^2 ) + \int_0^T\int_{\omega}  \rho_{\ast}^{2} |h|^2 \right) \\
		+C\left( \int_0^T\int_0^1 \rho _{0}^{2} (|g_1 |^2 + |g_2 |^2 ) + \| u_0 \| _{H_{a}^{1}}^{2} + \| v_0 \| _{H_{a}^{1}}^{2} \right)  \ .
		\end{multline}
	\end{prop}

	\begin{proof}
		Let us multiply the first equation in (\ref{pb-lin}) by $\hat{\rho} ^2 u$ and the second one by $\hat{\rho} ^2 v$, and let us integrate over $[0,1]$. In this case, we obtain
		\begin{multline}\label{eq2}
		\int_{0}^{1}\hat{\rho} ^2 [u_t u + v_t v ]-\int _{0}^{1} \hat{\rho} ^2 [(au_x) _x u + (av_x )_x v]\\
		= - \int_{0}^{1} \hat{\rho}^2[(b_{11} u + b_{12} v)u+ (b_{21} u + b_{22} v)v]  + \int _{0}^{1} \hat{\rho} ^2 h \chi _{\omega } u+ \int _{0}^{1} \hat{\rho} ^2 [g_1 u + g_2 v].
		\end{multline}
		Clearly, the terms in the right hand side of \eqref{eq2} can be estimated as follows:
		\begin{align*}
		& \displaystyle \left|  \int_{0}^{1} \hat{\rho}^2[(b_{11} u + b_{12} v)u+ (b_{21} u + b_{22} v)v]  \right| \leq C \int_{0}^{1} \hat{\rho} ^2 (|u|^2 + |v|^2)  ,\\
		&\displaystyle \int _{0}^{1} \hat{\rho} ^2 h \chi _{\omega } u   = \int _{\omega} (\rho _{0} u)(\rho _{\ast } h)
		\leq \frac{1}{2} \int _{\omega} \rho _{0}^{2} |u|^2
		+\frac{1}{2} \int _{\omega} \rho _{\ast }^{2} |h|^2
		\end{align*}
		and
		\begin{multline*}
		\displaystyle \int _{0}^{1} \hat{\rho} ^2 [g_1 u + g_2 v]
		\leq \frac{1}{2} \int _{0}^{1} \hat{\rho} ^2 (|g_1 |^2 + |u|^2 +|g_2 |^2 + |v|^2 ) \\ \leq C\left( \int_{0}^{1} \rho_{0}^{2} (|g_1 |^2 + |g_2 |^2)  + \int_{0}^{1} \rho_{0}^{2} (|u|^2 + |v|^2) \right) .
		\end{multline*}

		Now, let us deal with the left hand side of \eqref{eq2}. Notice that
		\begin{align*}
		\displaystyle \int _{0}^{1} \hat{\rho} ^2 [u_t u + v_t v ]
		&=\frac{1}{2} \frac{d}{dt} \int_{0}^{1} \hat{\rho} ^2 (|u|^2 + |v|^2)
		-\int_{0}^{1} \hat{\rho} (\hat{\rho} )_t  (|u|^2 + |v|^2 )  \\
		&:= \frac{1}{2} \frac{d}{dt} \int_{0}^{1} \hat{\rho} ^2 (|u|^2 + |v|^2)
		-\mathcal{I}
		\end{align*}
		and
		\begin{multline*}
		-\int _{0}^{1} \hat{\rho} ^2 [(au_x) _x u + (av_x )_x v]
		=\int_{0}^{1} [(\hat{\rho} ^2 u)_x (au_x ) + (\hat{\rho} v )_x (av_x )]  \\
		= 2\int_{0}^{1} \hat{\rho} \hat{\rho} _x a (uu_x + vv_x )  + \int_{0}^{1} \hat{\rho} ^2 a (u_{x}^{2} + v_{x}^{2 })   := \mathcal J + \int_{0}^{1} \hat{\rho} ^2 a (u_{x}^{2} + v_{x}^{2 }).
		\end{multline*}

		Summing up, we have just checked that
		\begin{multline}\label{eq3}
		\frac{1}{2} \frac{d}{dt} \int_{0}^{1} \hat{\rho} ^2 (|u|^2 + |v|^2)
		+ \int_{0}^{1} \hat{\rho} ^2 a (|u_{x}|^{2} + |v_{x}|^{2 })\\ \leq C\left(   \int_{0}^{1} \rho_{0}^{2} (|u|^2 + |v|^2)+\int _{\omega} \rho _{\ast }^{2} |h|^2+\int_{0}^{1} \rho_{0}^{2} (|g_1 |^2 + |g_2 |^2)  \right) +|\mathcal{I}|+|\mathcal J|
		\end{multline}

		Next, we will estimate $ \mathcal I$. Firstly, we put
		\[A(t,x)=\tau(t)(e^{\lambda(|\psi|_\infty+\psi)}-e^{3\lambda|\psi|_\infty}):= \zeta (t,x) \bar{\eta } (x),\]
		where $\displaystyle \bar{\eta} (x):= \frac{e^{\lambda (|\psi |_{\infty} + \psi )} - e^{3\lambda |\psi |_{\infty }}}{\eta (x)} $ is a bounded function on $[0,1]$. Secondly, we observe that
		\begin{align*}
		\displaystyle \hat{\rho} (\hat{\rho} )_t
		&=e^{-sA} \zeta ^{-5/2} \left[ -se^{-sA} A_t \zeta ^{-5/2} + e^{-sA} \left( \frac{-5}{2} \right)  \zeta ^{-7/2} \zeta _t \right]  \\
		&=-se^{-2sA} \bar{\eta } (x) \zeta ^{-5} \zeta _t - \frac{5}{2} e^{-2sA} \zeta ^{-6} \zeta _{t} \\
		&=e^{-2sA} \zeta ^{-2} \zeta _t \left[ -s\zeta ^{-3} \bar{\eta} (x) - \frac{5}{2} \zeta ^{-4} \right] \\
		&=\rho_{0}^{2} \zeta _t \left[ -s\zeta ^{-3} \bar{\eta} (x) - \frac{5}{2} \zeta ^{-4} \right]
		\end{align*}
		and, for any $t\in [0,T)$, we have
		\begin{align*}
		\displaystyle |\hat{\rho} (\hat{\rho} )_t |
		&\leq C \rho_{0}^{2} \tau ^2 \left|s\zeta ^{-3} \bar{\eta } (x) - \frac{5}{2} \zeta ^{-4} \right|
		= C \rho_{0}^{2} \zeta ^2 \frac{\left|s\zeta ^{-3} \bar{\eta } (x) - \frac{5}{2} \zeta ^{-4} \right| }{\eta ^2 } \\
		&= C \rho_{0}^{2} \frac{\left|s\zeta ^{-1} \bar{\eta } (x) - \frac{5}{2} \zeta ^{-2} \right| }{\eta ^2 }
		\leq C\rho_{0}^{2} .
		\end{align*}
		From this, we obtain
		\[|\mathcal I| \leq C\int_{0}^{1} \rho_{0}^{2} (|u|^2 + |v|^2 ).\]

		Now, in order to deal with $\displaystyle \mathcal J $, we consider the estimate
		\begin{equation*}
		\displaystyle |\mathcal J|
		\leq \frac{1}{2} \int_{0}^{1} [\hat{\rho} ^2 a (|u_x |^2 + |v_x |^2 ) + \hat{\rho} _{x}^{2} a(|u|^2 +|v|^2 )] := \frac{1}{2} \int_{0}^{1} \hat{\rho} ^2 a (|u_x |^2 + |v_x |^2 )   + \tilde{\mathcal J}
		\end{equation*}
		and recall  that $A_x = \zeta _x = \lambda \zeta \psi _x $ in $[0,T]\times [0,1]$. Hence,
		\begin{align*}
		\displaystyle \hat{\rho}_{x}^{2}
		&=\left( -e^{-sA} sA_x  \zeta ^{-5/2} -\frac{5}{2} e^{-sA} \zeta ^{-7/2} \zeta_{x} \right) ^2
		\leq C \left( e^{-2sA} s^2 A_{x}^{2}  \zeta ^{-5} +\frac{25}{4} e^{-2sA} \zeta ^{-7} \zeta_{x}^{2} \right) \\
		&=Ce^{-2sA} \zeta ^{-2} \psi _{x}^{2} \left( s^2 \lambda ^2 \zeta ^{-1} + \frac{25}{4} \lambda ^2 \zeta ^{-3} \right)
		=C\rho _{0}^{2} \psi _{x}^{2} \left( s^2 \lambda ^2 \zeta ^{-1} + \frac{25}{4} \lambda ^2 \zeta ^{-3} \right) \\
		&\leq C\rho_{0}^{2}
		\end{align*}
		and we get $\displaystyle \tilde{\mathcal J} \leq C\int_{0}^{1} \rho_{0}^{2} a (|u|^2 + |v|^2 ) $.

		Recalling inequality \eqref{eq3}, we  conclude that
		\begin{multline*}
		\displaystyle  \frac{1}{2} \frac{d}{dt} \int_{0}^{1} \hat{\rho} ^{2} (|u|^2 + |v|^2 )   + \frac{1}{2} \int_{0}^{1} \hat{\rho} ^2 a (|u_x |^2 + |v_x |^2 )   \\
		\leq C\left( \int_{0}^{1} \rho_{0}^{2} (|u|^2 + |v|^2 )    + \int_{\omega} \rho_{\ast }^2 |h|^2    + \int_{0}^{1} \rho_{0}^{2} (|g_1 |^2 + |g_2 |^2 )  \right)
		\end{multline*}
		and, integrating in time, we obtain the desired result.
	\end{proof}

	\begin{prop} \label{ad2}
		Assume the hyphotesis of Theorem $\ref{linearcontrol}$ and suppose that $h$ and $(u,v)$ satisfy ($\ref{eq25}$). Then
		\begin{multline*}
		\displaystyle \int_0^T\int_0^1 \rho_{\ast}^{2} (|u_t |^2 + |v_t |^2 + |(au_x )_x |^2 +|(av_x )_x |^2 )
		\leq C\int_0^T\int_0^1 \rho_{0}^{2} (|u|^2 + |v|^2 ) dt \\
		+C \left( \int_0^T\int_{\omega} \rho_{\ast }^{2} |h|^2 + \int_0^T\int_0^1 \rho_{0}^{2} (|g_1 |^2 + |g_2 |^2 ) + \| u_0 \|_{H_{a}^{1}}^{2} + \| v_0 \|_{H_{a}^{1}}^{2} \right) .
		\end{multline*}
	\end{prop}

	\begin{proof}
		Multiplying the first equation in (\ref{pb-lin}) by $\rho_{\ast}^{2} u_t$ and the second one by $\rho_{\ast}^{2} v_t$, ,   we take
		\begin{multline*}
		\int_{0}^{1} \rho_{\ast}^{2} (|u_t|^2 + |v_t|^2 )
		-\int_{0}^{1} \rho_{\ast }^{2} [(au_x )_x u_t +(av_x )_x v_t ]  \\
		=-\int_{0}^{1} \rho_{\ast }^{2} [(b_{11} u + b_{12} v)u_t + (b_{21} u + b_{22} v )v_t ]+\int_{0}^{1} \rho_{\ast}^{2} h \chi_{\omega } u_t  +\int_{0}^{1} \rho_{\ast}^{2} (g_1 u_t + g_2 v_t ) .
		\end{multline*}
		Notice that,
		\begin{align}\label{twice}
		\displaystyle -\int_{0}^{1} \rho_{\ast }^{2} [(au_x )_x u_t +(av_x )_x v_t ]
		=&\frac{1}{2} \frac{d}{dt} \int_{0}^{1} \rho_{\ast}^{2} a (|u_x|^2 + |v_x|^2)\\ & -\frac{1}{2} \int_{0}^{1} (\rho_{\ast}^{2} )_t a (|u_x|^2 + |v_x|^2)
		+\int_{0}^{1} (\rho_{\ast}^{2} )_x a(u_t u_x + v_t v_x ) \nonumber \\
		=:& \frac{1}{2} \frac{d}{dt} \int_{0}^{1} \rho_{\ast}^{2} a (|u_x|^2 + |v_x|^2) -\mathcal K ,
		\end{align}
		Proceeding as in the proof of Proposition \ref{ad1},
		\begin{multline}\label{eq4}
		\int_{0}^{1} \rho_{\ast}^{2} (|u_t|^2 + |v_t|^2 )
		+\frac{1}{2} \frac{d}{dt} \int_{0}^{1} \rho_{\ast}^{2} a (|u_x|^2 + |v_x|^2 )  \\
		\leq C\left(\int_{0}^{1} \rho_{\ast}^{2} (|u|^2 + |v|^2 )
		+\int_{\omega} \rho_{\ast}^{2} |h|^2
		+\int_{0}^{1} \rho_{\ast}^{2} (|g_1|^2 + |g_2|^2 ) \right)\\+\frac{3}{8} \int_{0}^{1} \rho_{\ast}^{2} (|u_t|^2 + |v_t|^2 )+|\mathcal{K}|.
		\end{multline}

		Using Young's inequality with $\varepsilon$, we have
		\begin{align}\label{ca}
		\displaystyle \mathcal K
		&=\frac{1}{2} \int_{0}^{1} (\rho_{\ast}^{2} )_t a (|u_x|^2 + |v_x|^2)
		-\int_{0}^{1} [((\rho_{\ast}^{2} )_x \rho_{\ast}^{-1} a u_x ) (\rho_{\ast} u_t ) + ((\rho_{\ast}^{2} )_x \rho_{\ast}^{-1} a v_x ) (\rho_{\ast} v_t )] \nonumber \\
		&\leq C \int_{0}^{1} [ |(\rho_{\ast}^{2} )_t | + |(\rho_{\ast}^{2})_x |^{2} \rho_{\ast}^{-2} ] a (|u_x|^2 + |v_x|^2)
		+\frac{1}{8} \int_{0}^{1} \rho_{\ast}^{2} (|u_t|^2 + |v_t|^2 ).
		\end{align}
		Since $|\zeta_t |\leq C \zeta ^2$ and $\zeta_{x} = A_x = \lambda \zeta \psi _{x}$ in $[0,T)\times [0,1]$, we have
		\[\displaystyle |(\rho_{\ast}^{2} )_t|
		=\hat{\rho} ^2 |2sA_t \zeta ^{-3} + 8\zeta ^{-4} \zeta_t |
		\leq C\hat{\rho}^2 (2s\zeta^{-1} + 8\zeta^{-2} )\]
		and
		\begin{align*}
		\displaystyle |(\rho_{\ast}^{2} )_x |^2 \rho_{\ast }^{-2}
		&\leq C(4s^2 e^{-4sA} A_{x}^{2} \zeta ^{-16} + 64 e^{-4sA} \zeta ^{-18} \zeta_{x}^{2})\rho_{\ast }^{-2} \\
		&=C\hat{\rho}^2 (4s^2 A_{x}^{2} \zeta ^{-3} + 64 \zeta ^{-5} \zeta_{x}^{2}) \\
		&=C\hat{\rho}^2 (4s^2 \lambda ^2 \psi_{x}^{2} \zeta ^{-1} + 64 \zeta ^{-3}\lambda^2 \psi_{x}^{2} ),
		\end{align*}
		which imply
		\[\displaystyle \mathcal K \leq  C \int_{0}^{1} \hat{\rho}^2 a (|u_x|^2 + |v_x|^2)
		+\frac{1}{8} \int_{0}^{1} \rho_{\ast}^{2} (|u_t|^2 + |v_t|^2 ).\]

		%

		Thus, from $\eqref{eq4}$,  applying Proposition $\ref{ad1}$ and using $\rho_{\ast } \leq C \rho_{0}$, we get
		\begin{multline}\label{part1}
		\int_0^T\int_0^1 \rho_{\ast}^{2} (|u_t |^2 + |v_t |^2 )
		\leq C\left(\int_0^T\int_0^1 \rho_{0}^{2} (|u|^2 + |v|^2 )  + \int_0^T\int_{\omega} \rho_{\ast }^{2} |h|^2\right.\\ + \left.\int_0^T\int_0^1 \rho_{0}^{2} (|g_1 |^2 + |g_2 |^2 ) + \| u_0 \|_{H_{a}^{1}}^{2} + \| v_0 \|_{H_{a}^{1}}^{2} \right) .
		\end{multline}

		In order to conclude the proof, it remains to estimate  $\int_0^T\int_0^1 \rho_{\ast}^{2} (|(au_x )_x |^2 +|(av_x )_x |^2 )$. In fact, it is enough to multiply  the first equation  in ($\ref{pb-lin}$) by $-\rho_{\ast }^{2} (au_x)_x$ and the second one by  $-\rho_{\ast }^{2} (av_x)_x$, and proceed as in the first part of this proof.
	\end{proof}

	\section{Main Result } \label{sec-main}

	In this section, our goal is to prove Theorem \ref{main}.	Let us define the functions spaces
	\begin{multline*}
	E:= \bigg\{  (u,v,h)\in\mathbf{ \left[L^2((0,T)\times (0,1))\right]^2}\times L^2((0,T)\times \omega):\\
	u(t,\cdot),v(t,\cdot)  \text{ are absolutely continuous in } [0,1], \text{ a.e. in } [0,T], \\
	u_t, u_x, (au_x)_x, \rho_{\ast} h\in L^2 ((0,T)\times (0,1)),\ \
	v_t, v_x, (av_x)_x \in L^2 ((0,T)\times (0,1)) ,\\
	\rho_{0} u, \rho_{0} [u_t - (au_x)_x -h\chi_{\omega}],
	\rho_{0} v, \rho_{0} [v_t - (av_x)_x ] \in L^2 ((0,T)\times (0,1)),\\
	u(t,1)\equiv v(t,1)\equiv u(t,0)\equiv v(t,0)\equiv 0 \text{ a.e in } [0,T], \emph{ and }
	u(0,\cdot), v(0,\cdot)\in H_{a}^{1}\bigg\},
	\end{multline*}
	and $F:=G \times G \times H_{a}^{1} \times H_{a}^{1}$, where
	\[\displaystyle G := \left\lbrace  g\in L^2((0,T)\times (0,1)): \rho_{0} g \in L^2 ((0,T)\times (0,1))   \right\rbrace .\]

	We also consider the Hilbertian norm
	\begin{multline*}
	\n{(u,v,h)}{_E}^2:=
	\int_0^T\int_0^1 \rho_{0}^{2} (|u|^2 +|v|^2)
	+\int_0^T\int_{\omega} \rho_{\ast}^{2} |h|^2 \\
	\hspace{0.5cm} +\int_0^T\int_0^1 \rho_{0}^{2} |u_t - (au_x)_x -h\chi_{\omega}|^2
	+\int_0^T\int_0^1 \rho_{0}^{2} |v_t - (av_x)_x |^2 \\
	\hspace{0.5cm}+\|u(0,\cdot)\|_{H_{a}^{1}}^{2}
	+\|v(0,\cdot)\|_{H_{a}^{1}}^{2} .
	\end{multline*}
	The proof that $E$ is a Hilbert space is given in Appendix \ref{hilbert}.

	Now, we set the mapping $H: E\to F$, given by
	\[H (u,v,h) = (H_{1} (u,v,h), H_{2} (u,v,h), u(\cdot ,0), v(\cdot ,0) ),\]
	where
	\[\displaystyle H_{1} (u,v,h) = u_t-\left(\mu_1 \left(x,\int_0^1u \   \right)u_x \right)_x+f_1(t,x,u,v) -h\chi_\omega\]
	and
	\[\displaystyle H_{2} (u,v,h) = v_t-\left(\mu_2 \left(x,\int_0^1v \   \right)v_x \right)_x+f_2(t,x,u,v) .
	\]
	Applying \textit{Lyusternik's Inverse Mapping Theorem}, see \cite{alekseev1987optimal}, we will prove that $H$ has a right inverse mapping defined in a small ball contained in  $F$.
	Due to the choice of the spaces $E$ and $F$, the existence of that inverse mapping will imply the local null controllability of ($\ref{pb1}$). Before doing it, we will establish some results which will guarantee that $H$ satisfies the hypotheses of Lyusternik's Theorem.

	\begin{lemma}\label{supremo}
		Define $\displaystyle \beta (x) =e^{\lambda (|\psi|_{\infty} +\psi)}-e^{3\lambda |\psi|_{\infty}} $ and $\bar{\beta} =\displaystyle \max_{x\in[0,1]} \beta(x)$. There exists $s>0$ such that, if $s\bar{\beta} <M<0$, then
		\[
		\displaystyle \sup_{t\in[0,T]} \left\lbrace e^{\frac{-2M}{m(t)}} \left[ \left( \int_{0}^{1} u \right)^2 +\left( \int_{0}^{1} v \right)^2 \right] \right\rbrace
		\leq C\|(u,v,h)\|_{_E}^{2},\]
		for all $(u,v,h)\in E$.
	\end{lemma}

	\begin{proof}
		In fact, for each $(u,v,h)\in E$, consider $q_1 : [0,T]\longrightarrow \R$ and $q_2 : [0,T]\longrightarrow \R$, given by
		\[
		\displaystyle q_1 (t) = e^{\frac{-M}{m(t)}} \int_{0}^{1} u(t,x)
		\emph{  and  }
		q_2 (t) = e^{\frac{-M}{m(t)}} \int_{0}^{1} v(t,x) .
		\]
		Taking $k>0$, we quickly get $\displaystyle e^{\frac{-k}{m(t)}} \leq 8! [m(t)]^8/k^8$,
		for any $t\in [0,T]$. Since $A=\tau \beta$, taking $s>0$ such that $2s (\bar{\beta} - \beta) >k$, we have
		\[\displaystyle -\frac{2M}{m(t)} +2sA =-\frac{2M}{m(t)} +\frac{2s\beta}{m(t)} \leq \frac{-2s(\bar{\beta} -\beta )}{m(t)} <\frac{-k}{m(t)}\]
		and
		\[
		\displaystyle e^{\frac{-2M}{m(t)}} < e^{-2sA} e^{\frac{-k}{m(t)}} \leq Ce^{-2sA} \tau^{-8} \leq C\rho_{\ast}^{2} ,
		\]
		for any $t\in [0,T]$. From this point, we may argue as in \cite{jrl2016} (see Lemma 4.4, on page 533), in order to check that $q_1,q_2\in H^1(0,T)$ and
		\[
		\|q_1\|_{H^1(0,T)} +\|q_2\|_{H^1(0,T)}
		\leq C\|(u,v,h)\|_{E}^{2}.
		\]
		Therefore, the desired result is a consequence of the continuous embedding $ H^1(0,T)\hookrightarrow C(0,T)$.
	\end{proof}

	As a consequence of Lemma $\ref{supremo}$, we deduce the useful result below:

	\begin{cor}\label{conseq}
		There exists $C>0$ such that

		\begin{equation*}
		\int_0^T\int_0^1 \rho_{0}^{2} \left( \int_{0}^{1} \bar{u} \right) ^2 |(au_x)_x|^2
		+\int_0^T\int_0^1 \rho_{0}^{2} \left( \int_{0}^{1} \bar{v} \right) ^2 |(av_x)_x|^2 \leq  C\|(u,v,h)\|_{_E}^{2}  \|(\bar{u},\bar{v},\bar{h})\|_{_E}^{2},
		\end{equation*}
		for any $(u,v,h),(\bar{u},\bar{v},\bar{h})\in E$.
	\end{cor}

	\begin{proof}
		Take $(u,v,h),(\bar{u},\bar{v},\bar{h})\in E$ and let $M<0$ be the constant mentioned in Lemma $\ref{supremo}$. Since $\rho_{0}^{2} \rho_{\ast}^{-2} = \zeta^6 $ and $\tau^6 \leq \frac{45}{4M^6} e^{\frac{-2M}{m}}$, applying Proposition $\ref{ad2}$ and Lemma $\ref{supremo}$, we obtain

		\begin{align*}
		\displaystyle \int_0^T\int_0^1 \rho_{0}^{2} \left( \int_{0}^{1} \bar{u} \right) ^2 |(au_x)_x|^2
		&\leq \sup_{t\in[0,T]} \left\lbrace e^{\frac{-2M}{m}} \left(  \int_{0}^{1} \bar{u} \right) ^2 \right\rbrace
		\int_0^T\int_0^1 \rho_{\ast}^{2} |(au_x)_x|^2 \\
		&\leq C\|(u,v,h)\|_{E}^{2}  \|(\bar{u},\bar{v},\bar{h})\|_{E}^{2} .
		\end{align*}
		Analogously, a similar estimate also holds to $\int_0^T\int_0^1 \rho_{0}^{2} \left( \int_{0}^{1} \bar{v} \right) ^2 |(av_x)_x|^2$.
	\end{proof}

	\begin{prop}\label{2parts}
		The mapping $H: E\to F$ has the following properties:
		\begin{itemize}
			\item[(a)] $H$ is well defined;
			\item[(b)] For each $(u,v,h)\in E$, let us define $\bar{f}_{1}^{i} = D_i f_{1}(t,x,u,v)$ and $\bar{f}_{2}^{i} = D_i f_{2}(t,x,u,v)$, with $i=3,4$. Then, the linear mapping $T:E\longrightarrow G$ and $S:E\longrightarrow G$, given by
			\begin{align*}
			\displaystyle T(\bar{u},\bar{v},\bar{h})
			&= \bar{u}_t
			-\ell_{1}^{\prime} \left( \int_{0}^{1} u\right)
			\left( \int_{0}^{1} \bar{u}\right)(au_x)_x
			-\ell_1 \left(\int_{0}^{1} u\right) (a\bar{u}_x)_x \\
			&\hspace{1cm}+\bar{f}_{1}^{3} \bar{u}
			+ \bar{f}_{1}^{4} \bar{v}
			-\bar{h}\chi_{\omega},
			\end{align*}
			and
			\begin{align*}
			\displaystyle S(\bar{u},\bar{v},\bar{h})
			&= \bar{v}_t
			-\ell_{2}^{\prime} \left( \int_{0}^{1} v\right)
			\left( \int_{0}^{1} \bar{v}\right)(av_x)_x
			-\ell_2 \left(\int_{0}^{1} v\right) (a\bar{v}_x)_x \\
			&\hspace{1cm}+\bar{f}_{2}^{3} \bar{u}
			+ \bar{f}_{2}^{4} \bar{v},
			\end{align*}
			are the Gateaux derivative of $H_1$ and $H_2$ at $(u,v,h)\in E$, respectively.
		\end{itemize}

		\begin{proof} \noindent
			\begin{itemize}
				\item[(a)] For each $(u,v,h)\in E$, we must check that $H(u,v,h)\in F$. Of course, $H_3 (u,v,h)=u(\cdot,0)\in H_{a}^{1} (0,1)$ and $H_4 (u,v,h)=v(\cdot,0)\in H_{a}^{1} (0,1)$. Besides, recalling the Assumptions $\eqref{hyp_a}$ and $\eqref{hyp_f}$, and using Corollary $\ref{conseq}$, we take
				\begin{align*}
				\displaystyle &\int_0^T\int_0^1 \rho_{0}^{2} |H_1(u,v,h)|^2
				=\int_0^T\int_0^1 \rho_{0}^{2} \left| u_t-\ell_1 \left(\int_0^1 u \right)(au_x)_x+f_1(t,x,u,v) -h\chi_\omega\right| ^2 \\
				&\hspace{0.5cm} \leq 3 \int_0^T\int_0^1 \rho_{0}^{2} |u_t - (au_x)_x -h\chi_{\omega}|^2
				+3\int_0^T\int_0^1 \rho_{0}^{2} \left| \ell_1 \left( \int_{0}^{1} u\right) -\ell_1 (0) \right| ^2 |(au_x)_x|^2 \\
				&\hspace{1cm}+ 3 \int_0^T\int_0^1 \rho_{0}^{2} |f_1 (t,x,u,v) - f_1 (t,x,0,0)|^2 \\
				&\hspace{0.5cm} \leq 3 \int_0^T\int_0^1 \rho_{0}^{2} |u_t - (au_x)_x -h\chi_{\omega}|^2
				+C\int_0^T\int_0^1 \rho_{0}^{2} \left( \int_{0}^{1} u\right) ^2 |(au_x)_x|^2 \\
				&\hspace{1cm}+C \int_0^T\int_0^1 \rho_{0}^{2} (|u|^2+|v|^2) \\
				&\hspace{0.5cm}\leq C(\|(u,v,h)\|_{E}^{2} +\|(u,v,h)\|_{E}^{4} )
				\end{align*}
				Therefore, $H_1 (u,v,h)\in G$ and, in a similar way, we also have $H_2 (u,v,h)\in G$.

				\item[(b)] In this part, fix $(u,v,h)\in E$. Thus, for any $(\bar{u},\bar{v},\bar{h})\in E$ and $\lambda \neq 0$, we take
				\begin{align*}
				\displaystyle &\frac{1}{\lambda} [H_1(u+\lambda \bar{u},v+\lambda \bar{v},h+\lambda \bar{h}) - H_1(u,v,h)] -T(\bar{u},\bar{v},\bar{h}) \\
				&\hspace{0.5cm} =-\left[ \frac{1}{\lambda} \left( \ell_1 \left( \int_{0}^{1} (u+\lambda \bar{u} )\right)  - \ell_1 \left( \int_{0}^{1} u\right)  \right)
				-\ell_{1}^{\prime} \left( \int_{0}^{1} u\right) \left( \int_{0}^{1} \bar{u} \right) \right] (au_x)_x \\
				&\hspace{1cm} -\left[ \ell_1 \left( \int_{0}^{1} (u+\lambda \bar{u} )\right)  - \ell_1 \left( \int_{0}^{1} u\right) \right] (a\bar{u} _x)_x \\
				&\hspace{1cm} +\left[ \left( \frac{f_1(t,x,u+\lambda \bar{u},v+\lambda \bar{v})-f_1 (t,x,u,v)}{\lambda} \right) - (\bar{f}_{1}^{3} \bar{u} + \bar{f}_{1}^{4} \bar{v} )  \right] \\
				&\hspace{0.5cm}:=A_{\lambda}+B_{\lambda}+C_{\lambda}.
				\end{align*}
				We will see that $A_{\lambda}, B_{\lambda} \emph{ and } C_{\lambda}$ converge to zero in $G$, as $\lambda \to 0$. Indeed, taking into account  $\eqref{hyp_a}$ and Mean Value Theorem, for each $\lambda \neq 0$, there exists
				$u_{\lambda}=u_{\lambda}(t)$ such that $u_{\lambda} \to \int_{0}^{1} u$ for any $t\in[0,T]$,
				\[
				\displaystyle \int_0^T\int_0^1 \rho_{0}^{2} |A_{\lambda}|^2 = \int_0^T\int_0^1 \rho_{0}^{2}
				\left|  \ell_{1}^{\prime} (u_\lambda)  - \ell_{1}^{\prime} \left( \int_{0}^{1} u\right) \right|^{2}
				\left( \int_{0}^{1} \bar{u} \right)^2 |(au_x)_x|^2
				\to 0
				\]
				and
				\begin{align*}
				\displaystyle \int_0^T\int_0^1 \rho_{0}^{2} |B_{\lambda}|^2
				&=\lambda^2 \int_0^T\int_0^1 \rho_{0}^{2} |\ell_{^1}^{\prime} (u_{\lambda})|^2 \left( \int_{0}^{1} \bar{u} \right)^2 |(au_x)_x|^2 \\
				&\leq C\lambda^2 \|(u,v,h)\|_{_E}^{2}  \|(\bar{u},\bar{v},\bar{h})\|_{_E}^{2} \to 0,
				\end{align*}
				as $\lambda\to 0$.
				On the other hand, for each $\lambda\neq 0$, we can apply again Mean Value and Lebesgue's Theorem, in order to obtain $w_{\lambda}=w_{\lambda} (t,x)$ satisfying: $w_{\lambda} \to (t,x,u,v)$, for any $(t,x)\in (0,T)\times (0,1)$, and
				\begin{align*}
				\displaystyle \int_0^T\int_0^1 \rho_{0}^{2} |C_{\lambda}|^2
				&=\int_0^T\int_0^1 \rho_{0}^{2} |[D_3 f_1 (w_{\lambda})\bar{u} + D_{4} f_1 (w_{\lambda})\bar{v}]-[\bar{f}_{1}^{3} \bar{u}+\bar{f}_{1}^{4} \bar{v}]|^2 \\
				&\hspace{-2cm}\leq \int_0^T\int_0^1 \rho_{0}^{2} |D_3 f_1 (w_{\lambda})-\bar{f}_{1}^{3}|^2 |\bar{u}|^2
				+\int_0^T\int_0^1 \rho_{0}^{2} |D_4 f_1 (w_{\lambda})-\bar{f}_{1}^{4}|^2 |\bar{v}|^2
				\to 0,
				\end{align*}
				as $\lambda \to 0$.

				As a consequence, $T$ is the Gateaux derivative of $H_1$ at $(u,v,h)\in E$. Likewise, $S$ is the Gateaux derivative of $H_2$ at $(u,v,h)\in E$.
			\end{itemize}
		\end{proof}

		\begin{prop}\label{c1}
			The mapping $H:E\to  F$ is continuously differentiable.
		\end{prop}

		\begin{proof}
			Clearly, $H_3, H_4 \in C^1 (E,H_{a}^{1})$. Now, take $(u,v,h)\in E$ and let $((u^n,v^n,h^n))_{n=1}^{\infty}$ be a sequence which converges to $(u,v,h)$ in $E$. For each $(\bar{u},\bar{v},\bar{h})\in \bar{B}_{1}(0)\subset E$, we have proved in Proposition $\ref{2parts}$ that
			\begin{multline*}
			\displaystyle H_{1}^{\prime} (u,v,h)(\bar{u},\bar{v},\bar{h})
			= \bar{u}_t
			-\ell_{1}^{\prime} \left( \int_{0}^{1} u \right)
			\left( \int_{0}^{1} \bar{u}\right)(au_{x})_x
			-\ell_1 \left(\int_{0}^{1} u \right) (a\bar{u}_x)_x \\
			+\bar{f}_{1}^{3} \bar{u}
			+ \bar{f}_{1}^{4} \bar{v}
			-\bar{h}\chi_{\omega}
			\end{multline*}
			and
			\begin{multline*}
			\displaystyle H_{1}^{\prime} (u^n,v^n,h^n)(\bar{u},\bar{v},\bar{h})
			= \bar{u}_t
			-\ell_{1}^{\prime} \left( \int_{0}^{1} u^n \right)
			\left( \int_{0}^{1} \bar{u}\right)(au_{x}^{n})_x
			-\ell_1 \left(\int_{0}^{1} u^n \right) (a\bar{u}_x)_x \\
			+D_3 f_1 (t,x,u^n,v^n) \bar{u}
			+ D_4 f_1 (t,x,u^n,v^n) \bar{v}
			-\bar{h}\chi_{\omega}.
			\end{multline*}
			Thus,
			\begin{align*}
			\displaystyle &(H^{\prime}_{1} (u^n,v^n,h^n) - H^{\prime}_{1} (u,v,h))(\bar{u},\bar{v},\bar{h})
			=-\ell_{1}^{\prime} \left( \int_{0}^{1} u^n\right) \left( \int_{0}^{1} \bar{u} \right) [a(u^n-u)_x]_x \\
			&\hspace{.5cm} -\left[ \ell_{1}^{\prime} \left( \int_{0}^{1} u^n\right)-\ell^{\prime}_{1} \left( \int_{0}^{1} u \right) \right] \left( \int_{0}^{1} \bar{u} \right)   (au_x)_x
			-\left[ \ell_{1} \left( \int_{0}^{1} u^n\right)-\ell_{1} \left( \int_{0}^{1} u \right) \right] (a\bar{u}_x)_x \\
			&\hspace{.5cm}+[D_3 f_1 (t,x,u^n,v^n)-D_3 f_1 (t,x,u,v)]\bar{u}
			+[D_4 f_1 (t,x,u^n,v^n)-D_4 f_1 (t,x,u,v)]\bar{v} \\
			&\hspace{.5cm}:=X_{1}^{n}+X_{2}^{n}+X_{3}^{n}+X_{4}^{n}+X_{5}^{n}.
			\end{align*}
			From assumption $\ref{hyp_a}$ and Corollary $\ref{conseq}$, we get
			\begin{align*}
			\displaystyle \int_0^T\int_0^1 \rho_{0}^{2} |X_{1}^{n}|^2 &
			\leq C\int_0^T\int_0^1 \rho_{0}^{2} \left( \int_{0}^{1} \bar{u} \right) ^2 |[a(u^n-u)_x]_x|^2
			\\
			& \leq C\|(u^n-u,v^n-v,h^n-h)\|_{E}^{2}
			\to 0
			\end{align*}
			and
			\begin{align*}
			\int_0^T\int_0^1 \rho_{0}^{2} |X_{3}^{n}|^2
			& \leq C\int_0^T\int_0^1 \rho_{0}^{2} \left( \int_{0}^{1} (u^n-u)\right)^2 |(a\bar{u}_x)_x|^2\\
			& \leq C\|(u^n-u,v^n-v,h^n-h)\|_{E}^{2}
			\to 0,
			\end{align*}
			as $n\to +\infty$.

			On the other hand, due to Lemma $\ref{supremo}$ and assumptions $\ref{hyp_a}$ and $\ref{hyp_f}$, we obtain
			\begin{align*}
			\int_0^T\int_0^1 \rho_{0}^{2} |X_{2}^{n}|^2 & \leq\sup_{t\in[0,T]} \left\lbrace e^{\frac{-2M}{m}} \left( \int_{0}^{1} \bar{u} \right)^2 \right\rbrace \int_0^T\int_0^1 \rho_{\ast}^{2}
			\left| \ell_{1}^{\prime} \left( \int_{0}^{1} u^n\right)-\ell^{\prime}_{1} \left( \int_{0}^{1} u \right) \right|^2 |(au_x)_x|^2 \\
			&\hspace{.5cm}\leq \int_0^T\int_0^1 \rho_{\ast}^{2}
			\left| \ell_{1}^{\prime} \left( \int_{0}^{1} u^n\right)-\ell^{\prime}_{1} \left( \int_{0}^{1} u \right) \right|^2 |(au_x)_x|^2 \to 0,
			\end{align*}
			and
			\begin{align*}
			\int_0^T\int_0^1 \rho_{0}^{2} |X_{4}^{n}|^2&
			\leq C\left( \int_0^T\int_0^1 \rho_{0}^{2} |D_3 f_1 (t,x,u^n,v^n) - D_3 f_1 (t,x,u,v)|^2 |\bar{u}|^2 \right)^{\frac{1}{2}}  \left( \int_0^T\int_0^1 \rho_{0}^{2} |\bar{u}|^2 \right)^{\frac{1}{2}} \\
			&\leq C\left( \int_0^T\int_0^1 \rho_{0}^{2} |D_3 f_1 (t,x,u^n,v^n) - D_3 f_1 (t,x,u,v)|^2 |\bar{u}|^2 \right)^{\frac{1}{2}} \to 0,
			\end{align*}
			as $n\to +\infty$, where we have also applied Lebegue's Theorem. Clearly, $\int_0^T\int_0^1 \rho_{0}^{2} |X_{5}^{n}|^2$ is similar to $\int_0^T\int_0^1 \rho_{0}^{2} |X_{4}^{n}|^2$ and we conclude that $H^{\prime}_{1}$ is a continuous mapping. Analogously, this conclusion remains valid to $H^{\prime}_{2}$. In this case, the proof is complete.
		\end{proof}
	\end{prop}

	\begin{proof}[\textbf{Proof of Theorem $\ref{main}$}]
		From Propositions $\ref{2parts}$ and $\ref{c1}$, we already know that the $H\in C^1(E,F)$. We state that $H^{\prime} (0,0,0):E\longrightarrow F$ is onto. In fact, consider $b_{11} (t,x)=D_3 f_1 (t,x,0,0)$, $b_{12} (t,x)=D_4 f_1 (t,x,0,0)$, $b_{21} (t,x)=D_3 f_2 (t,x,0,0)$ and $b_{22} (t,x)=D_4 f_2 (t,x,0,0)$ in ($\ref{pb-lin}$). Thus, given $(g_1,g_2,u_0,v_0) \in F$, we apply Theorem $\ref{linearcontrol}$ in order to obtain $(u,v,h)$ which solves ($\ref{pb-lin}$) and satisfies the relations in ($\ref{eq25}$). As a result, $(u,v,h)\in E$ and $H^{\prime}(0,0,0)(u,v,h)=(g_1, g_2,u_0,v_0)$, as we were supposed to check.

		Hence, by \textit{Lyusternik's Inverse Mapping Theorem} (Theorem \ref{lyusternik}) , there exist $r>0$ and a mapping $\tilde{H} :B_{r} (0) \subset F \longrightarrow E$ such that
		\[
		\displaystyle H(\tilde{H} (y)) = y \emph{  for each  } y\in B_{r}(0)\subset F .
		\]
		In particular, if $(\bar{u} _0,\bar{v} _0)\in H_{a}^{1} \times H_{a}^{1}$ and $\| (\bar{u}_0,\bar{v}_0)\|_{H_{a}^{1} \times H_{a}^{1}} < r$, we conclude that $(\bar{u},\bar{v},\bar{h})=\tilde{H} (0,0,\bar{u}_0,\bar{v}_0)\in E$ solves  $H(\bar{u},\bar{v},\bar{h}) = (0,0,\bar{u}_0,\bar{v}_0)$. Since $\int_0^T\int_0^1 \rho_{0}^{2} (|\bar{u}|^2+|\bar{v}|^2) <+\infty$, we get $\bar{u}(T,x)=\bar{v}(T,x)=0$ for any $x\in[0,1]$, following the result.
	\end{proof}

	\section{Some Additional Comments}\label{sec6}

	As a first comment, we note that, in assumption \ref{hyp_a}, we have taken a weak type of degeneracy and so that Dirichlet boundary conditions are required in \eqref{pb1}. However, if we had chosen strong type degeneracy, see \cite{alabau2006carleman}, \eqref{pb1} it would be treated with Neumann conditions. In this context, we believe that analogous results can be obtained.

	Another interesting question is concerned with global null controlability to \eqref{pb1}, which does not seem to be simple. Perhaps, this kind of result relies on a global inverse mapping theorem, see \cite{de1994global}, but much more refined estimates are necessary.

	Under some changes in the Lemma \ref{supremo} and following the arguments presented here, Theorem \ref{main} can also be obtained if we consider \eqref{pb1} with  the diffusion coefficients
	\[\left(\mu_1 \left(x,u\right)u_x \right)_x \text{ and } \left(\mu_1 \left(x,v\right)v_x \right)_x.\]

	Other important topics arrise from our current research:

	\begin{itemize}
		\item It would be very nice to obtain Theorem \ref{main} without imposing $\mu_1$ and $\mu_2$ have separated variables. Nevertheless, it is still an open problem.

		\item In the  system $\eqref{pb1}$, we can replace each nonlinearity  $f_i(t,x,u,v)$ by
		$f_i (t,x,u,v,u_x,v_x)$, with $i\in \{1,2\}$, in order to analyse whether it is possible to prove results about null controllability.

		\item Previously, in \cite{jrl2016}, we have obtained a local null controlability result for degenarate parabolic equations with nonlocal tems, which implies, throughout standard arguments, a local null boundary controllability result. However, the same fact can not be directly deduced for systems with a reduced number of controls, see \cite{ammar2011recent}. In other words, the boundary controllability of
		\begin{equation*}
		\begin{cases}
		u_t-\left(\mu_1 \left(x,\int_{0}^{1} u\right)u_x \right)_x+f_1(t,x,u,v) =0, & (t,x)\in (0,T)\times (0,1),\\
		v_t-\left(\mu_2 \left(x,\int_{0}^{1} v \right)v_x \right)_x+f_2(t,x,u,v)=0, & (t,x) \in (0,T)\times (0,1),\\
		u(t,1)=v(t,0)=v(t,1)=0,\ u(t,0)= h(t) \  &  t\in (0,T),\\
		u(0,x)=u_0(x) \ \ \mbox{ and }\ \ v(0,x)=v_0(x), & x\in  (0,1).
		\end{cases}
		\end{equation*}
		is a very interesting unknown issue.

	\end{itemize}

redpajama_set_name: RedPajamaArXiv