text:
\section{Introduction}

In this article we prove first and second order optimality conditions for a control problem subject to the Fokker-Planck equation
\begin{equation}
\label{FPintro}
\partial_t \rho - \nu \Delta \rho - {\rm div } \big(\rho B[u]\big) = 0,
\end{equation}
where $u$ is a time-dependent possibly multidimensional control, and $B[u](x)=c(x) + b(x)\otimes u$ is defined by given vector functions $c$ and $b$.

Fokker-Planck equations arise in many situations in which a large number of agents is involved. In particular,
these equations are present in models of mass opinion dynamics \cite{AlbiPareschiZanella2014}, tumor growth \cite{BoseTrimper2009,Herty2012}, bird flocks movement \cite{DuanFornasierToscani2010} and various biological events \cite{SchienbeinGruler1993,Chavanis2008nonlinear}, among others. The recent survey \cite{Furioli2017} describes several applications of the Fokker-Planck equation to different socio-economic phenomena. It is worth mentioning that some of these articles already analyze control problems associated to their models, as {\em e.g.}  \cite{Herty2012,AlbiPareschiZanella2014}.

Another motivation for investigating  optimal control problems governed by Fok\-ker-Plack equations comes from Mean Field Games (MFG) theory \cite{LasryLions2007,GomesSaude2014}. MFG' systems consist in a Hamilton-Jacobi equation coupled with a Fokker-Planck one. It has been observed \cite{LasryLions2007,cardaliaguet2010notes,Ryzhik2018} that the MFG' system is an optimality condition for the control of Fokker-Planck equations, but the rigorous proof of this fact remains an open problem, since the theory of optimal control of this equation is far from being complete.

Optimal control problems governed by  \eqref{FPintro} have been recently studied in a few articles.  A problem with one-dimensional piecewise constant controls is treated in \cite{AnnunziatoBorzi2010}, where the authors propose control strategies based on a receding horizon model predictive control framework.
That work was extended later in \cite{AnnunziatoBorzi2013} to the case of multidimensional constant control. A Fokker-Planck equation with a space-dependent time-independent control is investigated in \cite{RoyAnnunziatoBorzi2016fokker}. There, the authors prove existence of an optimal control, derive a first-order optimality system, and propose a numerical scheme along with simulations.
For a more general framework, with time- and space-dependent multidimensional control, an existence result and first order optimality conditions are given in \cite{FleigGuglielmi2017}.
See also \cite{BreitenKunischPfeiffer2018}, where the authors provide results on stabilization of \eqref{FPintro} through linearization (more details in Remark \ref{RemarkBiblio}).

Optimization problems associated to \eqref{FPintro} belong to the class of {\em bilinear optimal control}. This framework has been considered in {\em e.g.} \cite{CasasWachsmuthWachsmuth2018} for elliptic equations. On the other hand, the work \cite{AronnaBonnansKroener2018} dealt with infinite dimensional bilinear dynamical systems, but their results do not apply here.

In this paper we provide first and second order optimality conditions for an optimal control problem associated to \eqref{FPintro}, under quite mild regularity assumptions on the data functions and the spatial domain
of the state equation. The cost functional is of tracking type and includes a {\em Tikhonov regularizing term} on the control.  The control is only time-dependent. We obtain our second order necessary and sufficient conditions by application of  results proved in \cite{CasasTroeltzsch2012} in an abstract framework.

Our paper is organized as follows. In Section \ref{S2}\textcolor{cyan}, we present the equation, the basic assumptions, show well-posedness and some other properties of the state equation. The optimal control problem is introduced in Section \ref{S3}, where existence of optimal solution is proved. In Section \ref{S4}, we discuss properties of the control-to-state mapping, while Section \ref{S5} presents the adjoint system and first order optimality conditions. Finally, Section \ref{S6} is devoted to the second order analysis. Our main theorems are contained there.

\subsection*{Notation} \label{S1}
Given a real interval $[0,T]$, a normed space $X$ and $p\in [1,\infty],$  we let $L^p(0,T;X)$ denote the Lebesgue space of $L^p$-functions and  write $C([0,T];X)$ for the space of continuous functions, both with domain $[0,T]$ and values in $X.$
When $X=\mathbb{R}$ we just write $L^p(0,T)$ and, for any $m$, we let $\|\cdot\|_{p}$ denote the norm in $L^p(0,T;\mathbb{R}^m)$.
 Analogously, for a set $\Omega \subseteq \mathbb{R}^n,$ $L^p(\Omega;\mathbb{R}^m)$ is the Lebesgue space of $L^p$-functions with domain $\Omega$ and values in $\mathbb{R}^m$.  We omit $\mathbb{R}^m$ when $m=1$ and the values range in $\mathbb{R},$  and we  use the  form $\|\cdot\|_{L^p(\Omega)^m}$ to denote the norm in $L^p(\Omega;\mathbb{R}^n)$.

Throughout the article, we consider the real Hilbert spaces $L^2(\Omega)$, $H^1(\Omega)$ and $H^1(\Omega)^*.$ We let $(\cdot,\cdot)$ denote the scalar product in $L^2(\Omega)$ and  $\langle \cdot,\cdot\rangle$ the dual pairing between $H^1(\Omega)^*$ and $H^1(\Omega)$. Other scalar products and pairings will be distinguished by specifying the spaces as subindexes. For two vectors $u, v \in \mathbb{R}^n$, the result of the componentwise multiplication $u\otimes v$ is defined by
the vector $ w \in \mathbb{R}^n$ with components $w_i = u_i v_i$, for $i = 1,\ldots,n$.

\section{The controlled Fokker-Planck equation} \label{S2}

We consider the Fokker-Planck equation with initial and boundary conditions given by
\begin{align}
\label{FP}
\partial_t \rho(x,t) - \nu \Delta \rho(x,t) - {\rm div } \big(\rho(x,t) B[u(t)](x)\big) &= 0\quad \text{in } Q,\\
\label{FPinit} \rho(x,0) &= \rho_0(x)\quad \text{in } \Omega,\\
\label{PBNeumann}\Big(\nu \nabla \rho(x,t) + \rho(x,t) B[u(t)](x)\Big)\cdot n(x) &=0\quad \text{on } \Sigma,
\end{align}
where $\nu>0,$ $\rho_0 \in L^2(\Omega)$ stands for  an initial probability distribution; therefore, it is assumed to be nonnegative and to satisfy
\[
\int_\Omega \rho_0(x)dx =1;
\]
$\Omega \subset \mathbb{R}^n$ is a bounded domain with Lipschitz boundary $\Gamma:= \partial \Omega,$ and we set $\Sigma := \Gamma \times (0,T),$ $Q:=\Omega \times (0,T).$
The control is $u=(u_1,\dots,u_{n}) \in L^\infty(0,T;\mathbb{R}^{n})$,
and the function $B \colon \mathbb{R}^n \times \Omega \to \mathbb{R}^n$ is given by
\[
{B}[u](x) := c(x) + {b}(x)\otimes {u},
\]
where $c,b \in L^\infty(\Omega;\mathbb{R}^n)$ are fixed.
In \eqref{FP}, the differential operators $\Delta$ and ${{\rm div} }$ act only with respect to the spatial coordinate $x$. For definition and basic properties of the ${\rm div}$-operator we refer to \cite{GiraultRaviart2012}. In \eqref{PBNeumann}, $n$ denotes the outward normal unit vector on $\Gamma$.

\begin{remark}
\label{RemarkBiblio} Let us compare this equation with others considered in the literature. We  mention that
Breiten {\em et al.} \cite{BreitenKunischPfeiffer2018} assume $B$ to take the form $\nabla_x V$ for a potential $V$ given by
\[
V(x,t) = G(x) + \gamma(x) u(t)
\]
for a scalar control $u$.  They  study the stabilization problem (in infinite horizon) rather than the associated optimal control problem.

Fleig and Guglielmi \cite{FleigGuglielmi2017} considered a distributed multidimensional control and adopt a vanishing condition on $\Sigma$. In this case, the measure $\int_\Omega \rho (t,x)dx $ is not preserved. For that problem, they showed first order optimality conditions.

Equation \eqref{FP}, for a function $u$ that depends both on time and space, appears in second order Mean Field Games: see {\em e.g.} Gomes and Sa\'ude \cite[Theorem 2]{GomesSaude2014}, Lasry and Lions \cite[Section 2.6]{LasryLions2007}.
\end{remark}

\subsection{Existence and uniqueness of the solution of the Fokker-Planck equation}

We next make a formal derivation of the weak formulation of \eqref{FP}-\eqref{PBNeumann}.
Consider $\varphi$ sufficiently smooth. From \eqref{FP} we get
\begin{equation}
 \int_\Omega \big(\partial_t \rho  - \nu\Delta\rho - {\rm div}(\rho B[u]) \big) \varphi dx  = 0.
\end{equation}
We have
\begin{equation}
\label{var2}
\begin{split}
- \int_\Omega  \nu(\Delta\rho) \varphi
& =  -\nu  \int_{\Gamma} \varphi \nabla \rho \cdot n ds  + \nu \int_\Omega \nabla \rho \cdot \nabla\varphi dx ,
\end{split}
\end{equation}
\begin{equation}
\label{var3}
\begin{split}
- \int_\Omega{\rm div} (\rho B[u])\varphi dx
& = - \int_\Gamma \varphi \rho B[u] \cdot n ds  +   \int_\Omega \rho B[u] \cdot \nabla \varphi dx .
\end{split}
\end{equation}
Adding  \eqref{var2} and \eqref{var3}, in view of the boundary condition \eqref{PBNeumann} we get
\[
\int_\Omega \partial_t \rho\varphi dx + a[u(t)](\rho,\varphi) =0,
\]
where, for each $u\in \mathbb{R}^n,$  $a[u](\cdot,\cdot)$ is a bilinear mapping that to each pair $\psi,\varphi \in H^1(\Omega)$ associates the value
\begin{equation}
\label{adef}
a[u](\psi,\varphi):= \int_\Omega\big(\nu \nabla \psi + \psi B[u]\big) \cdot \nabla \varphi dx.
\end{equation}
We will work in the space
\begin{equation}
\label{Wdef}
W(0,T) :=\left\{ \rho \in L^2(0,T;H^1(\Omega)) : \partial_t \rho\in L^2(0,T;H^{1}(\Omega)^*)\right\},
\end{equation}
equipped with the norm
\[
\|\rho\|_{W(0,T)} := \left(\|\rho\|^2_{L^2(0,T;H^1(\Omega))} + \left\|\partial_t \rho\right\|_{L^2(0,T;H^1(\Omega)^*)}^2 \right)^{1/2}.
\]
It is known that $W(0,T)$ is a Hilbert space with the scalar product
\[
(\rho,\varphi)_{W(0,T)} := \int_0^T(\rho(t),\varphi(t))_{H^1(\Omega)} dt +
\int_0^T \big( \partial_t \rho(t) , \partial_t \varphi(t) \big)_{H^1(\Omega)^*} dt,
\]
which induces the norm $\|\cdot\|_{W(0,T)}.$ We shall also recall the following continuous embedding (see {\em e.g.} Chipot \cite[Theorem 11.4]{Chipot2012} or Dautray-Lions \cite[Theorem 1, page 473]{DautrayLions1992})
$$
W(0,T) \hookrightarrow C \big( [0,T];L^2(\Omega) \big)
$$
that will be of use throughout this article.

The weak formulation of \eqref{FP} can be rewritten as to find $\rho \in W(0,T)$ such that
\begin{align}
\label{varfor1} \frac{d}{dt} (\rho (\cdot), \varphi) + a[u(\cdot)](\rho(\cdot),\varphi) &=0\quad \text{on } \mathcal{D}'(0,T)\text{ for all } \varphi \in H^1(\Omega),\\
\label{initcond} \rho(0) &= \rho_0\quad \text{in } L^2(\Omega).
\end{align}

For convenience, let us consider a general right-hand side $f\in L^2(0,T;H^1(\Omega)^*)$ in \eqref{varfor1} and study the existence and uniqueness of the solution of the equation
 \begin{equation}
\label{varforf} \frac{d}{dt} ( \rho(\cdot), \varphi) + a[u(\cdot)](\rho(\cdot),\varphi) =\langle f(\cdot),\varphi\rangle\quad \text{on } \mathcal{D}'(0,T)\text{ for all } \varphi \in H^1(\Omega),
\end{equation}
with initial condition \eqref{initcond}.

\begin{remark}
Note that \eqref{varforf} can also be expressed as
\begin{equation}
\partial_t \rho+ a[u](\rho,\cdot ) = f \quad \text{in } L^2(0,T;H^1(\Omega)^*).
\end{equation}
\end{remark}

\begin{definition}
We call a function $\rho \in W(0,T)$ a {\em weak solution} of the Fokker-Planck equation \eqref{FP}, if it satisfies \eqref{varfor1}-\eqref{initcond}. The definition is formulated analogously for \eqref{varforf}-\eqref{initcond}, where the r.h.s. of the differential equation in \eqref{FP} was replaced by an arbitrary element $f$ of the space $L^2(0,T;H^1(\Omega)^*).$
\end{definition}

\begin{remark}
For practical matters, we mention the equivalent variational formulation of \eqref{FP}-\eqref{PBNeumann} according to  Ladyzhenskaya \cite{Ladyzhenskaya_etal1968} that is as  follows: $\rho$ is a function in $W_2^{1,0}(\Omega)$ such that
\begin{multline}
\label{varfor}
\iint_Q \Big\{ -\rho \partial_t \varphi+  \big( \nu\nabla \rho + \rho B[u]\big)\cdot\nabla \varphi \Big\} dx\, dt =\int_\Omega \rho_0 \varphi(\cdot ,0) d x \\
\text{for all }\varphi \in W^{1,1}_2(Q) \text{ with } \varphi(\cdot,T)=0.
\end{multline}
\end{remark}

\if{
\begin{proof}
From \eqref{FP} we get
\begin{equation}
\int_0^T \int_\Omega \big(\partial_t \rho  - \nu\Delta\rho - {\rm div}(\rho B) \big) \varphi dx dt = 0.
\end{equation}
We have
\begin{equation}
\label{var1}
\int_0^T \int_\Omega \partial_t \rho\varphi = -\int_0^T \int_\Omega \rho \dot\varphi dxdt -\int_\Omega \rho_0(x) \varphi(x,0) dx,
\end{equation}
\begin{equation}
\label{var2}
\begin{split}
-\int_0^T \int_\Omega  \nu\Delta\rho \varphi &= -\nu \int_0^T \int_\Omega {\rm div} (\varphi \nabla \rho) dx dt + \nu\int_0^T \int_\Omega \nabla \rho \cdot \nabla\varphi dx dt\\
& =  -\nu \int_0^T \int_{\Gamma} \varphi \nabla \rho \cdot n ds ds + \nu\int_0^T \int_\Omega \nabla \rho \cdot \nabla\varphi dx dt,
\end{split}
\end{equation}
\begin{equation}
\label{var3}
\begin{split}
-\int_0^T \int_\Omega{\rm div} (\rho B)\varphi dx dt & = -\int_0^T \int_\Omega {\rm div}(\rho B\varphi)dx dt + \int_0^T \int_\Omega \rho B \cdot \nabla \varphi dx dt \\
& = -\int_0^T\int_\Omega \varphi \rho B \cdot n ds dt +  \int_0^T \int_\Omega \rho B \cdot \nabla \varphi dx dt.
\end{split}
\end{equation}
Adding up \eqref{var1}, \eqref{var2} and \eqref{var3} yields the desired result.
\end{proof}
}\fi

\begin{theorem}
\label{ThmExistSol}
Given $\rho_0\in L^2(\Omega),$ $u\in L^\infty(0,T;\cR^n)$,
and $f\in L^2(0,T;H^{1}(\Omega)^*),$
there exists a unique weak solution $\rho$ of \eqref{varforf}-\eqref{initcond}
and it belongs to the space $W(0,T)$.
\end{theorem}

\begin{proof}
We apply Lions-Magenes \cite[Theorem 4.1 and Remark 4.3, pp. 238-239]{LioMag68a}; see also Chipot \cite[Theorem 11.7]{Chipot2012}.

We first verify the assumptions of continuity and coercivity of the bilinear form $a[u]$ necessary to apply \cite[Theorem 4.1 and Remark 4.3, pages 238-9]{LioMag68a}. For $\rho,\varphi \in H^1(\Omega),$ we have
\begin{equation}
\label{estarhophi}
\begin{split}
\big|a[u(t)](\rho,\varphi)\big| &\leq \left| \nu \int_\Omega \nabla \rho \cdot \nabla \varphi dx \right|  + \left| \int_\Omega  \rho B[u(t)]\cdot \nabla \varphi dx \right| \\
&\leq \nu\|\nabla \rho\|_{L^2(\Omega){^n}} \|\nabla \varphi\|_{L^2(\Omega){^n}} + \|B[u]\|_{L^\infty(Q){^n}} \|\rho\|_{L^2(\Omega)}\|\nabla \varphi\|_{L^2(\Omega)^n}\\
&\leq C\|\rho\|_{H^1(\Omega)} \|\varphi\|_{H^1(\Omega)},
\end{split}
\end{equation}
where $C := \nu +  \|B[u]\|_{L^\infty(Q)^n}$. From  the latter estimate we also obtain, for any positive $M,$ By Young's inequality, we also obtain
\begin{equation*}
 \left| \int_\Omega  \rho B[u]\cdot \nabla \varphi ds \right|
 \leq \|B[u]\|_{L^\infty(Q)^n}  \left( \frac{M}{2}\|\rho\|_{L^2(\Omega)}^2 + \frac{1}{{2}M} \|\nabla\varphi\|^2_{L^2(\Omega)^n}\right).
 \end{equation*}
 Thus
\begin{equation}
\label{acoercive}
\begin{split}
a[u(t)]&(\varphi,\varphi)
= \nu  \int_\Omega |\nabla\varphi|^2 dx +  \int_\Omega  \varphi B[u(t)]\cdot \nabla \varphi dx\\
&\geq \nu \|\varphi\|^2_{H^1(\Omega)}  -\nu \|\varphi\|^2_{L^2(\Omega)}   -\|B[u]\|_{L^\infty(Q)^n}\left( \frac{M}{ 2}\| \varphi \|_{L^2(\Omega)}^2 + \frac{1}{{2}M} \|\nabla\varphi\|^2_{L^2(\Omega)^n}\right) \\
&\geq \gamma \|\varphi\|^2_{H^1(\Omega)} - \lambda \|\varphi\|_{L^2(\Omega)}^2,
\end{split}
\end{equation}
if $M$ is chosen such that  $\displaystyle\gamma := \nu -\frac{\|B[u]\|_{L^\infty(Q)^n}}{2M}$ is positive, and we set $\lambda := \nu+ \displaystyle \frac{M}{2}\, \|B[u]\|_{L^\infty(Q)^n.} $

We conclude from \eqref{estarhophi} and \eqref{acoercive} that the hypotheses of boundedness and coercivity of $a[u]$  in Lions-Magenes \cite[Theorem 4.1 and Remark 4.3, pp. 238-239]{LioMag68a} are satisfied and, therefore, their result can be applied.
\end{proof}

\begin{proposition}
Given $\rho_0,u,f$ as in Theorem \ref{ThmExistSol}, the unique weak solution $\rho$ of \eqref{FP}-\eqref{PBNeumann} satisfies
\begin{equation}
\label{intrho1}
\int_\Omega \rho(x,t) dx = \int_\Omega \rho_0(x) ds\quad \text{for a.a. } t\in (0,T).
\end{equation}
\end{proposition}

\begin{proof}
This follows from \eqref{varfor1} by choosing $\varphi \equiv 1.$

\if{If not, the following proof holds.

To obtain the same result from \eqref{varfor}, we can do as follows. Fix $s\in (0,T).$ Consider a test function $\varphi:=\varphi(t)$, with $\varphi \in W^{1,1}_2(Q)$ and such that
\begin{equation}
\varphi(t) =
\left\{
\begin{split}
1,\quad &t\in [0,s],\\
-\frac{1}{\epsilon}(t-\epsilon-s),\quad & t\in (s,s+\epsilon],\\
0,\quad & t\in (s+\epsilon,T].
\end{split}
\right.
\end{equation}
Then, the variational formulation \eqref{varfor} gives
\begin{equation}
-\int_0^T\int_\Omega \rho(x,t) \dot\varphi(t) dx dt = \int_\Omega \varphi_0(x)\varphi(0)dx,
\end{equation}
which implies
\begin{equation}
\int_s^{s+\epsilon} \frac{1}{\epsilon} \int_\Omega \rho(x,t) dx dt = \int_\Omega \varphi_0(x) dx.
\end{equation}
Setting $R(t):= \int_\Omega \rho(x,t) dx,$ we can rewrite latter display as
\begin{equation}
\frac{1}{\epsilon} \int_s^{s+\epsilon} R(t) dt = \int_\Omega \varphi_0(x) dx.
\end{equation}
For a.e. $s\in (0,T)$, the l.h.s. of latter equation converges to $R(s),$ thus we get \eqref{intrho1} as desired.
}\fi
\end{proof}

\begin{theorem}[Estimates]
\label{ThmEstimates}
Given $\rho_0,u,f$ as in Theorem \ref{ThmExistSol}, the unique weak solution $\rho$ of \eqref{FP}-\eqref{PBNeumann} obeys the following estimates:
\begin{equation}
\label{estH}
\|\rho\|^2_{C([0,T];L^2(\Omega))} \leq C_0 \Big(\|\rho_0\|^2_{L^2(\Omega)} +\|f\|_{L^2(0,T;H^1(\Omega)^*)}^2\Big),
\end{equation}
\begin{equation}
\label{estV}
\|\rho\|^2_{L^2(0,T;H^1(\Omega))} \leq C_1 \Big(\|\rho_0\|^2_{L^2(\Omega)} +\|f\|_{L^2(0,T;H^1(\Omega)^*)}^2\Big),
\end{equation}
\begin{equation}
\label{estrhodot}
\left\|\partial_t \rho\right\|_{L^2(0,T;H^1(\Omega)^*)}^2 \leq C_2 \Big(\|\rho_0\|^2_{L^2(\Omega)} + \|f\|_{L^2(0,T;H^1(\Omega)^*)}^2\Big),
\end{equation}
the constants $C_i$ for $i=0,1,2$ depending continuously on $\|u\|_\infty.$
\end{theorem}

\begin{proof}
For a.a. $t\in [0,T]$, we have
\begin{equation}
\label{varforest}
\frac{d}{dt}\frac12 \| \rho(t)\|_{L^2(\Omega)}^2 + a[u(t)](\rho(t),\rho(t))=\langle f(t) , \rho(t) \rangle.
\end{equation}
Thus,
\begin{equation*}
\frac{d}{dt} \frac12 \| \rho(t)\|_{L^2(\Omega)}^2 + \int_\Omega \nu|\nabla\rho(t)|^2 dx + \int_\Omega \rho(t)B[u(t)]\cdot\nabla\rho(t)dx = \langle f(t),\rho(t) \rangle,
\end{equation*}
which implies
\begin{equation}
\label{bigestimate}
\begin{split}
\frac12 \frac{d}{dt} & \|\rho(t)\|^2_{L^2(\Omega)} + \nu \|\nabla\rho(t)\|_{L^2(\Omega)^n}^{2}\\
&\leq \|B[u(t)]\|_{L^\infty(\Omega)^n}  \|\rho(t)\|_{L^2(\Omega)} \|\nabla\rho(t)\|_{L^2(\Omega)^n} + \|f(t)\|_{H^1(\Omega)^*} \|\rho(t)\|_{H^1(\Omega)}  \\
&\leq \|B[u(t)]\|_{L^\infty(\Omega)^n} \frac12 \left( M \|\rho(t)\|_{L^2(\Omega)}^2+\frac1M\|\nabla\rho(t)\|_{L^2(\Omega)^n}^2\right) \\
&\quad+
{\frac{N}{2}} \|\rho(t)\|_{H^1(\Omega)}^2+{\frac{1}{2N}}\|f(t)\|_{H^1(\Omega)^*}^2\\
&={\frac12} \left(M \|B[u(t)]\|_{L^\infty(\Omega)^n} +N\right) \|\rho(t)\|_{L^2(\Omega)}^2
\\
&\quad+ {\frac12}\left(\frac{\|B[u(t)]\|_{L^\infty(\Omega)^n}}{M} + N \right)\|\nabla\rho(t)\|_{L^2(\Omega)^n}^2 + {\frac{1}{2N}}\|f(t)\|_{H^1(\Omega)^*}^2.
\end{split}
\end{equation}
The latter inequality holds for any positive $M$ and $N.$
Let us set
\[
M := \frac{2}{\nu} \|B[u]\|_{L^\infty(Q)^n}, \quad N :=  \nu,
\]
Then, from \eqref{bigestimate}, we get
\[
\frac{d}{dt} \|\rho(t)\|^2_{L^2(\Omega)} \leq
{\frac12}\left(M \|B[u(t)]\|_{L^\infty(\Omega)^n} +N\right)\|\rho(t)\|_{L^2(\Omega)}^2+{\frac{1}{2N}}\|f(t)\|_{H^1(\Omega)^*}^2.
\]
By Gronwall's Lemma applied to the latter inequality, we deduce that
\[
\sup_{t\in [0,T]}   \|\rho(t)\|^2_{L^2(\Omega)} \leq e^{D_0} \Big(\|\rho_0\|^2_{L^2(\Omega)}+{\frac{1}{2N}} \|f\|_{L^2(0,T;H^1(\Omega)^*)}^2\Big),
\]
where
\[D_0:=  \frac{1}{2} \int_0^T\big(M\|B[u(t)]\|_{L^\infty(\Omega)^n}+{\frac{1}{N}}\big)dt.
\]
This yields \eqref{estH}.

Now, going back to \eqref{bigestimate},  integrating both sides of it with respect to time, and  setting
$K:=\nu-\frac{1}{2}\left(\frac{\|B[u]\|_{L^\infty(Q)}}{M}+N\right) =\frac{\nu}{4} >0,$
we infer
\begin{equation*}
\begin{split}
K\int_0^T & \|\nabla\rho(t)\|_{L^2(\Omega)^n}^2 dt \leq
{\frac12}\|\rho(0)\|^2_{L^2(\Omega)} - {\frac12}\|\rho(T)\|^2_{L^2(\Omega)} \\
&+
\frac{1}{2}\int_0^T \left( M \|B[u(t)]\|_{L^\infty(\Omega)^n} + N\right)  \|\rho(t)\|_{L^2(\Omega)}^2 dt + \frac{1}{2N} \|f\|_{L^2(0,T;H^1(\Omega)^*)}^2.
\end{split}
\end{equation*}
This implies \eqref{estV}.

To get \eqref{estrhodot}, let us take $\varphi \in H^1(\Omega)$ with $\|\varphi\|_{H^1(\Omega)} = 1$,
and write
\[
\left\langle \partial_t \rho(t),\varphi \right\rangle =- a[u(t)](\rho(t),\varphi) + \langle f(t),\varphi\rangle.
\]
Hence
\[
\left|\left\langle \partial_t \rho(t),\varphi \right\rangle \right| \leq K \left(\|\rho(t)\|_{H^1(\Omega)} + \|f(t)\|_{H^1(\Omega)^*} \right) \|\varphi\|_{H^1(\Omega)}
\]
that for a.a. $t \in (0,T)$ implies
\[
\left\|\partial_t \rho(t)\right\|_{H^1(\Omega)^*}\leq K \left(\|\rho(t)\|_{H^1(\Omega)} + \|f(t)\|_{H^1(\Omega)^*} \right).
\]
This yields the estimate \eqref{estrhodot}.
\end{proof}

\begin{remark}
The fact that the constants $C_i,$ for $i=0,1,2,$ depend on $\|u\|_\infty$ is not written explicitly in the estimates \eqref{estH}-\eqref{estrhodot} since it will not play a role in the forthcoming estimates, invoking that the set  $\mathcal{U}_{\rm ad}$ of admissible controls (given in \eqref{uad}) is bounded in $L^\infty(0,T;\mathbb{R}^n).$
\end{remark}

\begin{proposition}[Weak maximum principle]
\label{Th11.9-Chipot}
Fix $u\in L^\infty(0,T;\cR^n),$ $\rho_0^1,\rho_0^2 \in L^2(\Omega)$ and $f^1,f^2 \in L^2(0,T;H^1(\Omega)^*).$ Let $\rho_1,\rho_2$ be the weak solutions of \eqref{FP}-\eqref{PBNeumann} with control $u,$  associated to $\rho_0^1, f^1$ and $\rho_0^2,f^2,$ respectively. Suppose further that
\[
\rho_0^1 \leq \rho_0^2\quad \text{and} \quad f^1 \leq f^2,
\]
where the second inequality means $\langle f_1,\varphi \rangle \leq \langle f_2,\varphi \rangle$ for every $\varphi \in L^2(0,T;H^1(\Omega))$ such that  $\varphi \geq 0$ a.e. in $Q.$
Then we have that
\[
\rho_1\leq \rho_2\quad \text{a.e. in }Q.
\]
\end{proposition}

\begin{proof} This result follows straightforwardly from Chipot \cite[Theorem 11.9, p. 202]{Chipot2012}.

\if{This proof is taken from Chipot {\cite[Theorem 11.9, p. 202]{Chipot2012}}.
We have
\begin{equation}
\frac{d}{dt} \big( (\rho_1-\rho_2)(t),\varphi \big) + a[u](t;(\rho_1-\rho_2)(t),\varphi) = \left\langle \textcolor{red}{(f^1-f^2)(t)},\varphi \right\rangle,\quad
\end{equation}
\textcolor{red}{for all  $\varphi \in H^1(\Omega)$.}
We know, by Riesz Representation Theorem, that we can write
\begin{equation}
a[u](t;\rho,\cdot) = (A(t)\rho)(\cdot),
\end{equation}
\textcolor{red}{where $A(t)$ is ...}

We get
\begin{equation}
\frac{d}{dt} (\rho_1-\rho_2)\textcolor{red}{(t)} +A(t)(\rho_1-\rho_2)\textcolor{red}{(t)} =  \textcolor{red}{(f^1-f^2)(t)},\quad \text{in } L^2(0,T;H^1(\Omega)^*).
\end{equation}
In view of \cite[Theorem 2.8, p28]{Chipot2012} we have that $(\rho_1-\rho_2)^+ \in L^2(0,T;H^1(\Omega)),$ then
\begin{equation}
\langle (\rho_1-\rho_2)_t,(\rho_1-\rho_2)^+ \rangle + a[u](\deleted{t}\textcolor{red}{\cdot};\rho_1-\rho_2, (\rho_1-\rho_2)^+ ) = \langle f^1-f^2,(\rho_1-\rho_2)^+ \rangle.
\end{equation}
By \cite[Lemma 11.2]{Chipot2012},  $\displaystyle\frac{d}{dt} \|(\rho_1-\rho_2)^+\|^2_{L^2(\Omega)} = \frac{d}{dt} \langle \rho_1-\rho_2,(\rho_1-\rho_2)^+ \rangle,$ thus
\begin{equation}
\label{ddtrhoi}
\frac{d}{dt} \|(\rho_1-\rho_2)^+\|^2_{L^2(\Omega)}  + a[u](\deleted{t} \textcolor{red}{\cdot};(\rho_1-\rho_2)^+, (\rho_1-\rho_2)^+ ) = \langle f^1-f^2,(\rho_1-\rho_2)^+ \rangle,
\end{equation}
where we used $a[u](t;\rho_1-\rho_2, (\rho_1-\rho_2)^+ ) =a[u](t;(\rho_1-\rho_2)^+, (\rho_1-\rho_2)^+ )$ ({\bf XXX: Justify this}).
Using the property \eqref{acoercive} and $f^1\leq f^2,$ we obtain from \eqref{ddtrhoi} that
\begin{equation}
\frac{d}{dt} \|(\rho_1-\rho_2)^+\|^2_{L^2(\Omega)}  - \lambda \|(\rho_1-\rho_2)^+\|^2_{L^2(\Omega)} \leq 0.
\end{equation}
Setting $z(t) : = \|(\rho_1-\rho_2)^+(t)\|^2_{L^2(\Omega)},$ we have
$z'(t) - \lambda z(t) \leq 0, z(0)=0.$ Thus $z(t) \leq 0$ for a.e. $t\in [0,T].$
We obtain that $\|(\rho_1-\rho_2)^+(t)\|^2_{L^2(\Omega)} = 0$ for a.e. $t\in [0,T].$ This concludes the proof.
}\fi
\end{proof}

\begin{corollary}
The weak  solutions $\rho$ of \eqref{FP}-\eqref{PBNeumann} are nonnegative whenever the initial condition $\rho_0$ is nonnegative.
\end{corollary}

\begin{proof}
The claim follows from Proposition \ref{Th11.9-Chipot} by setting $\rho_0^1\equiv0{;}$
see also Breiten {\em et al.} \cite{BreitenKunischPfeiffer2018}.
\end{proof}

\section{The optimal control problem} \label{S3}
\setcounter{equation}{0}

In our optimal control problem, we minimize the cost functional
\begin{equation}
\label{cost}
\begin{split}
J(\rho,u) :=&\frac{\alpha_Q}{2} \int_0^T \|\rho(t)-\rho_Q(t)\|^2_{L^2(\Omega)} dt \\
&+ \frac{\alpha_{\Omega}}{2}\|\rho(T)-\rho_\Omega\|_{L^2(\Omega)}^2 + \sum_{i=1}^n \left(\gamma_i\|u_i\|^2_{2} + \beta_i \int_0^T u_i(t) dt \right),
\end{split}
\end{equation}
with {$\rho_Q \in L^2(Q), \rho_\Omega \in L^2(\Omega),$ and} $\beta_i,\gamma_i\geq 0$ for $i=1,\dots,n,$
subject to  the control constraints
\begin{equation}
\label{uconstraint}
{u^{\min}}(t) \leq u(t) \leq u^{\max}(t)\qquad \text{for a.e. } t\in [0,T],
\end{equation}
where the inequalities are defined componentwise, and the bounds $u^{\min},u^{\max}$ belong to $L^\infty(0,T;\mathbb{R}^n).$ The set of admissible controls is
\begin{equation} \label{uad}
\mathcal{U}_{\rm ad} : = \{u\in L^\infty(0,T;\mathbb{R}^n): \, \eqref{uconstraint}\, \text{holds }\}.
\end{equation}
\begin{definition}\label{DefG}
Let us define the {\em control-to-state mapping}
\[
G\colon L^\infty(0,T;\mathbb{R}^n) \to W(0,T),
\]
that associates to each $u \in L^\infty(0,T;\mathbb{R}^n) $ the unique weak solution $\rho \in W(0,T)$ of \eqref{FP}-\eqref{PBNeumann}.
When necessary, we may write $G(u)$ to denote the state $\rho$ corresponding to $u$.
\end{definition}
Our optimal control problem can be rewritten as
\begin{equation}
\label{P}\tag{P}
\min_{u \in \mathcal{U}_{\rm ad}} \,\, J(G(u),u).
\end{equation}

We study two types of solutions, that we define next. We say that $\bar{u}\in \mathcal{U}_{\rm ad}$ is an {\em $L^\infty$-local solution}  (resp., {\em $L^2$-local solution}) of \eqref{P} if there exists $\varepsilon>0$ such that $J(G(\bar{u}),\bar{u})\leq J(G(u),u)$ holds for every $u \in \mathcal{U}_{\rm ad} \cap B^\infty_\varepsilon(\bar{u})$ (resp., $u \in \mathcal{U}_{\rm ad} \cap B^2_\varepsilon(\bar{u})$), where $B^\infty_\varepsilon(\bar{u})$ denotes the open ball in $L^\infty(0,T;\cR^n)$ (resp., in $L^2(0,T;\cR^n)$) of radius $\varepsilon$ centered at $\bar u$.

\subsection{Existence of optimal controls}

\begin{theorem}[Existence of an optimal control]
\label{ThmExist}
There exists (at least) one optimal control for \eqref{P}.
\end{theorem}
\begin{proof}
In view of the control constraints \eqref{uconstraint}, the set of admissible controls $\mathcal{U}_{\rm ad}$ is bounded in $L^\infty(0,T;\mathbb{R}^n).$  Hence, thanks to the estimates \eqref{estV}-\eqref{estrhodot}, the set of states associated to admissible controls is bounded in $W(0,T).$ Therefore, the cost functional $J$ is bounded from below on the set of admissible state-control pairs. Thus, there exists a minimizing sequence $\{(\rho_k,u_k)\} \subset W(0,T) \times \mathcal{U}_{\rm ad},$ where $\rho_k:=G(u_k)$, such that
\[
J(\rho_k,u_k ) \longrightarrow \inf_{u\in \mathcal{U}_{\rm ad}} J(G(u) , u).
\]
Since $\{u_k\}$ is bounded in $L^\infty(0,T;\mathbb{R}^n),$ it contains a weakly$^*$  converging subsequence, thus, keeping the same index, we have
\[
u_k \overset{*}{\rightharpoonup} \bar{u}\quad \text{in } L^\infty(0,T;\mathbb{R}^n)
\]
for some $\bar{u} \in L^\infty(0,T;\mathbb{R}^n).$
The corresponding sequence of states $\{\rho_k\}$ forms a bounded sequence in $W(0,T).$ Thus, there exists $\bar\rho\in W(0,T)$ such that (extracting if necessary a subsequence)
\begin{equation}
\label{limrhok}
\rho_k \rightharpoonup \bar\rho\quad \text{in } W(0,T).
\end{equation}

The main work of the proof is to show that $\bar{\rho}$ is the state associated to $\bar{u},$ and then $J(\bar{\rho},\bar{u}) = \inf_{u\in \mathcal{U}_{\rm ad}} J(G(u) , u).$ For this, we prove that we can pass to the limit in \eqref{varfor}.

Using Aubin-Lions' Lemma \cite{Aubin1963}, we can deduce from \eqref{limrhok} that $\{\rho_k\}$ has a subsequence converging  to $\bar{\rho}$ strongly in $L^2(Q)$. This is, keeping the same index for the subsequence,
\[
\rho_k \longrightarrow \bar\rho\quad \text{(strongly) in } L^2(Q).
\]
Putting all together, we have that $\rho_k \rightharpoonup \bar{\rho}$ weakly in $L^2(0,T;H^1(\Omega))$,  $\displaystyle \partial_t \rho_k \rightharpoonup \partial_t \bar{\rho}$ weakly in $L^2(0,T;H^1(\Omega)^*),$ and $\rho_k\to\bar{\rho}$ strongly in $L^2(Q)$.

We next show that $\bar\rho = G(\bar{u}).$

Given $\varphi \in W^{1,1}_2(Q),$ with $\varphi(\cdot,T)=0,$ we have
\begin{multline*}
\int_0^T\left(\left\langle\partial_t \rho_k (t), \varphi(t) \right\rangle + \int_\Omega \nu\nabla\rho_k(t)\cdot \nabla\varphi(t)dx\right)dt \\
\to
\int_0^T\left( \left\langle \partial_t \bar\rho (t), \varphi(t) \right\rangle + \int_\Omega \nu\nabla\bar{\rho}(t)\cdot \nabla\varphi(t)dx \right)dt.
\end{multline*}
In the weak formulation of the state equation, it only remains to check the convergence of the part containing $B[u_k].$ To this aim,
for $i=1,\dots,n,$ we consider the terms
\begin{equation}
\label{termurho}
\rho_k  B_i[u_k] \frac{\partial \varphi}{\partial x_i}
= \rho_k \big(c_i +  b_i u_{k,i} \big) \frac{\partial \varphi}{\partial x_i}.
\end{equation}
It is easy to see that
$\displaystyle \rho_k c_i \frac{\partial \varphi}{\partial x_i} \to \bar{\rho} c_i \frac{\partial \varphi}{\partial x_i}$ and  $\displaystyle \rho_k b_i \frac{\partial \varphi}{\partial x_i} \to \bar{\rho} b_i \frac{\partial \varphi}{\partial x_i}$ strongly in $L^1(Q)$,
since $\rho_k \to \bar{\rho}$ strongly in $L^2(Q)$,  $c_i,b_i\in L^\infty(\Omega)$ and $\displaystyle\frac{\partial \varphi}{\partial x_i} \in L^2(Q).$
Therefore, in view of the weak$^*$ convergence of $\{u_k\}$ to $\bar{u}$, we get that
\[
\iint_Q u_k  \rho_k b_i \frac{\partial \varphi}{\partial x_i} dx\, dt \to
\iint_Q \bar{u} \bar{\rho} b_i \frac{\partial \varphi}{\partial x_i} dx \, dt, \quad \text{for } k \to \infty.
\]
 Thus, $\bar{\rho}$ is the state associated to $\bar{u}.$  This concludes the proof.
\end{proof}

\section{Properties of the control-to-state mapping} \label{S4}
\setcounter{equation}{0}
In this section, we prove Lipschitz continuity and Fr\'echet differentiability of the control-to-state mapping $G$  and the first order necessary optimality condition of Proposition \ref{PropFirst} as a corollary.

\subsection{Lipschitz continuity of the control-to-state mapping}
We start by proving that $G$ is Lipschitz continuous. This result is then applied in Lemma \ref{LemDeltap} to obtain the analogous property for the associated adjoint state, and both properties are used  later on in the discussion of the second order sufficient optimality condition of Theorem \ref{SOSC}.

\if{
For the proof of the following proposition and for other occasions throughout the remainder of the article, we will use ${\rm div} (\rho b\otimes v)$ to denote the element of $L^2(0,T;H^1(\Omega)^*)$, induced by any pair of functions $(\rho,v)$ with $\rho \in W(0,T)$ and $v\in L^2(0,T;\mathbb{R}^n)$, and defined by
\begin{equation}\label{d}
{\rm div} (\rho b\otimes v)(\varphi):=\int_\Gamma \rho(b\otimes v)\cdot [\bf n}ds- \int_\Omega \rho(t) b \otimes v(t)\cdot \nabla\varphi dx,\quad \text{for }\varphi\in H^1(\Omega),\, \textcolor{blue}{t \in [0,T]}.
\end{equation}
}\fi

For the proof of the next proposition and for other occasions throughout the remainder of the article, we will use the following type of functional: for  given $w \in L^2(\Omega;\mathbb{R}^n)$, we introduce $d[w] \in H^1(\Omega)^*$ defined by
\[
d[w](\varphi) := - \int_\Omega w \cdot \nabla \varphi\, dx
\]
for all $\varphi\in H^1(\Omega).$

\begin{proposition}
\label{Propdeltarho}
For any pair $u_1,u_2 \in L^\infty(0,T;\mathbb{R}^n)$ with associated states $\rho_1:=G(u_1)$, $\rho_2:=G(u_2)$, one has
\begin{equation}
\label{estdeltarho}
\|\rho_2-\rho_1\|^2_{C([0,T];L^2(\Omega))} + \|\rho_2-\rho_1 \|^2_{W(0,T)} \leq C\|\rho_1\|_{C([0,T];L^2(\Omega))}^2\|u_2-u_1\|_2^2,
\end{equation}
where  $C$ depends continuously on $\|u_2\|_\infty$.
\end{proposition}
\begin{proof}
Let us consider two weak solutions $\rho_1, \rho_2$ of \eqref{FP}-\eqref{PBNeumann}, associated to $u_1,u_2$, respectively.
Setting $\delta\rho:=\rho_2-\rho_1,$  we have, for any $\varphi\in H^1(\Omega)$,
\[
\frac{d}{dt} \int_\Omega (\delta\rho) \varphi dx + \int_\Omega\big( \nabla(\delta\rho)+(\delta\rho)B[u_2]\big)\cdot \nabla\varphi dx ={-\int_\Omega \rho_1  b\otimes (u_2-u_1)}\cdot \nabla\varphi dx.
\]
Applying the estimates \eqref{estH}-\eqref{estrhodot} of Theorem \ref{ThmEstimates}, we get
\begin{equation}
\label{rho2rho1}
\|\delta\rho\|^2_{C ([0,T];L^2(\Omega))} + \|\delta\rho\|^2_{W(0,T)} \leq  \frac{C_{u_2}}{N} \|d[\rho_1b\otimes (u_2-u_1)]\|_{L^2(0,T;H^1(\Omega)^*)}^2,
\end{equation}
where $C_{u_2}$ depends on $\|u_2\|_\infty.$
We have, for a.a. $t\in [0,T],$
\begin{equation*}
\label{fest}
\begin{split}
\left\|d\big[\rho_1(t) \right.&\left.b\otimes(u_2(t)-u_1(t))\big]\right\|_{H^1(\Omega)^*}\\
&= \sup_{\varphi \in H^1(\Omega),\|\varphi\|_{H^1(\Omega)}=1} \left| \int_\Omega \rho_1(t)\,b \otimes (u_2-u_1)(t) \cdot\nabla\varphi \right| dx \\
&\leq \|\rho_1(t)  \,b \otimes (u_2-u_1)(t)\|_{L^2(\Omega)^n}.
\end{split}
\end{equation*}
Thus,
\begin{equation}
\label{estdrho1}
\begin{split}
\|d[\rho_1b\otimes &(u_2-u_1)]\|_{L^2(0,T;H^1(\Omega)^*)}^2 \leq \| \rho_1\,b \otimes (u_2-u_1) \|_{L^2(0,T;L^2(\Omega)^n)}^2 \\
&=  \iint_Q \Big[\rho_1(x,t) \sum_i  b_i(x)(u_{2,i}(t)-u_{1,i}(t))\Big]^2 dx dt\\
&\leq \|b\|_{L^\infty(\Omega)^n}^2 \int_0^T \Big[ \sum_i |u_{2,i}(t)-u_{1,i}(t)|\Big]^2 \int_\Omega \rho_1 (x,t)^2 dx dt\\
&\leq n \| b\|_{L^\infty(\Omega)^n}^2 \|\rho_1 \|_{L^\infty(0,T;L^2(\Omega))}^2 \|u_2-u_1\|_2^2.
\end{split}
\end{equation}
Hence, in view of \eqref{rho2rho1} and \eqref{estdrho1}, we get the desired estimate \eqref{estdeltarho} for $C:= \frac{nC_{u_2}}{N} \| b\|_{L^\infty(\Omega)^n}^2.$ This concludes the proof.
\end{proof}

\subsection{Fr\'echet differentiability of the control-to-state mapping}
Our next  step is to prove the differentiability of $G$ by application of the Implicit Function Theorem.
To this aim, we define the mapping
\begin{equation}
\label{calG}
(\rho,u) \longmapsto  \mathcal{G}(\rho,u):= \left(\partial_t \rho+a[u](\rho,\cdot),\rho(0)-\rho_0\right)
\end{equation}
from $W(0,T)\times L^\infty(0,T;\cR^n)$ to $L^2(0,T;H^1(\Omega)^*)\times L^2(\Omega).$ The state equation can be viewed as the equation
\begin{equation}
\mathcal{G}(\rho,u) =0.
\end{equation}

\begin{proposition}
\label{GCinfty}
The mapping $\mathcal{G}$ is of class $C^\infty.$
\end{proposition}

\begin{proof}
The first component $\mathcal{G}_1$ of $\mathcal{G}$ has the form
\[
\mathcal{G}_1(\rho,u)(\varphi) = \left\langle \partial_t \rho , \varphi\right\rangle + \int_\Omega \nabla \rho \cdot \nabla \varphi \, dxdt + \int_\Omega \rho B[u] \cdot \nabla \varphi \, dxdt.
\]
Its first two summands clearly define linear and continuous mappings from $W(0,T)$ to $L^2(0,T;H^1(\Omega)^*)$, hence they are of class $C^\infty$. Therefore, it suffices to confirm that the operator
\[
(\rho,u) \mapsto \rho B[u]
\]
is of class $C^\infty$ from  $W(0,T)\times L^\infty(0,T;\cR^n)$ to $L^2(Q)$. This, however, follows easily from the quadratic nature of $\rho B[u]$. For increments $\sigma \in W(0,T)$ and $v \in L^\infty(0,T;\cR^n)$, we have
\[
\begin{aligned}
(\rho + \sigma)B[u + v] &= (\rho + \sigma)\big( c + b \otimes(u+v)\big)\\
& = \rho B[u] + \rho b \otimes v + \sigma(c + b \otimes u) + \sigma b \otimes v.
\end{aligned}
\]
This is a second order Taylor expansion of the mapping $(\rho,u) \mapsto \rho B[u]$ with continuous linear and quadratic parts, and vanishing remainder term. Let us exemplarily confirm the continuity of the quadratic term. We have
\[
\|\sigma b \otimes v\|_{L^2(Q)^n} \le \| b \|_{L^\infty(\Omega)^n} \|\sigma\|_{C([0,T];L^2(\Omega))}\|v\|_2
\le c \|\sigma\|_{W(0,T)}\|v\|_2,
\]
hence the continuity of the quadratic form $\sigma b \otimes v$. Therefore $\mathcal{G}_1$ is of class
$C^\infty$. The second component of $\mathcal{G}$ is obviously of class $C^\infty$.
\end{proof}

\begin{corollary}
\label{CorGCinfty}
The control-to-state mapping $G$ is of class $C^\infty.$
\end{corollary}

\begin{proof}
In view of Proposition \ref{GCinfty}, we have that $\mathcal{G}$ is of class $C^\infty$ and that
\[
\partial_\rho \mathcal{G}(\rho,u)z = \big( \partial_t z+a[u](z,\cdot),z(0)\big).
\]
Given $f\in L^2(0,T;H^1(\Omega)^*)$ and $z_0\in L^2(\Omega),$ the system
\begin{equation}
\begin{split}
\partial_t z+a[u](z,\cdot) &= f,\\
z(0) &= z_0,
\end{split}
\end{equation}
is again a Fokker-Planck equation as the state equation and, therefore, it has a unique (weak) solution $z[z_0,f]$  that belongs to $W(0,T),$ and depends continuously on $z_0 \in L^2(\Omega)$ and on $f \in L^2(0,T;H^1(\Omega)^*).$
Therefore, thanks to the existence and uniqueness Theorem \ref{ThmExistSol}, $\displaystyle \partial_\rho \mathcal{G}(\rho,u)$ is an isomorphism from
$W(0,T)$ to $L^2(0,T;H^1(\Omega)^*) \times L^2(\Omega)$.
Thus, the hypotheses of the Implicit Function Theorem are satisfied, and then $\mathcal{G}(\rho,u)=0$ implicitly defines the control-to-state operator $G: u \mapsto \rho$ that is itself of class $C^\infty.$
\end{proof}

\begin{proposition}
\label{PropG'}
Let $u\in L^\infty(0,T;\cR^n)$ and $\rho:=G(u).$
For any $v\in L^\infty(0,T;\mathbb{R}^n),$ we have that $G'(u)v=z,$ where $z$ belongs to $W(0,T)$ and is the weak solution of the {\em linearized state equation} at $(\rho,u)$ that is given by
\begin{equation}
\label{lineq}
\begin{split}
\partial_t z+a[u](z,\cdot) &= d[\rho b\otimes v],\\
z(0) &= 0.
\end{split}
\end{equation}
\end{proposition}

\begin{proof}
{The fact that $z$ is in $W(0,T)$ follows from Theorem \ref{ThmExistSol} as already observed in the proof of Corollary \ref{CorGCinfty} above.}
The representation $G'(u)v=z$ follows from a direct application of the Implicit Function Theorem by differentiating
the state equation \eqref{FP}-\eqref{PBNeumann} with respect to the control.\end{proof}

The linearized state equation \eqref{lineq} at $(\rho,u)$ in the direction of $v,$ is to be understood as
\begin{equation}
\label{lineareq2}
\begin{split}
\int_\Omega \big( \partial_t z\varphi + (\nu\nabla z+zB[u])\cdot \nabla \varphi \big) dx &= - \int_\Omega \rho ( b\otimes v)\cdot \nabla \varphi dx \,\, \text{ for all } \varphi \in H^1(\Omega),\\
z(0) &= 0,
\end{split}
\end{equation}
or, in the strong form:
\begin{equation}
\label{lineareq3}
\begin{split}
\partial_t z-\nu\Delta z-{\rm div} (zB[u]) &= {\rm div} \big( \rho (b \otimes v) \big),\quad \text{in } Q,\\
\big(\nu\nabla z+zB[u]\big)\cdot n &= - \big(\rho  b\otimes v\big)\cdot n,\quad \text{in } \Sigma,\\
z(0) &= 0\quad \text{in } \Omega.
\end{split}
\end{equation}
In fact, integration by parts in the last two terms in the l.h.s. and in the term in the r.h.s of \eqref{lineareq2} yields, respectively,
\begin{gather*}
\int_\Omega \nu\nabla z \cdot \nabla \varphi dx = \nu \int_\Gamma \varphi\nabla z \cdot n ds - \nu \int_\Omega\varphi\Delta z dx,\\
\int_\Omega zB[u] \cdot \nabla \varphi dx = \int_\Gamma \varphi z B[u]\cdot nds - \int_\Omega\varphi {\rm div} (zB[u]) dx,\\
 -\int_\Omega \rho ( b\otimes v)\cdot \nabla \varphi dx = -\int_\Gamma\varphi\rho ( b\otimes v)\cdot n ds + \int_\Omega  \varphi {\rm div}(\rho( b\otimes v)) dx ,
\end{gather*}
from which \eqref{lineareq3} follows.

\begin{remark}
Following Ladyzhenskaya et al. \cite{Ladyzhenskaya_etal1968},  an equivalent variational formulation of the linearized state equation \eqref{lineareq2} is given as follows: $z$ is a function in $W_2^{1,0}(\Omega)$ such that
\begin{multline}
\label{lineareqLady}
\iint_Q \left( - z \partial_t \varphi + (\nu\nabla z+zB[u])\cdot \nabla \varphi \right) dx\, dt = - \iint_Q \rho (b\otimes v) \cdot \nabla \varphi dx\,dt\\
\text{for all }\varphi \in W^{1,1}_2(Q) \text{ with } \varphi(\cdot,T)=0.
\end{multline}
\end{remark}

Let us define the {\em reduced cost functional} as
\[
F(u):=J(G(u),u).
\]
By the chain rule, the reduced cost functional $F$ is continuously Fr\'echet differentiable, since $J$ and $G$ have this property.

\begin{proposition}[First order necessary condition]
\label{PropFirst}
If $\bar{u}$ is an $L^\infty$-local minimum for \eqref{P}, then
\begin{equation}
\label{FirstAbstract}
F'(\bar{u})(u-\bar{u}) \geq 0\quad \text{for every } u\in \mathcal{U}_{\rm ad}.
\end{equation}
\end{proposition}

\begin{proof}
Given that $\bar{u}$ is an $L^\infty$-local minimum, there exists $\bar\varepsilon\in (0,1)$ such that, if $0<\varepsilon<\bar\varepsilon,$ then
\[
F\big(\bar{u}+\varepsilon(u-\bar{u})\big)-F(\bar{u}) \geq 0\quad \text{for every } u \in \mathcal{U}_{\rm ad}.
\]
Dividing by $\varepsilon$ and passing to the limit yields the desired standard result.
\end{proof}

\section{The adjoint equation} \label{S5}
\setcounter{equation}{0}

\subsection{Definition of the adjoint equation}

By an adjoint state, the variational inequality \eqref{FirstAbstract} can be transformed to a more convenient form. To this aim, we introduce the following adjoint equation
for the adjoint state $p$ associated with $(\rho,u)$:
\begin{align}
\label{adjoint}
- \partial_t p - \nu \Delta p + B[u]\cdot \nabla p &= {\alpha_Q}(\rho-\rho_Q)\quad \text{in } Q,\\
 p(T) &={\alpha_\Omega}( \rho(T)-\rho_\Omega)\quad \text{ in } \Omega,\\
 \partial_np &= 0\quad \text{on } \Sigma.
\end{align}
The form of the adjoint equation can be found, e.g. by application of a formal Lagrangian technique,
cf. \cite{Troltzsch2010}[chpt. 2.6].
 In weak formulation, the adjoint equation at $(\rho,u)$ is defined by
\begin{equation}
\label{adjointeq}
\begin{split}
 \int_\Omega  \Big(- p_t\varphi  + \nu \nabla p \cdot \nabla \varphi  + {\varphi }\, B[u]\cdot \nabla p \Big)dx &= {\alpha_Q} \int_\Omega (\rho-\rho_Q)\varphi dx\\
 &\quad \text{on } \mathcal{D}'(0,T)\,\,\text{for all } \varphi\in H^1(\Omega),\\
  p(T) &= {\alpha_\Omega} (\rho(T)-\rho_\Omega).
  \end{split}
 \end{equation}
  Its unique (weak) solution
$p$ is called the {\em adjoint state} associated with $(\rho,u)$.
Note that, with  $f:= {\alpha_Q}(\rho-\rho_Q),$ \eqref{adjointeq} can be rewritten as
 \begin{equation}
 \label{adjointeqgen}
\begin{split}
  -\partial_t p + a[u(\cdot)](\cdot,p) &= f\quad \text{in } L^2(0,T;H^1(\Omega)^*),\\
  p(T) &= {\alpha_\Omega}(\rho(T)-\rho_\Omega).
  \end{split}
 \end{equation}

 %

 \begin{proposition} \label{P5.1}
 Given any $f\in L^2(0,T;H^1(\Omega)^*)$,  equation \eqref{adjointeqgen} has a unique weak solution $p$, it belongs to $W(0,T)$ and the following estimate holds,
 \begin{equation}
\|p\|_{W(0,T)}+\|p\|_{C([0,T];L^2(\Omega))} \leq C \big( \| p(T)\|_{L^2(\Omega)} + \|f\|_{L^2(0,T;H^1(\Omega)^*)}  \big).
 \end{equation}
 In particular, for the adjoint equation \eqref{adjointeq}, one has
  \begin{equation}
  \label{pest}
\|p\|_{W(0,T)}+\|p\|_{C([0,T];L^2(\Omega))}
 \leq C \big(\|\alpha_\Omega  (\rho(T) - \rho_\Omega) \|_{L^2(\Omega)} +
\| \alpha_Q(\rho-\rho_Q) \|_{L^2(Q)}  \big).
 \end{equation}
\end{proposition}
 %
 \begin{proof}
This follows from Theorem \ref{ThmExistSol} by the transformation of time $\tau = T - t$ and $\tilde p(\tau) = p(T-\tau)$. In particular, the terminal condition for $p$ becomes an initial condition for $\tilde p$.
\end{proof}

\subsection{First order necessary optimality conditions in terms of the adjoint state}
In this subsection we rewrite the first order condition of Proposition \ref{PropFirst} in terms of the adjoint state. For this, we show the following technical result.

\begin{lemma}
Assume that  $u \in L^\infty(0,T;\mathbb{R}^n)$ is given and let $\rho=G(u)$ be its associated state.
Let $z$ be the weak solution of the linearized state equation \eqref{lineq} corresponding to $v\in L^\infty(0,T;\mathbb{R}^n),$ and $p$ the weak solution of the adjoint equation \eqref{adjoint}. Then
\begin{equation}
\label{zp}
{\alpha_Q}\int_0^T\int_\Omega(\rho-\rho_Q)z dx dt + {\alpha_\Omega}\int_\Omega(\rho(T)-\rho_\Omega)z(T)dx = -\int_0^T \int_\Omega\rho( b\otimes v)\cdot \nabla p dx dt.
\end{equation}
\end{lemma}

\if{
or, equivalently,
\begin{multline}
\int_0^T\int_\Omega(\rho-\rho_Q)z dx dt + \int_\Omega(\rho(T)-\rho_\Omega)z(T)dx \\= \int_0^T \int_\Omega p\, {\rm div}(\bar{\rho}( b\otimes v)) dt - \int_0^T\int_{\Gamma} p\bar{\rho}( b\otimes v) \cdot n ds dt\\
=\int_0^T \int_\Omega p\, {\rm div}(\bar{\rho}( b\otimes v)) dt + \int_0^T\int_{\Gamma} p(\nu\nabla z+z B[\bar{u}]) \cdot n ds dt,
\end{multline}
where the second equality follows from the boundary condition of $z.$
}\fi

\begin{proof}
We write down the weak formulation of the linearized equation \eqref{lineareq2} with test function $\varphi := p$, where $p$ is the adjoint state associated with $u$, i.e. the solution of the adjoint equation \eqref{adjointeq}. Moreover, we insert in the adjoint equation \eqref{adjointeq} the test function $\varphi := z$. This yields
 \begin{gather}
 \label{ztp}
  \int_0^T \left\langle \partial_t z , p \right\rangle\, dt + \int_0^T \int_\Omega \Big(\nu\nabla z\cdot \nabla p +zB[u]\cdot \nabla p\Big) {dxdt}= -  \int_0^T \int_\Omega \rho ( b\otimes v)\cdot \nabla p dxdt\\
\label{ptz}
\int_0^T \left\langle - \partial_t p , z \right\rangle\, dt + \int_0^T \int_\Omega \Big(\nu \nabla p \cdot \nabla z  + zB[u]\cdot \nabla p \Big)dxdt = {\alpha_Q}\int_0^T \int_\Omega (\rho-\rho_Q)\,zdxdt.
\end{gather}
Performing an integration by parts in \eqref{ptz},  invoking the conditions $z(0) = 0$, and
$p(T) = {\alpha_\Omega}(\rho(T)-\rho_\Omega)$, equation \eqref{ptz} is transformed to
\[
\begin{aligned}
& \int_0^T\left\langle p,\partial_t z \right\rangle dt+ \int_0^T \int_\Omega \nu \nabla p \cdot \nabla z  + zB[u]\cdot \nabla p\Big)dxdt \\
&\qquad \qquad =  {\alpha_\Omega}  \int_\Omega (\rho(T)-\rho_\Omega) z(T)\,dx + {\alpha_Q} \int_0^T \int_\Omega (\rho-\rho_Q)\,zdxdt.
\end{aligned}
\]
The left-hand side of latter equation is equal to the left-hand side of \eqref{ztp}. Therefore,
also the associated right-hand sides must be equal. This, however, confirms the claim of the lemma.
\end{proof}

Let us introduce the notation
\begin{equation}
\label{Phi}
\Phi_i(t):=-\int_\Omega \rho(t)  b_i\frac{\partial {p}(t)}{\partial x_i} dx + \gamma_i u_i(t)+ \beta_i \quad \text{for } i=1,\dots,n.
\end{equation}
\begin{theorem}
For any $u,v\in L^\infty(0,T;\mathbb{R}^n),$ one has
\begin{equation}
\label{F'exp}
F'(u)v=\sum_{i=1}^n \int_0^T \Phi_i(t) v_i(t)dt,
\end{equation}
where $\Phi$ is defined in \eqref{Phi}, with $\rho = G(u)$ and $p$ being the associated adjoint state.
\end{theorem}

\begin{proof}
Let $z$ be the solution of the state equation linearized at $(\rho,u).$
Then one has
\begin{equation}\label{F'}
\begin{split}
&F'(u)v= D_{\rho}J(G(u),u)G'(u)v + D_uJ(G(u),u)v \\
& = {\alpha_Q}  \int_0^T\int_\Omega(\rho-\rho_Q)z dx dt + {\alpha_\Omega}  \int_\Omega(\rho(T)-\rho_\Omega)z(T)dx +\sum_{i=1}^n \int_0^T \left(\gamma_iu_iv_i+ \beta_i v_i\right)dt\\
&= -\int_0^T \int_\Omega\rho( b\otimes v)\cdot \nabla p dx dt+\sum_{i=1}^n \int_0^T \left(\gamma_i u_iv_i+ \beta_i v_i\right)dt\\
&=\sum_{i=1}^n \int_0^T \Big[-\int_\Omega \rho  b_i\frac{\partial p}{\partial x_i} dx + \gamma_i u_i+ \beta_i \Big]v_i\,dt,
\end{split}
\end{equation}
where we used \eqref{zp} in the second equality. This proves the result.
\end{proof}

By means of the expression \eqref{F'exp}, we can reformulate the first order optimality condition of Proposition \ref{PropFirst} as follows:
\begin{corollary}
If $\bar{u}$ is an $L^\infty$-local minimum for \eqref{P}, then
\begin{equation}
\label{First}
\sum_{i=1}^n \int_0^T \bar\Phi_i(t)\big(u_i(t)-\bar{u}_i(t)\big)\,dt \geq 0\quad \text{for all } u\in \mathcal{U}_{\rm ad},
\end{equation}
where {$\bar\Phi$ is the function given in \eqref{Phi} associated to $\bar{u}$.
Consequently, we have
\begin{equation}
\label{signPhi}
\left\{
\begin{split}
& \displaystyle \bar{\Phi}_i(t) > 0   \Longrightarrow  \bar{u}_i(t) = u_i^{\min}(t),\\
& \displaystyle \bar{\Phi}_i(t) < 0   \Longrightarrow  \bar{u}_i(t) = u_i^{\max}(t),\\
& \displaystyle u_i^{\min}(t) < \bar{u}_i (t) < u_i^{\max}(t)  \Longrightarrow \bar{\Phi}_i(t) =0,
\end{split}
\right.
\end{equation}
a.e. on $[0,T]$ and for all $i=1,\dots,n.$}
\end{corollary}

\section{Second order analysis}\label{S6}
\setcounter{equation}{0}

To establish second order optimality conditions, we will apply general results by Casas and Tr\"oltzsch,   \cite[Theorems 2.2 and 3.3]{CasasTroeltzsch2012}.
For this purpose, we have to verify that the conditions {\rm \ref{C1}}-{\rm \ref{C3}} below  are satisfied for problem \eqref{P}. More precisely, we will need {\rm \ref{C1}} in the second order necessary condition of Theorem \ref{SONC} below, while {\rm \ref{C1}}-{\rm \ref{C3}} are used in the sufficient one of Theorem \ref{SOSC}.

\subsection{Second order conditions for an optimization problem in Banach spaces}

We consider a Banach space $U_\infty$ and a Hilbert space $U_2,$ endowed with the norms $|\cdot|_\infty$ and $|\cdot|_2$, respectively, and such that $U_\infty$ is continuously embedded in $U_2.$ Let us introduce the abstract optimization problem
\begin{equation}
\label{calP}\tag{$\mathcal{P}$}
\min_{u\in \mathcal{K}} \mathcal{J}(u),
\end{equation}
where $\mathcal{K}\subseteq U_\infty$ is a given nonempty convex set and $\mathcal{J} \colon \mathcal{A}\to \mathbb{R}$ is the objective function, defined  and twice continuously differentiable in an open subset $\mathcal{A} \subset U_\infty$ that covers $\mathcal{K}.$
We say that $\bar{u}$ is an {\em $U_\infty$-local solution} of \eqref{calP} if there exists $\varepsilon>0$ such that $\mathcal{J}(\bar{u}) \leq \mathcal{J}(u)$ holds for all $u\in \mathcal{K} \cap \{u\in U_\infty: \|u-\bar{u}\|_\infty<\varepsilon\}.$
\\
\indent If $\bar{u}$ is an $U_\infty$-local solution of \eqref{calP}, then the following first order necessary condition is satisfied:
\begin{equation}
\label{FONCJ}
\mathcal{J}'(\bar{u})(u-\bar{u})\geq 0 \quad \text{for all } u\in \mathcal{K}.
\end{equation}

 Let us fix $\bar{u}$ in $\mathcal{K}$. We consider the following conditions for problem \eqref{calP}. All the notions of differentiability of $\mathcal{J}$ are to be understood in the sense of $U_\infty.$
\begin{itemize}
\item[\namedlabel{C1}{(C1)}] The functional $\mathcal{J}$ is of class $C^2$ in
in $\mathcal{A}$. For every $u \in \mathcal{K},$ there exist continuous extensions
\[
\mathcal{J}'(u) \in \mathcal{L}(U_2;\mathbb{R}),\quad \mathcal{J}''(u)  \in \mathcal{B}(U_2;\mathbb{R})
\]
of $\mathcal{J}'(u)$ and $\mathcal{J}''(u)$, where $\mathcal{B}(U_2;\mathbb{R})$ denotes the Banach space of continuous bilinear real functionals on  $U_2\times U_2$.
\item[\namedlabel{C2}{(C2)}] For any sequence $\{ (u_k,v_k)\} \subset \mathcal{K} \times U_2 $ with $u_k\to \bar{u}$ in $U_2$ and $v_k \rightharpoonup v$ weakly in $U_2$ as $k \to \infty$, there holds
\begin{gather}
\label{limF'}
\mathcal{J}'(\bar{u}) v = \lim_{k\to\infty} \mathcal{J}'(u_k)v_k.
\end{gather}
\item[\namedlabel{C3}{(C3)}] For any sequence defined as in {\rm (C2),} the following two properties are satisfied for some $\Lambda>0:$ it holds
\begin{gather}
\label{liminfD2J} \mathcal{J}''(\bar{u}) v^2 \leq \liminf_{k\to \infty} \mathcal{J}''(u_k) v_k^2,\\
\label{liminfD2J0} \text{and, if } v=0, \text{ then } \,\,\, \Lambda \liminf_{k\to\infty} |v_k|_2^2 \leq \liminf_{k\to \infty} \mathcal{J}''(u_k)v_k^2.
\end{gather}
\end{itemize}

For a fixed control $\bar{u} \in \mathcal{K}$, let us define the following sets
\begin{equation}
\label{cones}
 \begin{split}
S(\bar{u})&:= \big\{v\in U_\infty: v=\lambda(u-\bar{u})\, \text{for some } \lambda>0 \text{ and } u\in \mathcal{K} \big\},\\
C(\bar{u})&:= {\rm cl}_{U_2}(S(\bar{u})) \cap\{ v\in U_2: \mathcal{J}'(\bar{u})v=0\},\\
D(\bar{u}) &:= \big\{v\in S(\bar{u}) : \mathcal{J}'(\bar{u})v=0 \big\}.
\end{split}
\end{equation}
The set $S(\bar{u})$ is called {\em cone of feasible directions}, while $C(\bar{u})$ is the {\em critical cone.}
\\
\vspace{1.5ex}
\\
\indent
We first state Theorems 2.2 and 2.3 from \cite{CasasTroeltzsch2012} that we will apply to obtain the second order necessary and sufficient optimality conditions in Theorems \ref{SONC} and \ref{SOSC} below.

\begin{theorem}[Casas-Tr\"oltzsch \cite{CasasTroeltzsch2012}]
\label{SONC2012}
Let $\bar{u}$ be an $U_\infty$-local solution for \eqref{P}. Assume that {\rm \ref{C1}} and the {\em regularity condition} $C(\bar{u}) = {\rm cl}_{U_2}D(\bar{u})$ hold. Then
\[
\mathcal{J}''(\bar{u})v^2 \geq 0\quad \text{for all } v\in C(\bar{u}).
\]
\end{theorem}

\begin{theorem}[Casas-Tr\"oltzsch \cite{CasasTroeltzsch2012}]
\label{SOSCCT}
Suppose that {\rm \ref{C1}-\ref{C3}} are fulfilled. Let $\bar{u} \in \mathcal{K}$ satisfy  the
first order necessary condition \eqref{FONCJ} along with
\begin{equation}
\mathcal{J}''(\bar{u})v^2 >0\quad \text{for all } v \in C(\bar{u}) \backslash \{0\}.
\end{equation}
Then, there exist $\varepsilon>0$ and $\delta >0$ such that
\begin{equation}
\mathcal{J}(\bar{u}) + \frac{\delta}{2} |u-\bar{u}|_2^2 \leq \mathcal{J}(u)\quad  \text{for all } u \in \mathcal{K} \cap B_{2,\varepsilon}(\bar{u}),
\end{equation}
where $B_{2,\varepsilon}(\bar{u})$ is the open ball in $U_2,$ of radius $\varepsilon$ and centered in the origin.
\end{theorem}
\vspace{1.5ex}
In the remainder of this paper, we will confirm the three conditions {\rm \ref{C1}-\ref{C3}} for our optimal control problem.
For this purpose, we consider
\[
U_\infty=\mathcal{A}:= L^\infty(0,T;\cR^n),\quad U_2:= L^2(0,T;\cR^n), \quad \mathcal{J}:= F,\quad \mathcal{K}:=\mathcal{U}_{\rm ad}.
\]
Some continuity and extension properties of the second derivative of the reduced cost $F$
will be crucial for our next arguments.

\subsection{Second derivative of the reduced cost functional}

For a second order analysis of our optimal control problem, we have to
extend the first- and (later) second-order derivatives with respect to the control to the Hilbert space $L^2(0,T;\mathbb{R}^n).$ To this aim, we first show that that the linearized equation \eqref{lineq} is well-posed when $v\in L^2(0,T;\mathbb{R}^n).$
\begin{proposition}
\label{Propestimatez}
For every $v\in L^2(0,T;\mathbb{R}^n)$, the linearized equation \eqref{lineq}  at $u$ with state $\rho = G(u)$ has a unique solution $z$ in $W(0,T).$ The following estimate holds true,
\begin{equation}
\label{estimatez}
\|z\|_{{C([0,T];L^2(\Omega))}} + \|z\|_{W(0,T)}
\leq K \| b\|_{L^\infty(\Omega)^n} \|\rho\|_{{C([0,T];L^2(\Omega))}} \|v\|_2,
\end{equation}
where $K$ depends continuously on $\|u\|_\infty$ but does not depend on $v$.
\end{proposition}

\begin{proof}
Existence and uniqueness of a solution to \eqref{lineq} for $v\in L^2(0,T;\mathbb{R}^n)$ follow straightforwardly from Theorem \ref{ThmExistSol}. It remains to estimate the $L^2(0,T;H^1(\Omega)^*)$-norm of the r.h.s. function $d[\rho b\otimes v]$, in order to apply the estimates of Theorem \ref{ThmEstimates}. Take any $\varphi \in H^1(\Omega);$ then for a.a. $t \in [0,T]$ we have
\begin{equation*}
\begin{split}
\big| d[\rho(t) b\otimes v(t)](\varphi)\big| &=
\left| \int_\Omega \rho(t) ( b\otimes v(t)) \cdot \nabla\varphi dx \right|\\
&\leq \sum_{i=1}^n \left| v_i(t) \int_\Omega \rho(x,t) b_i(x)\frac{\partial\varphi}{\partial x_i}(x) dx\right|\\
&\leq \sum_{i=1}^n \left| v_i(t)\right|  \| b_i\|_{L^\infty(\Omega)} \|\rho(t)\|_{L^2(\Omega)}  \|\varphi\|_{H^1(\Omega)}.
\end{split}
\end{equation*}
Thus, for a.e. $t\in [0,T],$
$$
\left\| d[\rho(t) b\otimes v(t)] \right\|_{H^1(\Omega)^*}\leq \| b\|_{L^\infty(\Omega)^n} \|\rho(t)\|_{L^2(\Omega)} \sum_{i=1}^n \left| v_i(t)\right|,
$$
and hence
$$
\left\|d[\rho b\otimes v]\right\|_{L^2(0,T;H^1(\Omega)^*)}\leq M \| b\|_{L^\infty(\Omega)^n} \|\rho\|_{L^2(0,T;L^2(\Omega))} \|v\|_2,
$$
for some constant $M$ depeding on the control dimension $n$.
The estimate \eqref{estimatez} follows from Theorem \ref{ThmEstimates} with $f:= d[\rho b\otimes v]$. This concludes the proof.
\end{proof}

Let us compute the second derivative $G''(u)[v_1,v_2]$ of the control-to-state operator. Notice that the differentiability of $G$ has already been proven. We have
\[
\partial_t G(u) - \nu \Delta G(u) - {\rm div} (B[u] G(u)) =0
\]
subject to associated initial and boundary conditions.
Then, differentiating with respect to $v_1 \in L^\infty(0,T;\mathbb{R}^n)$, we obtain
\[
\partial_t (G'(u) v_1) -\nu \Delta (G'(u)v_1)-{\rm div}( b\otimes v_1\, G(u)) - {\rm div} (B[u] G'(u)v_1)=0.
\]
Another differentiation with respect to $v_2 \in L^\infty(0,T;\mathbb{R}^n)$ yields
\begin{multline*}
\partial_t (G''(u)[v_1,v_2])- \nu\Delta(G''(u)[v_1,v_2])-{\rm div}( b \otimes v_1 G'(u)v_2)\\
-{\rm div}( b\otimes v_2 G'(u)v_1) - {\rm div}(B[u] G''(u)[v_1,v_2])=0.
\end{multline*}
Setting $z_i := G'(u)v_i,$ for $i=1,2,$ and $w:=G''(u)[v_1,v_2]$ for the unknown second-order derivative, we get the following equation for $w$, written in strong form:
\begin{equation*}
\begin{split}
w_t- \nu\Delta w  - {\rm div}(B[u] w) &= {\rm div} \, (z_1\, b \otimes v_2 + z_2\, b \otimes v_1),\\
(\nu\nabla w + wB[u])\cdot n &= -(z_1\, b \otimes v_2 + z_2\, b \otimes v_1) \cdot n ,\\
w(0)&=0.
\end{split}
\end{equation*}
The boundary condition is also obtained by implicit differentiation.
In weak formulation, the equation for $w = G''(u)[v_1,v_2]$ reads
\begin{equation*}
\begin{split}
\int_0^T \Big\langle \partial_t w , \varphi \Big\rangle\, dt &+ \iint_Q ( \nu\nabla w + w B[u]) \cdot \nabla \varphi \, dxdt\\
&=
- \iint_Q (z_1 \, b \otimes v_2 + z_2\, b \otimes v_1) \cdot \nabla \varphi \, dxdt, \quad \forall \varphi \in L^2(0,T;H^1(\Omega)),\\
w(0)&=0.
\end{split}
\end{equation*}
Notice that $v_1$ and $v_2$ depend only on $t$. Moreover,  we have $b \in L^{\infty}(\Omega;\mathbb{R}^n)$ and $z_i \in W(0,T)$, for
$i = 1,2$. In particular, $z_i$ belongs to $C([0,T];L^2(\Omega))$, for
$i = 1,2$.
Therefore, $z_1\, b \otimes v_2 + z_2\, b \otimes v_1$ are well defined
in $L^2(Q)$, since the $z_i$ are bounded w.r.t. $t$. The homogeneous initial condition follows obviously by differentiating the equation
$G(u)(\cdot,0) = \rho_0(\cdot)$ twice with respect to $u$.
The associated weak formulation is
\begin{align}
\label{varforw}
\frac{d}{dt} ( w, \varphi) + a[u](w,\varphi) &=\big(d[z_1 b\otimes v_2] + d[z_2 b\otimes v_1]\big)(\varphi)\\
&\qquad\quad \text{on } \mathcal{D}'(0,T)\,\, \text{for all } \varphi \in H^1(\Omega),\\
\label{initcondw} w(0) &= 0\quad \text{in } L^2(\Omega).
\end{align}

\begin{proposition}
Let $v_1,v_2 \in L^\infty(0,T;\cR^n)$ be control increments with associated linearized states $z_1,z_2$, respectively.
Then the second order derivative of the reduced cost functional $F$ at $u$ in the direction pair $(v_1,v_2)$  is given by
\begin{multline}
F''(u)[v_1,v_2] = \int_0^T\int_\Omega \Big({\alpha_Q} z_1z_2   - \nabla p \cdot ( z_2 b \otimes v_1 +z_1 b\otimes v_2 )\Big)dxdt \label{F''}
\\
+ \int_0^T\sum_{i=1}^n \gamma_i v_{1,i} v_{2,i} dt + {\alpha_{\Omega}} \int_\Omega z_1(T)z_2(T) dx,
\end{multline}
where $p$ is the adjoint state associated with $u$.
\end{proposition}

\begin{proof}
The second order derivative of $F$ is computed by the chain rule. One has
\begin{multline*}
F'(u)v_1={\alpha_{Q}}\big(G(u)-\rho_Q,G'(u)v_1\big)_{L^2(Q)}+{\alpha_{\Omega}}\big((G(u))(T)-\rho_\Omega,(G'(u)v_1)(T)\big)_{L^2(\Omega)}\\
+\sum_{i=1}^n \Big( \gamma_i(u_i,v_{1,i})_{L^2(0,T)} + \beta_i\int_0^T v_{1,i} dt \Big).
\end{multline*}
Then, differentiating this expression w.r.t. $v_2$,
\begin{equation}
\label{F''exp}
\begin{split}
&F''(u)[v_1,v_2] =  {\alpha_Q} \big( G'(u)v_2,G'(u)v_1\big)_{L^2(Q)} + {\alpha_Q} \big( G(u)-\rho_Q,G''(u)[v_1,v_2]\big)_{L^2(Q)} \\
& + {\alpha_\Omega} \big((G'(u)v_2)(T),(G'(u)v_1)(T)\big)_{L^2(\Omega)} \\
&+ {\alpha_\Omega} \big((G(u))(T)-\rho_\Omega,(G''(u)[v_1,v_2])(T)\big)_{L^2(\Omega)}+\sum_{i=1}^n \gamma_i(v_{2,i},v_{1,i})_{L^2(0,T)}\\
&= {\alpha_Q}(z_1,z_2)_{L^2(Q)} + {\alpha_Q}(G(u)-\rho_Q,w)_{L^2(Q)} + {\alpha_\Omega}(z_1(T),z_2(T))_{L^2(\Omega)}\\
&+{\alpha_\Omega}(p(T),w(T))_{L^2(\Omega)}+\sum_{i=1}^n \gamma_i(v_{1,i},v_{2,i})_{L^2(0,T)},
\end{split}
\end{equation}
where we used that $z_i = G'(u)v_i,$ for $i=1,2$, and $w=G''(u)[v_1,v_2]$. Moreover, we invoked the terminal condition for $p$ from the adjoint equation \eqref{adjoint}.
This form of $F''$ includes the solution $w$ that implicitly depends on the increments $v_1,\, v_2$ via
the partial differential equation \eqref{varforw}. Now we proceed in the same way that we used to show
equation \eqref{zp}. We insert $p$ as test function in equation \eqref{varforw} and  $w$ as test function in the adjoint equation \eqref{adjointeq}. After some integration by parts, we arrive at the relation
\begin{equation*}
{\alpha_Q} \big(G(u)-\rho_Q,w\big)_{L^2(Q)} + {\alpha_\Omega}\big(p(T),w(T)\big)_{L^2(\Omega)}=- \Big(\nabla p \, , \, b \otimes v_1 z_2
+ b\otimes v_2 z_1\Big)_{L^2(Q)}
\end{equation*}
that yields \eqref{zp}. Inserting this result in \eqref{F''exp}, we verify the claim.
\end{proof}

We will state our second order optimality conditions in terms of the {\em quadratic form}
\begin{equation}
\label{D2F}
F''(u)v^2 = \iint_Q \Big[{\alpha_Q}z^2-2\nabla p\cdot (z b \otimes v)\Big]dxdt + \sum_{i=1}^n \gamma_i\int_0^T v_i^2 dt + {\alpha_{\Omega}} \int_\Omega z(T)^2 dx,
\end{equation}
that is obtained from the general form \eqref{F''} of $F''$ by the choice $v_1=v_2=v$.

\subsection{Critical directions}

For the sets defined in \eqref{cones}, the following characterization is easily obtained.
\begin{proposition}[Characterization of the critical cone]
\begin{equation}
\label{CharaCu}
C(\bar{u})= \left\{
\begin{split}
&v\in L^2(0,T;\cR^n): \,\text{a.e. on } [0,T]  \text{ and for } i=1,\dots,n,\\
& v_i =
\left\{
\begin{split}
\geq 0\quad &\text{if } \bar{u}_i=u^{\min}_i\\
\leq 0\quad &\text{if } \bar{u}_i=u^{\max}_i\\
0\quad &\text{if } \bar\Phi_i\neq 0
\end{split}
\right.
\end{split}
\right\},
\end{equation}
where $\bar\Phi$ was defined in \eqref{Phi}.
\end{proposition}

\begin{proof}
We follow essentially the lines of the proof in \cite[p. 273]{CasasTroeltzsch2012}.
Let us use $K(\bar{u})$ to denote the set on the r.h.s. of \eqref{CharaCu}.
\\
\indent Take $v\in C(\bar{u})$. Then, there exists $\{v_k\} \subseteq S(\bar{u})$
such that $v_k \to v$ in $L^2(0,T;\cR^n).$ By definition of $S(\bar{u})$ one has, necessarily, that  $v_{k,i}\geq0$ if $\bar{u}_i = u^{\min}_i$ and $v_{k,i} \leq0$ if $\bar{u}_i=u^{\max}_i,$ a.e. on $[0,T]$ and for every $i=1,\dots,n.$
Clearly, this property is preserved for the limit in $L^2(0,T;\cR^n),$ so that it holds also for $v.$
From this fact and the first order necessary condition of \eqref{signPhi}, we deduce that
$$
0=F'(\bar{u})v = \sum_{i=1}^n \int_0^T \bar\Phi_i(t)v_i(t) dt = \sum_{i=1}^n \int_0^T \left| \bar\Phi_i(t)v_i(t) \right|dt,
$$
which implies that $v_i=0$ if $\bar\Phi_i \neq 0.$ This proves that $C(\bar{u}) \subseteq K(\bar{u})$.
In order to prove the converse inequality, take $v\in K(\bar{u})$, and define, for each positive integer $k$ and for each $i=1,\dots,n,$
\begin{equation}
\label{vki}
v_{k,i}(t):=
\left\{
\begin{split}
&0 \quad \text{if } u^{\min}_i(t)<\bar{u}_i(t)<u^{\min}_i(t)+\frac1k \, \text{ or } \, u^{\max}_i(t)-\frac1k < \bar{u}_i(t) < u^{\max}_i(t),\\
&\mathbb{P}_{[-k,k]}(v_i(t))\quad \text{otherwise},
\end{split}
\right.
\end{equation}
where $\mathbb{P}_{[-k,k]}$ denotes the pointwise projection onto the interval $[-k,k]$.
It easily follows that $u_k:=\bar{u}+\lambda_k v_k$ belongs to $\mathcal{U}_{\rm ad}$ for
$\displaystyle\lambda_k := \frac{1}{k^2}.$
Hence $v_k \in S(\bar{u}).$ Thus $K(\bar{u}) \subseteq  {\rm cl}_{L^2}(S(\bar{u})) \subseteq C(\bar{u}),$ which concludes the proof.
\end{proof}

\subsection{Second order necessary optimality conditions}
In this subsection we confirm the conditions \ref{C1} and {\rm \ref{C2}} for \eqref{P} and apply the second order necessary condition of Theorem \ref{SONC2012} to our optimal control problem.

We impose the following additional hypothesis, that is assumed to hold throughout the remainder of the article.
\begin{assumption}
\label{Hypb}
The function $b$ belongs to $W^{1,\infty}(\Omega;\mathbb{R}^n)$ and, for a.a. $x \in \Gamma,$ it holds
\[
\big(b(x) \otimes u\big) \cdot n(x) = 0 \quad  \text{for all } u \in \mathbb{R}^n.
\]
\end{assumption}
\begin{remark}
 Assumption \ref{Hypb} above is fulfilled in particular, if $b \in W^{1,\infty}_0(\Omega;\mathbb{R}^n)$.
 Moreover, it holds if $b(x) \cdot n(x) = 0$ a.e. on $\Gamma$ and the control is of the form $u(t)= \mathfrak{u}(t) \big(1,\dots,1\big) \in \mathbb{R}^n,$ with $\mathfrak{u}$ being a scalar function.
 \end{remark}

Let us first prove the following technical result.

\begin{lemma}
\label{LemDeltap}
Given $u_1,u_2\in L^\infty(0,T;\cR^n)$ with associated adjoint states $p_1,  p_2$, respectively, we have
\begin{equation}
\label{deltap}
\|p_2- p_1\|_{W(0,T)}+\|p_2- p_1\|_{C([0,T];L^2(\Omega))} \leq C\|p_1\|_{C([0,T];L^2(\Omega))}\|u_2-u_1\|_2,
\end{equation}
where $C$ depends continuously on $\|u_2\|_\infty.$
\end{lemma}
\begin{proof}
Setting $\delta p:= p_2-p_1,$ we obtain
\begin{equation}
\label{adjointdiff}
\begin{split}
- \partial_t \delta p - \nu \Delta \delta p+B[u_2]\cdot \nabla\delta p &=   {\alpha_Q}( \rho_2-\rho_1)  - \big({b \otimes (u_2- u_1)}\big)\cdot \nabla {p_1},\\
 \delta p(T) &= {\alpha_\Omega}( \rho_2(T)-\rho_1(T))\quad \text{on } \Omega,\\
 \partial_n(\delta p) &= 0\quad \text{on } \Sigma.
\end{split}
\end{equation}
We can write the first two equations as
\begin{equation*}
\begin{split}
- \partial_t \delta p + a[u_2](\cdot,\delta p) &= f,\\
 \delta p(T) &={\alpha_\Omega}( \rho_2(T)-\rho_1(T)),
 \end{split}
\end{equation*}
where $f:= {\alpha_Q} (\rho_2-\rho_1) - \big(b \otimes (u_2- u_1)\big)\cdot \nabla {p_1}$.
We next obtain  estimates that are analogous to \eqref{estH}-\eqref{estrhodot}. For $\varphi \in H^1(\Omega),$ we have
\[
\begin{split}
&\iint_Q f \varphi dxdt = \iint_Q   \Big(\alpha_Q(\rho_2-\rho_1) - \big(b \otimes (u_2- u_1)\big)\cdot \nabla {p_1}\Big) \varphi\,dxdt\\
&= \iint_Q  \Big( \alpha_Q(\rho_2-\rho_1)\varphi + p_1  {\rm div} \big(b \otimes (u_2- u_1)\varphi\big)\Big)  dxdt   -\iint_\Sigma  \varphi p_1 \big(b \otimes (u_2- u_1)\big) \cdot n ds dt\\
&= \iint_Q  \Big( \alpha_Q(\rho_2-\rho_1)\varphi + p_1  {\rm div} \big(b \otimes (u_2- u_1)\varphi\big)\Big)  dxdt  =: F(\varphi),
\end{split}
\]
where we used Assumption \ref{Hypb} in the second equality, and
where the functional $F$ belongs to $L^2(0,T;H^1(\Omega)^*)$.
The term
\[
q:= p_1 {\rm div} (b\otimes(u_2-u_1) \varphi)= p_1\big( \underbrace{{\rm div} (b\otimes(u_2-u_1))}_{L^2(0,T;L^\infty(\Omega))}\big)\varphi + p_1\underbrace{(b\otimes(u_2-u_1))}_{L^2(0,T;L^\infty(\Omega))} \cdot \nabla \varphi.
\]
belongs to $L^1(Q)$, because we have  $p_1 \in C([0,T];L^2(\Omega))$. It can be estimated by
\[
\|q\|_{L^1(Q)} \le C\, \|p_1\|_{C([0,T];L^2(\Omega))} \|u_2-u_1\|_2 \|\varphi\|_{H^1(\Omega)^*}
\]
Now the claim follows easily from Propositions \ref{Propdeltarho} and \ref{P5.1} by estimating $F$.
\end{proof}

\begin{proposition}
\label{PropA} The conditions {\rm \ref{C1}} and {\rm \ref{C2}} are fulfilled for problem  \eqref{P}.
\end{proposition}

\begin{proof}
 Recall the expressions for $F'$ and $F''$ given in \eqref{F'} and \eqref{D2F}.
 Let us first verify {\rm \ref{C1}}.
In view of the representation \eqref{F'exp},  the continuous extension of $F'$ is possible if  $ \Phi$
in \eqref{Phi} belongs to $L^2(0,T;\cR^n)$. Obviously, this holds true if the functions $ \rho \, b_i \partial_{x_i} p$ belong to $L^2(0,T;L^1(\Omega))$, for each $1,\dots,n.$ This, however follows from $\rho \in C([0,T];L^2(\Omega))$, $\partial_{x_i} p \in L^2(0,T;L^2(\Omega))$, and $b_i \in L^\infty(\Omega)$.

Next, we confirm the associated extension property of $F''(u)$. First of all, we mention that the mapping $v \mapsto z = G'(u)v$ is continuous from $L^2(0,T;\cR^n)$ to $W(0,T)$, hence also to $C([0,T];L^2(\Omega))\cap L^2(0,T;H^1(\Omega))$. Therefore, the only nontrivial term in the expression \eqref{F''} for $F''$ is the integral
\[
\iint_Q\nabla p \cdot (z_2 b \otimes v_1 + z_1 b\otimes v_2 )\, dxdt.
\]
In view of symmetry, it suffices to consider the term
\[
\iint_Q \nabla p \cdot (b \otimes v_1) z_2\, dxdt,
\]
and to show that the mapping
\[
(v_1,v_2) \longmapsto \iint_Q \nabla p \cdot (b \otimes v_1) z_2\, dxdt
\]
can be continuously extended to $L^2(0,T;\cR^n) \times L^2(0,T;\cR^n)$.
 This is true, since
 \[
 \begin{aligned}
\left|\iint_Q \nabla p \cdot( b \otimes v_1) z_2\, dxdt\right| &\le  \|\nabla p\|_{L^2(0,T;L^2(\Omega)^n)} \|b\|_{L^\infty(\Omega)^n} \|v_1\|_2 \|z_2\|_{{C([0,T];L^2(\Omega))}}\\
&\le C  \, \|p\|_{L^2(0,T;H^1(\Omega))}  \|v_1\|_2 \|v_2\|_2,
 \end{aligned}
 \]
 where we used the estimate \eqref{estimatez} for $z_2$. This is the desired extension property.
\vspace{1.5ex}

  Now we are going to confirm  {\rm \ref{C2}}. Take a sequence $\{ (u_k,v_k)\} \subset \mathcal{U}_{\rm ad} \times L^2(0,T;\cR^n)$ as it is required in {\rm \ref{C2}}, this is, with $u_k\to \bar{u}$ in $L^2(0,T;\cR^n)$ and $v_k \rightharpoonup v$ (weakly) in $L^2(0,T;\cR^n)$.
Let us prove that \eqref{limF'} holds.

We have
\begin{equation}
\label{C2der}
F'(u_k) v_k = \sum_{i=1}^n \int_0^T \underbrace{\left[ \int_\Omega \rho_k  b_i \frac{\partial p_k}{\partial x_i} dx + \gamma_i u_{k,i} + b_i \right]}_{\displaystyle\Phi_{k,i}} v_{k,i} dt,
\end{equation}
where $\rho_k := G(u_k)$ and $p_k$ is the solution of \eqref{adjointeq} associated to $(\rho_k,u_k).$
Define  $\bar{\rho}:=G(\bar{u})$, let $\bar p$ denote the corresponding adjoint state  and $\bar\Phi$ be the associated function defined in \eqref{Phi}.
In view of the estimate \eqref{deltap} applied to $p_k-\bar p$, we have  $p_k \to \bar p$ in $C([0,T];L^2(\Omega))$.  From \eqref{estdeltarho},  we get $\rho_k \to \bar{\rho}$ in $C([0,T];L^2(\Omega))$, then we deduce  that
$$
\left\|\int_\Omega \rho_k(\cdot)  b_i \frac{\partial p_k(\cdot)}{\partial x_i} dx-\int_\Omega \bar{\rho}(\cdot)  b_i \frac{\partial\bar p(\cdot)}{\partial x_i} dx\right\|_{L^2(0,T)} \to 0
$$
 for each $i=1,\dots,n.$
 Therefore, $\Phi_{k}$ given in \eqref{C2der} converges strongly in $L^2(0,T;\cR^n)$ to $\bar\Phi$.
  Consequently, thanks to the convergence of a pairing of weakly-strongly convergent sequences, we can pass to the limit in \eqref{C2der}, and thus \eqref{limF'} follows.
\end{proof}

By application of Theorem \ref{SONC2012}, we get the following result.

\begin{theorem}[Second order necessary optimality condition]
\label{SONC}
Let $\bar{u}$ be an $L^\infty$-local solution for \eqref{P}. Then
\begin{equation}
\label{inequalityNC}
F''(\bar{u})v^2 \geq 0\quad \text{for all } v\in C(\bar{u}).
\end{equation}
\end{theorem}

\begin{proof}
{In view of Theorem \ref{SONC2012} and Proposition \ref{PropA}}, it only remains to prove that $C(\bar{u}) = {\rm cl}_{L^2}D(\bar{u}).$ The inclusion $C(\bar{u}) \supseteq {\rm cl}_{L^2}D(\bar{u})$ follows easily {from the expression of $F'$ given in \eqref{F'exp}, since $\bar\Phi$ belongs to
$L^2(0,T;\cR^n)$ as observed in the proof of Proposition \ref{PropA}}.

Let us show the converse inclusion, $C(\bar{u}) \subseteq {\rm cl}_{L^2}D(\bar{u}).$ Take $v\in C(\bar{u})$, and define $v_k,$  for each positive integer $k$, as done in \eqref{vki}.
Hence, we have $v_k \in S(\bar{u}).$ By definition, we know that  $v_k \to v$ strongly in $L^2(0,T;\cR^n)$.
Moreover, since $|v_{k,i}(t)|\leq |v_i(t)|$ a.e. on $[0,T]$ for every $i=1,\dots,n,$ we have in addition that
$|\bar\Phi_i(t)v_{k,i}(t)|\leq |\bar\Phi_i(t)v_{i}(t)|=0$ a.e. on $[0,T].$
This implies $F'(\bar{u})v_k=0.$ Consequently, $v_k\in D(\bar{u})$ holds for every $k$ and, hence, $v\in {\rm cl}_{L^2}D(\bar{u}).$ This yields the desired inclusion.

Finally, the necessary condition \eqref{inequalityNC} follows from Theorem \ref{SONC2012}.
\end{proof}

\subsection{Second order sufficient conditions}

To apply the second order sufficient condition of the abstract Theorem \ref{SOSCCT}, we need higher regularity of $\rho$. For the remainder of the article, we additionally impose the following hypotheses, along with Assumption \ref{Hypb} introduced above.
\begin{assumption}[Requirements for higher regularity] \label{A6.1}
~

\begin{itemize}
\item[(i)]  The function $c$ has a potential $-V \in
W^{2,\infty}(\Omega)$ so that
\begin{equation}
\label{potential}
c = \nabla V
\end{equation}
or $c$ belongs to $W^{1,\infty}(\Omega,\mathbb{R}^n)$ and satisfies
\begin{equation}
\label{orthogonality}
c(x) \cdot n(x) = 0 \quad \mbox{a.e. on } \Gamma,
\end{equation}
\item[(ii)] the initial and the desired distributions $ \rho_0$ and $\rho_\Omega$ belong to $H^1(\Omega)$,
\item[(iii)]  and
$$
\min_{i=1,\dots,n} \gamma_i >0.
 $$
\end{itemize}
 \end{assumption}
 %
 \begin{remark} The assumption \eqref{potential} that $c$ has a potential  was imposed by Breiten {\em et al.} \cite{BreitenKunischPfeiffer2018} to gain higher regularity of $\rho$.
\end{remark}
Let us show that under the Assumptions \ref{Hypb} and \ref{A6.1} the state equation for $\rho$ can be transformed to a linear heat equation with homogeneous Neumann condition and right-hand side in $L^2(Q)$. This is the key for higher regularity of $\rho$.

To see this, we first consider the case where  Assumption \ref{Hypb} and condition \eqref{orthogonality} are fulfilled. Then,  for any control function $u \in L^\infty(0,T;\mathbb{R}^n),$  ${\rm div }(\rho(c+b\otimes u)) = {\rm div} (\rho B[u])\in L^2(Q)$ and
\[
- \int_\Omega \rho(c+ b\otimes u) \cdot  \nabla \varphi\, dx=\int_\Omega{\rm div }\big( \rho(c+b\otimes u) \big) \varphi \, dx,  \quad \forall
\varphi \in H^1(\Omega), \forall u \in \mathbb{R}^n,
\]
since  $\rho(c+b\otimes u) \cdot n = \rho B[u]\cdot n = 0$
on $\Gamma.$ Therefore, the state equation \eqref{FP},\eqref{initcond} reduces to
the linear heat equation
\begin{equation} \label{heateq}
\begin{split}
\partial_t \rho - \nu \Delta\rho & = \hat f,\\
\rho(0)&=\rho_0,\\
\nabla\rho\cdot n &=0,
\end{split}
\end{equation}
where
\begin{equation} \label{hatf1}
\hat f:= {\rm div} (\rho B[u]) \in L^2(Q).
\end{equation}

Next, we confirm the reduction to an equation of the form \eqref{heateq}, if $c$ does not
satisfy \eqref{orthogonality} but fulfils \eqref{potential}. In this case, we follow an idea of
\cite[proof of Proposition 2.1]{BreitenKunischPfeiffer2018}, and apply the transformation $w = \exp{(V/\nu)}\rho$. Notice that $w$ enjoys the same regularity as $\rho$, i.e. $w \in W(0,T)$. We find
\[
\nabla \rho = -\frac{1}{\nu} \exp{(-V/\nu)} \nabla V w +  \exp{(-V/\nu)}\nabla w=  \exp{(-V/\nu)} \, \left(-\frac{1}{\nu}c \, w +  \nabla w\right).
\]
Inserting this in the boundary condition of the state equation, we arrive at
\[
\begin{aligned}
0=\nu \nabla \rho \cdot n + \rho B[u] \cdot n =&\exp{(-V/\nu)} (-c \, w +  \nu \nabla w)\cdot {n}\\
&+ \exp{(-V/\nu)}w (c + b\otimes u)\cdot n = \nu \exp{(-V/\nu)}\nabla w\cdot {n}
\end{aligned}
\]
that holds true if, and only if, $\nabla w \cdot n= 0$ on $\Gamma$. Therefore, $w$ satisfies homogeneous Neumann boundary conditions. Moreover, inserting the transformation of $\rho$ in the state equation, we obtain, after some computations by the product rule applied to ${\rm div}\,\nabla  \big(\exp{(-V/\nu)}w)$ that $w$ satisfies the linear heat
equation \eqref{heateq} with right-hand side
\begin{equation}
\label{hatf2}
\hat f = - c \cdot \nabla w - \frac{1}{\nu} wc\cdot b\otimes u + {\rm div} (w b\otimes u).
\end{equation}
Since $w \in W(0,T)$ and $\nabla w \cdot n= 0$ on $\Gamma$, all summands above  belong to $L^2(Q)$.

Next, we provide a regularity result for equation \eqref{heateq} that can be found, for homogeneous Dirichlet boundary conditions and  a smooth boundary $\Gamma$, in  Evans \cite[Theorem 5, page 360]{Evans2010book}  or Ladyzhenskaya
{\em et al.} \cite{Ladyzhenskaya_etal1968}. For the convenience of the reader, we will prove the part on $C([0,T];H^1(\Omega))$-regularity in Lipschitz domains. It is difficult to find an associated reference, although the result is known. The proof was communicated to us by Joachim Rehberg (WIAS Berlin).

\begin{lemma}
\label{Lmaxparabreg}
If $\Omega$ is a bounded Lipschitz domain, $\hat f$ is an arbitrary function of $L^2(Q)$ and $\rho_0$ belongs to $H^1(\Omega)$,
then the unique weak solution $\rho$ of the linear heat equation \eqref{heateq} enjoys the higher regularity $\rho \in C([0,T];H^1(\Omega)).$
Moreover, the estimate
\begin{equation} \label{a_priori_est}
\|\rho\|_{W(0,T)} + \|\rho\|_{C([0,T];H^1(\Omega))} \le C\, \big(\|\hat f\|_{L^2(Q)} + \|\rho_0\|_{H^1(\Omega)} \big)
\end{equation}
is satisfied with some $C > 0$ not depending on $\hat f$ neither on $\rho_0$.

\if{
\textcolor{green}{
If $\Omega$ is in addition of class $C^{1,1}$ or convex, then
the higher regularity $\rho \in W(0,T;H^2(\Omega),L^2(\Omega))$ is fulfilled, and the associated estimate
\[
\|\rho\|_{L^2(0,T;H^2(\Omega))} \le C\, (\|\hat f\|_{L^2(Q)} + \|\rho_0\|_{H^1(\Omega)})
\]
holds with some $C > 0$ not depending on $\hat f$ neither on $\rho_0$.}
}\fi
\end{lemma}

\begin{proof} The Laplace operator $\Delta$ is generated in $L^2(\Omega)$ by the bilinear form
\[
(\psi,\varphi) \mapsto  \int_\Omega \nabla \psi \cdot \nabla \varphi \, dx,
\]
defined on $H^1(\Omega) \times H^1(\Omega)$, cf. \cite[chpt. VI]{Kato1980} or \cite[chpt. 1]{Ouhabaz2005}. This operator, induced in
$L^2(\Omega)$, is non-positive and self-adjoint, hence it generates an analytic semigroup in $L^2(\Omega)$.

It is known that the negative of any linear operator in a Hilbert space that generates an analytic semigroup has maximal parabolic regularity. In our setting, this means the following: if $S$  is an arbitrary bounded or unbounded real interval, for every $f \in L^2(S;L^2(\Omega))$ and each initial value $y_0$ from the real interpolation space $\big(L^2(\Omega),{\rm dom}_{L^2(\Omega)}\Delta\big)_{\frac{1}{2},2},$ there exists a unique $y \in W^{1,2}(S;L^2(\Omega)) \cap L^2(S;{\rm dom}_{L^2(\Omega)}\Delta)$ such that
\[
\partial_t y -  \Delta y = f, \quad  y(0) = y_0.
\]
For this existence result, we refer to \cite[chpt. 1.3]{Ashyralyev_Sobolevski1994}, in particular to chpt. 1.3.3. Moreover, one has (see e.g. \cite[chpt. III.4.10]{Amann1995})
\begin{equation} \label{embedding0}
W^{1,2}(S;L^2(\Omega)) \cap L^2(S,{\rm dom}_{L^2(\Omega)}\Delta) \hookrightarrow C\big(\bar{S}; (L^2(\Omega),{\rm dom}_{L^2(\Omega)}\Delta)_{\frac{1}{2},2}\big).
\end{equation}

The operator $I -  \Delta$ is positive and self-adjoint with lower spectral bound $1$. Therefore, thanks to the spectral theory, its pure imaginary powers exist as bounded operators in $L^2(\Omega)$. In view of this, one has (see \cite[chpt. 1.18.10]{Triebel1978})
\[
\begin{aligned}
\big(L^2(\Omega),{\rm dom}_{L^2(\Omega)}\Delta\big)_{\frac{1}{2},2}& = \big(L^2(\Omega),{\rm dom}_{L^2(\Omega)}(-\Delta+I)\big)_{\frac{1}{2},2}\\
& = \big[L^2(\Omega),{\rm dom}_{L^2(\Omega)}(-\Delta+I)\big]_{\frac{1}{2}} = {\rm dom}_{L^2(\Omega)}
(-\Delta+I)^{\frac{1}{2}}.
\end{aligned}
\]
Since $\Delta$ is self-adjoint, the space ${\rm dom}_{L^2(\Omega)}
(-\Delta+I)^{\frac{1}{2}}$ is equal to its form-domain, i.e. $H^1(\Omega)$, \cite[chpt. VI 2.6]{Kato1980}.
Therefore, \eqref{embedding0} implies
\begin{equation} \label{embedding}
W^{1,2}(S;L^2(\Omega)) \cap L^2(S,{\rm dom}_{L^2(\Omega)}\Delta) \hookrightarrow C(\bar S;H^1(\Omega)).
\end{equation}
Exactly the same arguments are true for $\nu \Delta$ instead of $\Delta$, except that the lower spectral
bound is now estimated from below by  $\nu$. These statements are valid in arbitrary domains $\Omega$.

The estimate \eqref{a_priori_est} is shown as usual by considering a graph norm.
\end{proof}

\begin{remark} \label{Rem6.1} In addition to the continuous embedding \eqref{embedding}, also the embedding
\[W^{1,2}(S;L^2(\Omega)) \cap L^2(S;{\rm dom}_{L^2(\Omega)}\Delta) \hookrightarrow C^{1/2}(\bar{S};L^2(\Omega))
\]
is continuous. Therefore, by interpolation arguments, cf. Triebel \cite{Triebel1978}  the same holds for the embedding
\[
W^{1,2}(S;L^2(\Omega)) \cap L^2(S;{\rm dom}_{L^2(\Omega)}\Delta) \hookrightarrow C^{\alpha}(\bar{S};H^{1/2-\alpha}(\Omega)),
\]
with $0 < \alpha < \frac{1}{2}$. In view of this, the embedding
\[
W^{1,2}(S;L^2(\Omega)) \cap L^2(S;{\rm dom}_{L^2(\Omega)}\Delta) \hookrightarrow C(\bar{S};L^2(\Omega)),
\]
is compact. Therefore, for the heat equation \eqref{heateq}, the mapping $\hat f \mapsto \rho$ is compact from $L^2(Q)$ to $C([0,T];L^2(\Omega)) $.
\end{remark}

The next result follows from Breiten, Kunisch and Pfeiffer \cite[Proposition 6.1]{BreitenKunischPfeiffer2018}, who proved $W(0,T;H^2(\Omega),L^2(\Omega))$-regularity in $C^2$ domains. We provide a different approach and obtain $C([0,T];H^1(\Omega))$-regularity for a Lipschitz domain $\Omega$.

\begin{theorem}
\label{higher_reg1}
If Assumptions \ref{Hypb} and \ref{A6.1} are fulfilled,
 then, given any control $u\in L^\infty(0,T;\mathbb{R}^n),$ the unique weak solution $\rho$ of the Fokker-Planck equation \eqref{FP},\eqref{initcond} belongs to $C([0,T];H^1(\Omega)).$
\end{theorem}
\begin{proof}We have pointed out above that the state $\rho \in W(0,T)$ solves the linear heat equation \eqref{heateq} with a right-hand side $\hat f \in L^2(Q)$. Depending on the particular Assumption
\ref{A6.1}-(i) on $c$, the function $\hat f$ is given by \eqref{hatf1} or by \eqref{hatf2}.  Now the claimed regularity follows from Lemma \ref{Lmaxparabreg}.
\end{proof}

\begin{proposition}
\label{PropC3}
Under Assumptions \ref{Hypb} and \ref{A6.1}, problem \eqref{P} satisfies condition {\rm \ref{C3}}.
\end{proposition}

\begin{proof}
Consider a sequence $\{(u_k,v_k)\} \subset \mathcal{U}_{\rm ad} \times L^2(0,T;\cR^n)$ with $u_k\to \bar{u}$ in $L^2(0,T;\cR^n)$ and $v_k \rightharpoonup v$ weakly in $L^2(0,T;\cR^n)$ as in  {\rm \ref{C3}}. We will proceed in several steps.

{\bfseries 1) Proof of \eqref{liminfD2J}.}

{\em 1a)  Convergence of ${\{z_k\}}$.}
First, we show that $\{z_k\}$ has a well-defined limit. To this aim, we write the linearized equation for $z_k$
in the following form:
\begin{multline}
\label{lineareqk1}
\iint_Q \left( - z_k \partial_t \varphi + (\nu\nabla z_k+z_k B[\bar{u}])\cdot \nabla \varphi \right) dx\, dt\\
\hspace{2cm}= \iint_Q z_k b \otimes (\bar{u}-u_k)\cdot \nabla \varphi \, dxdt - \iint_Q\rho_k (b\otimes v_k) \cdot \nabla \varphi \, dx\,dt,
\end{multline}
for all $\varphi \in W^{1,1}_2(Q)$ with $\varphi(\cdot,T)=0$. All controls $u_k$ belong to $\mathcal{U}_{\rm ad}$ and are therefore uniformly bounded in $L^\infty(0,T;\cR^n)$. Hence, thanks to Proposition \ref{Propestimatez}, the functions $z_k$ vary in a bounded set of $L^\infty(0,T;L^2(\Omega))$. In view of this and since $u_k \to \bar{u}$ in $L^2(0,T;\cR^n)$, the term  $z_k b \otimes (\bar{u}-u_k)$ under the first integral in the r.h.s. of \eqref{lineareqk1} tends to zero in $L^2(Q)^n$. Moreover, it is easy to confirm that
\[
\rho_k (b\otimes v_k) \rightharpoonup \bar\rho (b\otimes v) \mbox{ in } L^2(Q)^n,\quad \text{ when } k \to \infty,
\]
since the $\rho_k$ are bounded in $C([0,T];L^2(\Omega))$.
Therefore, the functional in the r.h.s. of \eqref{lineareqk1} generated by these terms converges weakly in $L^2(0,T;H^1(\Omega)^*)$. The solution operator associated with the differential operator of the l.h.s. is linear and continuous from $L^2(0,T;H^1(\Omega)^*)$
to $W(0,T)$, hence it preserves weak convergence, so that $\{z_k\}$ is weakly convergent to some $z$ in $W(0,T)$ that  satisfies the equation
\begin{equation}
\label{lineareqk2}
\iint_Q \left( - z \partial_t \varphi + (\nu\nabla z+zB[\bar{u}])\cdot \nabla \varphi \right) dx\, dt
= - \iint_Q\bar\rho (b\otimes v) \cdot \nabla \varphi \, dx\,dt \quad
\end{equation}
for all $\varphi$ mentioned above. Hence $z$ is the linearized state associated to $\bar\rho$ in the direction $v.$

In order to show \ref{C3}, let us take a subsequence $\{(u_{k_j},v_{k_j})\}$ such that
$$
\liminf_{k\to \infty} F''(u_k) v_k^2 = \lim_{j\to \infty} F''(u_{k_j})v_{k_j}^2.
$$
The corresponding subsequence $\{z_{k_j}\}$ converges weakly to $z$ in $W(0,T)$ and then, in view of Aubin-Lions' Lemma, it contains a subsequence $\{z_{k_{j_\ell}}\}$ that converges strongly to $z$, i.e.
\begin{equation}\label{aux1}
z_{k_{j_\ell}} \to z \quad \text{in } L^2(Q).
\end{equation}
In order to simplify the notation, let us use the subindex $k$ for this subsequence.
One has
\begin{equation}
\label{F2k}
F''(u_k) v_k^2 = \iint_Q \Big[ z_k^2-2\nabla p_k \cdot ( z_k  b\otimes v_k)\Big] dxdt + \sum_{i=1}^n \gamma_i\int_0^T v_{k,i}^2 dt + \int_\Omega z_k(T)^2 dx,
\end{equation}
for $z_k$  satisfying \eqref{lineareqk1}.

{\em 1b) Convergence of ${\{p_k\}}$.}
Next we prove that $\nabla p_k$ tends to $\nabla \bar{p}$ in $C([0,T];L^2(\Omega)^n).$ From the adjoint equation \eqref{adjoint}, we get, for a generic $p$ associated to some state-control pair $(\rho,u),$
\begin{equation}
\label{adjoint2}
\begin{split}
- \partial_t p - \nu \Delta p  &= {\alpha_Q}(\rho-\rho_Q)-B[u]\cdot \nabla p \quad \text{in } Q,\\
 p(T) &={\alpha_\Omega}(\rho(T)-\rho_\Omega) \quad \text{ in } \Omega,\\
 \partial_np &= 0 \quad \text{on } \Sigma.
\end{split}
\end{equation}
The r.h.s. ${\alpha_Q}(\rho-\rho_Q)-B[u]\cdot \nabla p$ belongs to $L^2(Q)$. Moreover, thanks to Assumption \ref{A6.1}, (ii), ${\alpha_\Omega}( \rho(T)-\rho_\Omega)$ is in  $H^1(\Omega).$ We can apply Lemma \ref{Lmaxparabreg} and get that $\bar p$ and $p_k$ are in $C([0,T];H^1(\Omega)).$  Set $\delta p_k := p_k-\bar{p}$. Then, from \eqref{adjointdiff}, we obtain,
\begin{equation}
\label{deltapk}
\begin{split}
- \frac{\partial\delta p_k}{\partial t} - \nu \Delta \delta p_k &=  \underbrace{{\alpha_Q}( \rho_k-\bar\rho)  - ({b \otimes (u_k- \bar u)}\big)\cdot \nabla \bar p - B[u_k]\cdot \nabla\delta p_k}_{\displaystyle \hat g_k},\\
 \delta p_k(T) &= {\alpha_\Omega}\big( \rho_k(T) - \bar\rho(T) \big),\quad \text{on } \Omega,\\
 \partial_n(\delta p_k) &= 0,\quad \text{on } \Sigma.
\end{split}
\end{equation}
This is again a linear heat equation with homogeneous Neumann condition and r.h.s $\hat{g}_k$ in $L^2(Q)$. Hence, Lemma \ref{Lmaxparabreg}, applied to the adjoint equation, guarantees that $\delta p_k$  satisfies the estimate
\begin{equation}
\label{estimatedeltapk}
\|\delta p_k\|_{W(0,T)} + \|\delta p_k\|_{C([0,T];H^1(\Omega))} \leq C\Big(\|\hat g_k\|_{L^2(Q)} + \| \rho_k(T) - \bar\rho(T)\|_{H^1(\Omega)}\Big).
\end{equation}
 For $\|\hat g_k\|_{L^2(Q)},$ we have
\begin{equation}
\label{hatgk}
\begin{split}
\|\hat g_k\|_{L^2(Q)} \leq &
\,\alpha_Q \|\rho_k-\bar\rho\|_{L^2(Q)} + \sum_{i=1}^n \|b_i\|_{L^\infty(\Omega)}\|u_{k,i}-\bar{u}\|_2\left\| \frac{\partial \bar{p}}{\partial x_i} \right\|_{C([0,T];L^2(\Omega))}\\
&+\sum_{i=1}^n \big( \|c_i\|_{L^\infty(\Omega)}+\|b_i\|_{L^\infty(\Omega)}\|u_{k,i}\|_\infty\big)\left\| \frac{\partial \delta p_k}{\partial x_i} \right\|_{L^2(Q)}
\end{split}
\end{equation}
From Lemma \ref{LemDeltap}, we know $\nabla\delta p_k \to 0$ in $L^2(Q)^n$ since $u_k\to \bar{u}$ in $L^2(0,T;\cR^n)$ hence, the r.h.s. of latter display tends to 0 and $\|\hat g_k\|_{L^2(Q)} \to 0.$ In this way, we have shown the desired convergence
\begin{equation}
\label{aux2}
\nabla p_k \to \nabla \bar p \quad \text{in } C([0,T];L^2(\Omega)^n).
\end{equation}

{\em 1c) Convergence of ${\{\rho_k\}}$ in $C([0,T];H^1(\Omega))$.}
It remains to show that $\| \rho_k(T) - \bar\rho(T)\|_{H^1(\Omega)}$ in \eqref{estimatedeltapk} converges to 0.

Let us first discuss the case where $c$ obeys the condition \eqref{orthogonality} of Assumption \ref{A6.1}. Here we get, from Lemma \ref{Lmaxparabreg} applied to $\rho_k$,
\[
\|\rho_k\|_{W(0,T)} + \|\rho_k\|_{C([0,T];H^1(\Omega))} \leq C \big( \|{\rm div}(\rho_k B[u_k])\|_{L^2(Q)} + \|\rho_0\|_{H^1(\Omega)}\big).
\]
Since $\{\rho_k\}$ is bounded in $W(0,T)$ and $\{u_k\}$ is bounded in $L^\infty(0,T;\cR^n),$ the sequence $\{{\rm div}(\rho_k B[u_k])\}$ is bounded in $L^2(Q)$.
Consequently, from previous display we deduce that $\{\rho_k\}$ is bounded in $C([0,T];H^1(\Omega)).$ For $\rho_k-\bar\rho$ we have the following equation
\[
\frac{\partial(\rho_k-\bar\rho)}{\partial t} - \nu \Delta(\rho_k-\bar\rho) =\underbrace{{\rm div}(\rho_k B[u_k])- {\rm div} (\bar\rho B[\bar{u}])}_{\displaystyle\hat f_k},
\]
with initial condition $(\rho_k-\bar\rho)(0)=0$ and the homogeneous Neumann condition $\nabla(\rho_k-\bar\rho) \cdot n =0.$ For $\hat f_k$ we have
\begin{equation*}
\begin{split}
&\|\hat f_k\|_{L^2(Q)} =  \big\|{\rm div}(\rho_k B[u_k])- {\rm div} (\bar\rho B[\bar{u}]) \pm {\rm div}(\rho_k B[\bar{u}])\big\|_{L^2(Q)} \\
& \leq  \left\|{\rm div}\big(\rho_k (B[u_k]-B[\bar{u}])\big)+ {\rm div} \big((\rho_k-\bar\rho) B[\bar{u}]\big) \right\|_{L^2(Q)} \\
& \leq \sum_{i=1}^n \left[ \left( \left\| \frac{\partial \rho_k}{\partial x_i} \right\|_{C([0,T];L^2(\Omega))} \|b_i\|_{L^\infty(\Omega)}
+ \|\rho_k\|_{C([0,T];L^2(\Omega))} \left\| \frac{\partial b_i}{\partial x_i} \right\|_{L^\infty(\Omega)} \right)\|u_{k,i}-\bar{u}_i\|_2 \right.\\
&+\left. \left\|\frac{\partial (\rho_k-\bar\rho)}{\partial x_i}\right\|_{L^2(Q)} \|c_i+b_i\bar{u}_{k,i}\|_{L^\infty(Q)} + \|\rho_k-\bar\rho\|_{L^2(Q)} \left\|\frac{\partial c_i}{\partial x_i}+\frac{\partial b_i}{\partial x_i} \bar{u}_{k,i}\right\|_{L^\infty(Q)}\right].
\end{split}
\end{equation*}
Since $\{\rho_k\}$ is bounded in $C([0,T];H^1(\Omega)),$ $u_k \to \bar{u}$ in $L^2(0,T;\mathbb{R}^n)$ and $\rho_k \to \bar\rho$ in $W(0,T),$ the r.h.s. of latter display tends to 0, and hence $\|\hat f_k\|_{L^2(Q)}\to 0.$ We then have that $\rho_k \to \bar\rho$ in $C([0,T];H^1(\Omega)).$ From \eqref{estimatedeltapk} and \eqref{hatgk}, we get that $\delta p_k \to 0$ in $C([0,T];H^1(\Omega))$ as desired.

Now we briefly discuss the other case, where $c$ fulfils condition \eqref{potential} of Assumption \ref{A6.1}.
Here, we invoke again the transformation $w_k = \exp{(V/\nu)}\rho_k$ that we used to transform the Fokker-Planck equation to one with homogeneous boundary condition. Thanks to \eqref{hatf2}, the functions $w_k$ satisfy the linear heat equation \eqref{heateq} with r.h.s.
\begin{equation*}
\label{hatf3}
\tilde f_k := - c \cdot \nabla w_k - \frac{1}{\nu} w_k c\cdot b\otimes u_k + {\rm div} (w_k b\otimes u_k).
\end{equation*}
It is obvious that $w_k$ converges to $\bar w = \exp{(V/\nu)}\bar \rho$ in $W(0,T)$. This yields
$\tilde f_k \to \bar f$ in $L^2(Q)$, where
\[
\bar f := - c \cdot \nabla \bar w - \frac{1}{\nu} \bar w c\cdot b\otimes \bar u + {\rm div} (\bar w b\otimes \bar u).
\]
Setting now $\hat f_k := \tilde f_k - \bar f$, we have $\hat f_k \to 0$ in $L^2(Q)$. In the same way as above, Lemma \ref{Lmaxparabreg} yields that $w_k - \bar w \to 0$ in  $C([0,T];H^1(\Omega))$. This transfers to the sequence $\rho_k - \bar \rho$.

\vspace{5pt}

{\em 1d) Proof of \eqref{liminfD2J}.}  We are now ready to check that \eqref{liminfD2J} is satisfied.  Let us recall the expression of $F''(u_k)v_k^2$ given in \eqref{F2k}. For the bilinear term $\nabla p_k \cdot (z_k b\otimes v_k),$ we deduce from \eqref{aux2} that $
\nabla p_k \to \nabla \bar{p}$ (strongly) in $L^\infty(0,T;L^2(\Omega)^n),$ and from \eqref{aux1} that $z_k \to z$ (strongly) in $L^2(Q)$. Moreover, we have $v_k \rightharpoonup v$ in $L^2(0,T;\cR^n).$ Therefore, the product $z_k b\otimes v_k$ converges weakly in $L^1(0,T;L^2(\Omega)^n)$.
Thus, we conclude that
\begin{equation*}
\begin{split}
&\iint_Q \Big(\nabla p_k \cdot (z_k b\otimes v_k) - \nabla \bar{p} \cdot (z b\otimes v) \Big)dx dt\\
&= \iint_Q \Big( \underbrace{(\nabla p_k-\nabla \bar{p})}_{L^\infty(0,T;L^2(\Omega)^n)} \cdot \underbrace{(z_k b\otimes v_k)}_{L^1(0,T;L^2(\Omega)^n)}\\
&\qquad \qquad  + \nabla \bar{p} \cdot \big((z_k-z) b\otimes v_k\big) + \nabla \bar{p} \cdot (z b\otimes (v_k-v)) \Big)dx dt
\end{split}
\end{equation*}
tends to 0.

It remains to pass to the limit in the other terms of $F''(u_k)v_k^2.$
Notice that one has
$$\iint_Q z_k^2 dxdt  \to \iint_Q z^2 dxdt$$
in view of the strong convergence $z_k \to z$ in $L^2(Q)$. Moreover, due to the weak convergence $v_k \rightharpoonup v$ in $L^2(0,T;\mathbb{R}^n),$ one gets
$$
\sum_{i=1}^n \gamma_i \int_0^T v_i^2 dt \leq \liminf_{k\to \infty} \sum_{i=1}^n \gamma_i \int_0^T v_{k,i}^2 dt.
$$
Since $z \mapsto z(\cdot,T)$ is linear and continuous from $W(0,T)$ to $L^2(\Omega),$ it is also weakly lower semicontinuous and, therefore,
$$
\int_\Omega z(T)^2 dx \leq \liminf_{k\to \infty} \int_\Omega z_k(T)^2 dx.
$$
This concludes the proof of \eqref{liminfD2J}.

\vspace{7pt}

{\bfseries 2) Proof of \eqref{liminfD2J0}.}
It remains to show \eqref{liminfD2J0}. So now we assume that $v=0,$ this is $v_k  \rightharpoonup 0$ in $L^2(0,T;\mathbb{R}^n).$ Thus, the corresponding limit $z$ of $\{z_k\}$ is also 0. We get
\begin{equation}
\begin{split}
\liminf_{k\to \infty} F''(u_k)v_k^2 & = \liminf_{k\to \infty} \left[ \sum_{i=1}^n\gamma_i \int_0^T v_{k,i}^2 dt + \int_\Omega z_k(T)^2 dx \right]\\
& \geq \left( \min_{i=1,\dots,n} \gamma_i \right) \liminf_{k\to \infty} \|v_k\|^2_2.
\end{split}
\end{equation}
This yields \eqref{liminfD2J0} for $\displaystyle\Lambda :=  \min_{i=1,\dots,n} \gamma_i > 0,$ and then the proof is concluded.

\if{We have
\begin{equation}
\iint_Q p_k{\rm div} ( z_k b \otimes v_k) dxdt = \int_0^T\int_{\Gamma} p_kz_k b \otimes v_k \cdot n ds dt - \iint_Q \nabla p_k \cdot (b\otimes v_k)z_k dx dt.
\end{equation}
}\fi

\end{proof}

\begin{theorem}[Second order sufficient condition]
\label{SOSC}
Let $\bar{u} \in \mathcal{U}_{\rm ad}$ be such that the first order condition \eqref{FirstAbstract} is satisfied, and
\[
F''(\bar{u})v^2 >0,\quad \text{for all } v \in C(\bar{u}) \backslash \{0\}.
\]
Under the assumptions \ref{Hypb} and \ref{A6.1}, there exist $\varepsilon>0$ and $\delta >0$ such that the quadratic growth condition
\[
J(\bar{u}) + \frac{\delta}{2} \|u-\bar{u}\|_2^2 \leq F(u),\quad  \text{for all } u \in \mathcal{U}_{\rm ad} \cap B^2_\varepsilon(\bar{u}),
\]
holds, hence $\bar u$ is locally optimal in the sense of $L^2(0,T;\mathbb{R}^n)$.
\end{theorem}

\begin{proof}
The result  follows from the application of  Theorem \ref{SOSCCT} to our optimal control problem.
The associated assumptions are fulfilled in view of Propositions \ref{PropA} and \ref{PropC3}.
\end{proof}

\begin{remark}We should mention that also  \cite[Corollary 2.6]{CasasTroeltzsch2012} on local uniqueness of $\bar u$ and \cite[Theorem 2.7]{CasasTroeltzsch2012} on using another critical cone can be directly transferred to our optimal control problem, because the general assumptions \ref{C1}-\ref{C3} are satisfied.
\end{remark}

\section*{Acknowledgements}

We thank Joachim Rehberg (WIAS Berlin) for communicating to us the proof of Lemma
\ref{Lmaxparabreg}. Moreover, we are grateful to Hannes Meinlschmidt (RICAM Linz) for his support in
proving the statements of Remark \ref{Rem6.1}l

\if{

\newpage

{\bf To do:
\begin{itemize}
\item Can we get improved regularity?
\item First order conditions for a tracking problem with constraints of the type
$$
a_i \leq u_i(x,t) \leq b_i.
$$
\item Second order conditions for the Tykhonov term $\nu>0$.
\item Second order conditions for the Tykhonov term $\nu=0$.
\item {\bf Can we do something for the bang-bang case? See Casas, Wachsmuth and Wachsmuth \cite{CasasWachsmuthWachsmuth2018}.}
\end{itemize}
}

\section{Some remarks on classic literature of Fokker-Planck equation}

Introduction to the Fokker-Planck equation: see Risken 1996 \href{https://link.springer.com/chapter/10.1007/978-3-642-61544-3_4}{link}

\bigskip

Derivation of the equation from the stochastic setting: see Crispin W.. Gardiner. Handbook of stochastic methods for physics, chemistry and the natural sciences. Springer-Verlag, 1994.

}\fi

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv


text:
Q: Unexpected estimate in Gamma GLM summary output I have a question. How on earth is it possible to have a negative estimate for a form of a nominal variable (two forms: "HM" and "LM") when it should be positive? I'm modelling a positive continuous variable (rarefied species richness) by using a single nominal predictor. See the output summary below and the box plot. How is this even possible?
P.S. This is my first gamma glm ever
Call:
glm(formula = Rarified.Richness ~ Macrophyte, family = "Gamma",
    data = beysehir.diversity.indices.macr)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-0.39143  -0.09652  -0.00911   0.07579   0.43107

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.082461   0.006693  12.320 2.28e-07 ***
MacrophyteLM -0.024028   0.010597  -2.267   0.0468 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.05929918)

    Null deviance: 0.87647  on 11  degrees of freedom
Residual deviance: 0.59465  on 10  degrees of freedom
AIC: 65.46

Number of Fisher Scoring iterations: 4

"LM" does have a higher mean value for the response. Shouldn't it has a positive estimate then?
P.S. KS-test and chisq tests both report a good fit for this model.
BONUS QUESTION: Do you think there is any other problem with this model?

A: You're using the inverse link function here, which is maybe not what you intended.  That explains why you get a negative coefficient for LM (inverse of a smaller linear predictor yields a larger expected value).
Inverse link is canonical for the gamma family. If you want to use log or identity link, you need to specify that explicitly.

redpajama_set_name: RedPajamaStackExchange


New Home in a Hurry available by McBride Homes. This popular Royal II floor plan will wow you from the moment you enter the beautiful two story foyer. The open main level floor plan provides the perfect amount of space to make your house the entertainment hub of your group. The Luxury Kitchen is a chief's dream and the spindled railing on the staircase allows light to flow freely throughout the home. This home offers other outstanding features such as Luxury Master Bath, 9' ceilings on the first floor, four window wall in the great room, Walkout basement, beautiful lot that backs to trees, upgraded flooring that runs throughout the home. This home is located in the beautiful Arbors of Rockwood Community where you will have access to 167 acers of common ground, a dog park, walking trails, fishing pond and several covered picnic areas. Our offices are open from 10-6 each day. Pix may be of similar home and not exact model. Some pix may have optional features not included in this home.