Optimal Control of Hughes’ Model for Pedestrian Flow via Local Attraction

Abstract

We discuss the control of a human crowd whose dynamics is governed by a regularized version of Hughes’ model, cf. Hughes (Transp Res Part B: Methodol 36(6):507–535, 2002. https://doi.org/10.1016/s0191-2615(01)00015-7). We assume that a finite number of agents act on the crowd and try to optimize their paths in a given time interval. The objective functional can be general and it can correspond, for instance, to the desire for fast evacuation or to maintain a single group of individuals. We provide an existence and regularity result for the coupled PDE-ODE forward model via an approximation argument, study differentiability properties of the control-to-state map, establish the existence of a globally optimal control and formulate optimality conditions.

1 Introduction

The starting point of this work is Hughes’ model for the movement of a (large) crowd of pedestrians introduced in [37]. Its unknowns are the density $\rho = \rho (t,x)$ and the potential $\phi = \phi (t,x)$ functions for $x \in \varOmega \subset \mathbb {R}^2$ and $t \in (0,T)$. With the space-time cylinder denoted by $Q_T :=(0,T) \times \varOmega$, the model reads

$$\begin{aligned}&\partial _t\rho - \nabla \cdot ({\rho \, f(\rho )^2 \, \nabla \phi }) = 0 \quad \text {in }Q_T, \end{aligned}$$

(1.1a)

$$\begin{aligned}&|\nabla \phi | = \frac{1}{f(\rho )}\quad \text {in } Q_T. \end{aligned}$$

(1.1b)

In the simplest case, the model is supplemented with homogeneous Dirichlet boundary conditions for $\phi$ and $\rho$. Due to the hyperbolic nature of the first equation, the boundary conditions for $\rho$ have to be posed in a suitable sense (i. e., on a generalized inflow part involving the function f). While many models for pedestrian dynamics are microscopic in the sense that they provide constitutive laws for the motion of each pedestrian (e. g., systems of ODEs or cellular automata or the social force model, cf. [5, 14, 35]), Hughes’ model starts from a macroscopic approach. It is based on the following three assumptions:

• (i) The velocity v of the pedestrians is determined by the density $\rho$ of the surrounding pedestrian flow and the behavioral characteristics of the pedestrians only. Denoting the movement direction by $u\in \mathbb {R}^2$ there holds

  $$\begin{aligned} v = f(\rho ) \, u , \quad |u| = 1 , \end{aligned}$$

  where $f(\rho )$ is a non-increasing function that attains the value one at $\rho =0$, is positive for $0< \rho < 1$ and zero at $\rho =1$. It models the effect that the larger the surrounding density, the slower the individual motion becomes.

• (ii) Pedestrians have a common sense of the task (called potential $\phi$), i. e., they aim to reach their common destination by

  $$\begin{aligned} u = -\frac{\nabla \phi }{|\nabla \phi |}. \end{aligned}$$

• (iii) Pedestrians seek to minimize their (accurately) estimated travel time, but modify their velocity to avoid high densities. The potential is thus a solution of the Eikonal equation

  $$\begin{aligned} |\nabla \phi |=\frac{1}{f(\rho )}. \end{aligned}$$

Combining these rules yields the model (1.1). For brevity we write

$$\begin{aligned} {\widetilde{\beta }(\rho ,\phi )} :=f(\rho )^2\,\nabla \phi . \end{aligned}$$

(1.2)

Due to the previous explanations it is clear that ${- {\widetilde{\beta }(\rho ,\phi )}}{|}_{(t,x)}$ is the direction of the individuals in the point $(t,x)\in Q_T$.

The analysis of this system is quite involved since the derivative of $\phi$, being the viscosity solution to (1.1b), has jump discontinuities on a set that depends on $\rho$ and is not known a priori. Since we are going to consider optimal control problems, we shall focus on a regularized version of the forward model (1.1), given by

$$\begin{aligned} \partial _t\rho - \nabla \cdot ({\rho \, {\widetilde{\beta }(\rho ,\phi )}})&= \varepsilon \, \varDelta \rho&\quad&\text {in }Q_T, \end{aligned}$$

(1.3a)

$$\begin{aligned} -\delta _1\,\varDelta \phi + |\nabla \phi |^2&= \frac{1}{f(\rho )^2+\delta _2}&\quad&\text {in }Q_T . \end{aligned}$$

(1.3b)

Here $\varepsilon$, $\delta _1$ and $\delta _2$ are positive and fixed regularization parameters. The model is supplemented with initial conditions

$$\begin{aligned} \rho (x,0) = \rho _0(x) \quad \text {in } \varOmega \end{aligned}$$

(1.4)

and mixed boundary conditions

$$\begin{aligned} \begin{aligned} -({\varepsilon \,\nabla \rho + \rho \, {\widetilde{\beta }(\rho ,\phi )}}) \cdot n&= \eta \,\rho ,&\phi&= 0{} & {} \text {on } {\varSigma _O } , \\ ({\varepsilon \,\nabla \rho + \rho \, {\widetilde{\beta }(\rho ,\phi )}}) \cdot n&= 0 ,&\nabla \phi \cdot n&= 0{} & {} \text {on } \varSigma _W , \end{aligned} \end{aligned}$$

(1.5)

where we assume that the boundary consists of parts which act as “walls” $\partial \varOmega _W$ not allowing the pedetrians to traverse, and a “truncation” or “outer boundary” $\partial \varOmega _O$, at which pedestrians can exit with a given outflow velocity $\eta >0$. We assume ${\partial \varOmega _O } \cup \partial \varOmega _W = \partial \varOmega$ and ${\partial \varOmega _O } \cap \partial \varOmega _W = \emptyset$ and set $\varSigma _W :=(0,T) \times \partial \varOmega _W$ and $\varSigma _O :=(0,T) \times \partial \varOmega _O$. An example domain is illustrated in Fig. 1.

Fig. 1

Illustration of a possible domain $\varOmega$ and the boundary parts $\partial \varOmega _{O }$ and $\partial \varOmega _{W }$. The gray region is excluded from $\varOmega$

Remark 1.1

(On the regularizations (I)) As multiple regularisations are carried out, let us briefly discuss their meaning in terms of modeling and mathematical necessity.

$\delta _1$::

The elliptic regularization in the Eikonal equation (1.3b) is such that in the limit $\delta _1 \rightarrow 0$, we would recover the unique viscosity solution. For $\delta _1=0$, $\nabla \phi$ becomes discontinuous on a set that depends on $\rho$. In this situation, existence for the original Hughes’ model is only known in one spatial dimension which would be a serious limitation for our application.

$\delta _2$::

While $\delta _2$ prevents blow-up of the right hand side in (1.3b), setting $\delta _1=0$, its effect on the modulus of the velocity $|{\widetilde{\beta }(\rho ,\phi )}|$ is as follows: For $\delta _2 = 0$, one obtains $| {\widetilde{\beta }(\rho ,\phi )}| = f(\rho )$ while otherwise one has

$$\begin{aligned} |\beta (\rho ,\phi )| = \frac{f(\rho )^2}{\delta _2 + f(\rho )} = f(\rho )\frac{f(\rho )}{\delta _2 + f(\rho )} . \end{aligned}$$

As $\delta _2$, the behavior of $|\beta (\rho ,\phi )|$ is essentially unchanged except for values of $\rho$ that are close to one (i.e. when $f(\rho )$ becomes very small). On the other hand, due to the properties of f, at $\rho =1$, $|\beta |$ vanishes and thus this modification is not effective in the dynamics.

$\varepsilon$::

The additional diffusion added to the evolution equation of the density serves two purposes. First, as we will see below, it allows for a continuous solution. This is essential as we will couple the system to an additional ODE later on that requires point evaluation of $\rho$. From the modeling point of view, the term introduces additional randomness in the motion of the pedestrians that seems a very reasonable addition to the directed motion of the convection term.

1.1 Optimal Control Problem

Our optimal control problem is based on the following scenario. We assume that there is a small, given number $M>0$ of agents (guides), who are able to locally influence the motion of pedestrians in their vicinity. Think, for instance, of tourist guides or marked security personnel at large sports events. The i-th agent’s position is described by a function $x_i(t) \in \mathbb {R}^2$, $i=1,\ldots , M$, and the interaction is modeled by a radially symmetric and decreasing kernel $K(x) = k(|x|)$, which enters (1.3a) as an additional potential given by

$$\begin{aligned} \phi _K({\varvec{x}};x) = \sum _{i=1}^M K({x - x_i(t)})= \sum _{i=1}^M k({|x - x_i(t)|}). \end{aligned}$$

(1.6)

Here and throughout, ${\varvec{x}}= {\varvec{x}}(t) = (x_1(t), \ldots , x_M(t))^T$ is the function collecting all agent positions. Furthermore, we need to modify the model to insure that the maximal velocity of the crowd is still normalized to the velocity model $f(\rho )$ despite the agents’ presence. Indeed, for the unregularized model,

$$\begin{aligned} f(\rho )^2|\nabla \phi | = f(\rho ) \end{aligned}$$

holds. When $\phi$ is replaced by $\phi + \phi _K$, this is no longer true. Thus we explicitly normalize the transport direction, i. e., instead of (1.2), we define the transport velocity by

$$\begin{aligned} \beta (\rho ,\phi ,{\varvec{x}}) :=f(\rho ) \, h ({\nabla (\phi + \phi _K({\varvec{x}};\cdot ))}) . \end{aligned}$$

(1.7)

Here, h is a smoothed projection onto the unit ball. Throughout this article we will use

$$\begin{aligned} h({\varvec{y}}) = \min \!{_\varepsilon }\{{1,|{\varvec{y}}|}\} \frac{{\varvec{y}}}{|{\varvec{y}}|} \end{aligned}$$

with $\min \!{_\varepsilon }$ a fixed smooth approximation of the minimum function.

Remark 1.2

(On the regularizations (II)) We emphasise that the regularisation of the function h serves no mathematical purpose (in fact, all of our analysis works without it) but is rather needed to preserve the crucial property of our model that the direction of motion is determined by a vector with unit length while the speed of motion is solely determined by the surrounding density via the function $f(\rho )$.

We assume that the agents move with maximal possible velocity towards a prescribed directions $u_i(t) \in \mathbb {R}^2$, $i=1,\ldots ,M$, which act as control in the system. Since the agents are also part of the crowd, their effective velocity will depend on the surrounding density in the same way as it does for all other individuals in the crowd. We therefore assume $|u_i(t)| \le 1$ and the law of motion for the agents will be

$$\begin{aligned} \dot{x}_i(t) = f({\overline{\rho }(t,x_i(t))})\,u_i(t),\quad i=1,\ldots , M, \end{aligned}$$

(1.8)

where $\overline{\rho }$ is a suitable extension of $\rho$ from $\varOmega$ to $\mathbb {R}^2$ which will be detailed later. While this extension is necessary to ensure the existence of solutions since the agents may leave the domain $\varOmega$ on which $\rho$ is defined, the precise choice of the extension is clearly irrelevant in terms of modeling. Notice also that the ODE (1.8) does not prevent agents from walking through walls or to attract people from behind a wall. The first issue can be avoided by imposing additional state constraints in an optimal control problem. The second issue can be avoided when assuming the wall thickness to be larger than the agents’ attraction radius.

The complete forward system which we are going to consider finally reads

$$\begin{aligned}&\partial _t\rho - \nabla \cdot ({\rho \,\beta (\rho ,\phi ,{\varvec{x}})}) = \varepsilon \, \varDelta \rho \quad \text {in } Q_T, \end{aligned}$$

(1.9a)

$$\begin{aligned}&-\delta _1\,\varDelta \phi + |\nabla \phi |^2 = \frac{1}{f(\rho )^2+\delta _2} \quad \text {in } Q_T, \end{aligned}$$

(1.9b)

$$\begin{aligned}&\dot{x}_i(t) = f({\overline{\rho }(t,x_i(t))})\, u_i(t) \quad \text {for } t\in (0,T), \quad i=1,\ldots ,M , \end{aligned}$$

(1.9c)

together with initial condition (1.4) on $\rho$ and

$$\begin{aligned} x_i(0) = x_{i,0} \quad \text { for } i = 1,\ldots , M, \end{aligned}$$

(1.10)

as well as boundary conditions (1.5).

The aim of the present paper is to investigate several optimal control problems for this coupled system. We seek an optimal control function ${\varvec{u}}$ such that the solution triple ${\varvec{y}}= (\rho ,\phi ,{\varvec{x}})$ is optimal in a certain sense. Depending on the application in mind it remains to define a suitable objective functional. We particularize the following two examples:

• Minimal evacuation time: In this case one seeks to minimize the time required for the evacuation of a room. The exits where individuals can leave the domain are located at the truncation boundary $\partial \varOmega _O$ which is sufficiently far away from the actual domain of interest. As time-optimal control problems with PDEs are rather challenging, see, e. g., [6, 40, 47, 52], we consider a simpler but closely related model. We fix a reasonably large final time $T>0$ and minimize

  $$\begin{aligned}&J(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) :=c_1 \int _\varOmega \rho (T,x) \mathop {}\!d x\nonumber \\&\quad + c_2 \int _0^T \int _\varOmega t \, \rho (t,x) \mathop {}\!d x\mathop {}\!d t+ \frac{\alpha }{2T} \sum _{i=1}^M \Vert u_i\Vert _{H^1(0,T)}^2 \end{aligned}$$

  (1.11)

  with weighting parameters $c_1,c_2>0$ and a regularization parameter $\alpha >0$. The first term in J penalizes individuals remaining in the room at time T. The second term encourages individuals to leave the room as early as possible. The last terms provides the required regularity for the control variables, so that the forward system (1.9) is well-defined. From the modeling point of view, these terms also avoid unrealistic trajectories of the agents.

• Optimal binding of a crowd: In some applications it might be desired to keep the group of individuals together, i. e., trying to maintain a single group during an evacuation. This is also motivated by a similar approach which has been used to model the repulsive interaction of dogs in a flock of sheeps, see [13]. To this end, we define the center of mass and variance of $\rho$ as

  $$\begin{aligned} E_\rho (t) :=\frac{1}{M(t)}\int _\varOmega x \rho (t,x) \mathop {}\!d x\quad \text {and} \quad V_\rho (t) :=\frac{1}{M(t)}\int _\varOmega \rho |x-E_\rho (t)|^2 \mathop {}\!d x, \end{aligned}$$

  with total mass $M(t) = \int _\varOmega \rho (t,x) \mathop {}\!d x$. A crowd is optimally kept together when the functional

  $$\begin{aligned} J(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) :=\frac{1}{2T} \int _0^T V_\rho (t) \mathop {}\!d t+ \frac{\alpha }{2T}\sum _{i=1}^M \Vert u_i\Vert _{H^1(0,T)}^2 \end{aligned}$$

  (1.12)

  is minimized.

Clearly, it is also possible to use a combination of the objective functionals (1.11) and (1.12).

1.2 Related Work

We briefly review the literature regarding the analysis of the original Hughes model and related optimal control problems.

A first contribution on existence of solutions for the Hughes model in the one-dimensional case is [26]. There it was shown, starting from the regularized version (1.3), that in the limit $\varepsilon \rightarrow 0$ a suitable unique entropy solution $\rho$ exists. The proof is based on a vanishing viscosity argument and Kruzkov’s doubling of variables technique to show uniqueness. The results were complemented by more detailed studies on the unregularized problem, also numerically. For instance generalizations to higher spatial dimensions can be found in [18], even for a slightly more general class of models. Further articles examine Riemann-type solutions to the unregularized problem; see [3, 4, 28] in one spatial dimension. As far as modeling is concerned, slightly different models were derived in [12] based on a mean field games approach. In [16], a modified approach using multiple local potentials $\phi _i$ instead of one global potential $\phi$ is introduced, removing the possibly unrealistic assumption that every pedestrian has complete information of the entire crowd. Moreover, in [15], a discrete pedestrian model in a graph network is studied. Further complex scenarios based on a regularized version of Hughes’ model are investigated in [19, 25].

From a broader perspective, the optimal control of (1.9) falls into the class of the optimal control of coupled ODE-PDE systems. Such problems, with models from a range of different applications have been analyzed, for instance, in [17, 36, 38, 39, 51]. We want to mention a few contributions dealing with optimal control of coupled PDE-ODE models in the context of pedestrian dynamics. In [30, 31] the feedback control in a multi-agent system is considered aiming at an optimization of the dynamics of a population. The individuals are either followers or leaders and the interaction is modeled by an ODE system. Moreover, the authors study also the limit case where the number of followers tends to infinity, which results in a couples system of ODEs and a PDE of Vlasov type. A similar approach is studied in [1]. Therein, the interaction between individuals is a short range retraction and a long range attraction. While leaders are not visible to the remaining crowd, they still influence it by taking part in these interactions. Then an optimal control problem arises as an external force acts on the leaders. The authors also consider, in the limit of many individuals (grazing interaction limit), macroscopic Boltzmann type equations for this interaction, while the number of leaders remains fixed and finite. Furthermore, in [13] external agents act as control. Again, they start from a microscopic ODE model and subsequently obtain a continuous model for the uncontrolled population by means of a mean field limit. For a general overview on interacting particle systems and control, we refer the reader to [46].

Closer to our approach is the work of [7]. There, a system of hyperbolic conservation laws for the density of different pedestrian groups is coupled to ODEs accounting for agents. Due to the low regularity of solutions to the hyperbolic equations, a regularization in the ODEs, similar to Lemma 3.6 in our case, is used. See also [8, 9] for a similar approach in different settings. Also related to our approach are the studies [20, 21] where evacuation is optimized by shape optimization of obstacles in a room. A more detailed survey on the control of several models for pedestrian dynamics can be found in [2]. Finally, we want to mention our proceeding article [45] where we develop a numerical solution approach for our model based on a structure preserving finite volume method and a projected gradient algorithm.

This paper is organized as follows. In Sect. 2 we collect the required notation, introduce some assumptions and state a precise existence and uniqueness result for the regularized system (1.9). The full forward system involving also the ordinary differential equation (1.8) is investigated in Sect. 3 and the linearized forward system in Sect. 4. The latter is required to establish the differentiability of the control-to-state map, which in turn is the basis of optimality conditions. Based on this we discuss the optimal control problem in Sect. 5 and derive first-order necessary optimality conditions. The presentation of numerical results will be postponed to a forthcoming publication.

2 Mathematical Preliminaries

Let us first state the assumptions on the domain and data.

• (A1) $\varOmega \subset \mathbb {R}^2$ is an open, bounded domain with $C^4$-boundary $\partial \varOmega$.

• (A2) There exist two measurable sets $\partial \varOmega _O , \partial \varOmega _W \subset \partial \varOmega$ s. t. $\partial \varOmega _O \cup \partial \varOmega _W = \partial \varOmega$ and $\partial \varOmega _O \cap \partial \varOmega _W = \emptyset$. Moreover, $\partial \varOmega _O$ has positive measure with respect to the Lebesgue measure on $\partial \varOmega$.

• (A3) The initial density satisfies $\rho _0\in W^{3/2,4}(\varOmega )$ and $0 \le \rho _0 \le 1$ a.e. in $\varOmega$.

• (A4) There holds $f \in W^{3,\infty }(\mathbb {R}) \cap C_c(\mathbb {R})$ with $f(0) = 1$, $f(1) = 0$ and $f(\rho )>0$ for all $\rho \in [0,1)$.

Moreover, we require some assumptions on the potential functions (1.6) of the agents, which depend on the kernel K.

• (K1) The kernel $K :\mathbb {R}^2 \rightarrow \mathbb {R}$ is radially symmetric, i. e., $K(x) = k(|x|)$, where $k :[0,\infty ) \rightarrow \mathbb {R}$ is nonnegative and decreasing.

• (K2) The kernel satisfies $K \in W^{3,\infty }(\mathbb {R}^2)$.

Finally, we consider an assumption on the velocity controls of the agents.

• (C1) There holds $u_i \in L^\infty (0,T;\mathbb {R}^2)$ and $\Vert u_i\Vert _{L^\infty (0,T)} \le 1$ for $i = 1, \ldots , M$.

Throughout this article we frequently exploit the boundedness and Lipschitz continuity of the functions K, f, h and $g(x) :=x\,f(x)$ and its derivatives. When doing so, we denote the bounds and Lipschitz constants by $C_K$, $C_f$, $C_h$, $C_g$ and $L_K$, $L_f$, $L_h$, $L_g$, respectively.

Remark 2.1

(Assumptions)

• (i) Our results extend to the case $d=3$ upon adopting the Sobolev embeddings used in several places.

• (ii) Assumption (A2) essentially means that the door and wall parts of the boundary “do not meet”, i. e., we do not consider a truly mixed boundary value problem in order to avoid technical conditions ensuring sufficient regularity of solutions. (A2) can be replaced, e. g., in two spatial dimensions, by suitable angle conditions on the points where the two parts of the boundary meet. The interested reader is referred to results in [33].

• (iii) The optimal regularity for the initial datum $\rho _0$ is the Besov space $B_{pp}^{2-2/p}(\varOmega )$; see [22].

• (iv) While in Assumption (A4) f is defined on all of $\mathbb {R}$, as far as the modeling is concerned, only ${f}{|}_{[0,1]}$ is relevant. Indeed, we will later see that the solution to (1.9a) satisfies $0\le \rho \le 1$.

• (v) A reasonable choice in (K1) and (K2) is the kernel function

  $$\begin{aligned} K(x-x_i(t)) = { {\left\{ \begin{array}{ll} s\exp \left( {-\frac{R^2}{R^2-|x-x_i(t)|^2}}\right) , &{}\text {if}\ |x-x_i(t)| < R,\\ 0, &{}\text {otherwise}, \end{array}\right. }} \end{aligned}$$

  (2.1)

  where $s > 0$ is an intensity factor and $R > 0$ is related to an attraction radius.

Notations. The space-time cylinder and its lateral surface are denoted by $Q_T :=(0,T) \times \varOmega$ and $\varSigma _T = (0,T) \times \partial \varOmega$, respectively. The boundary surface can be divided into

$$\begin{aligned} {\varSigma _{O ,T}} = (0,T) \times \partial \varOmega _O \text { and } \varSigma _{W ,T} = (0,T) \times \partial \varOmega _W . \end{aligned}$$

The (Frobenius) inner product of two matrices $A, B \in \mathbb {R}^{n \times n}$ is denoted by

$$\begin{aligned} A :B :={{\,\textrm{trace}\,}}(A^T B) = \sum _{i,j=1}^n a_{ij}b_{ij} . \end{aligned}$$

The Jacobian of a function $h :\mathbb {R}^2 \rightarrow \mathbb {R}^2$ is denoted by Dh and the Hessian of a function $u :\mathbb {R}^2 \rightarrow \mathbb {R}$ is denoted by $\nabla ^2u$. Furthermore, $\varphi _\gamma \in C^\infty _c(\mathbb {R}^2)$ is a standard mollifier, see, e. g., Ch. 4.4 in [11], i. e., a function satisfying

$$\begin{aligned} {{\,\textrm{supp}\,}}{\varphi _\gamma } \subset \overline{B_{\gamma }(0)} \quad \text {and} \quad \int _{\mathbb {R}^2} {\varphi _\gamma } \mathop {}\!d x= 1 . \end{aligned}$$

(2.2)

Note that

$$\begin{aligned} \Vert f - {\varphi _\gamma } *f\Vert _{C(\overline{\varOmega })} \rightarrow 0 \text { as } \gamma \rightarrow 0 \end{aligned}$$

(2.3)

for every continuous function $f \in C(\overline{\varOmega })$, cf. Prop. 4.21 in [11].

Finally, we introduce the subspace of $H^1(\varOmega )$ incorporating the Dirichlet boundary conditions as

$$\begin{aligned} H^1_D (\varOmega ) :={\{}{v \in H^1(\varOmega )}{|}{v = 0 \text { a.e. on } {\partial \varOmega _O }}{\}} \end{aligned}$$

and, for $p\in [1,\infty ]$, we denote the subspace of $W^{2,p}(\varOmega )$ fulfilling the boundary conditions (1.5) by

$$\begin{aligned} W^{2,p}_\text {ND}(\varOmega ) :={\{}{v\in W^{2,p}(\varOmega )}{|}{v = 0 \text { a.e. on } {\partial \varOmega _O }, \nabla v \cdot n = 0 \text { a.e. on } \partial \varOmega _W }{\}} . \end{aligned}$$

For time-dependent functions we introduce, for $p \in (1,\infty )$ and $r,s \in \mathbb {N}_0$, the spaces

$$\begin{aligned} W^{r,s}_p(Q_T) :=L^p(0,T;W^{r,p}(\varOmega )) \cap W^{s,p}(0,T;L^p(\varOmega )) , \end{aligned}$$

equipped with the natural norm $\left( {\Vert \cdot \Vert _{L^p(0,T;W^{r,p}(\varOmega ))}^p + \Vert \cdot \Vert _{W^{s,p}(0,T;L^p(\varOmega ))}^p}\right) ^{1/p}$. Spaces with non-integral r and s are defined, as usual, as (real) interpolation spaces. Of particular interest in our application is the space $W^{2,1}_p(Q_T)$ with $p > d = 2$, which fulfills the embedding

$$\begin{aligned} W^{2,1}_p(Q_T) \hookrightarrow C([0,T];W^{1,p}(\varOmega )) \hookrightarrow C(\overline{Q_T}) . \end{aligned}$$

(2.4)

This is needed in order to allow point evaluations of the density $\rho$, required in the ordinary differential equation (1.9c). Finally, for functions from the Hölder space $C^{1,\alpha }(\varOmega )$ we introduce the norm

$$\begin{aligned} \begin{aligned}&\Vert u\Vert _{C^{1,\alpha }(\varOmega )} :=\Vert u\Vert _{C^1(\varOmega )} + \max _{|\beta |=1} |D^{\beta } u|_{C^{0, \alpha }(\varOmega )} ,\\&\quad \text {where } |u|_{C^{0, \alpha }(\varOmega )} = \sup _{x \ne y \in \varOmega } \frac{|u(x)-u(y)|}{|x-y|^{\alpha }}. \end{aligned} \end{aligned}$$

In the following we collect some important properties of the function spaces used in this article.

Lemma 2.2

For each $\theta \in [0,1]$, $p \in (1,\infty )$ and $0\le s < r$, the continuous embedding

$$\begin{aligned} W^{1,p}(0,T;W^{s,p}(\varOmega ))\cap L^p(0,T;W^{r,p}(\varOmega )) \hookrightarrow W^{\theta ,p}(0,T; W^{\theta \,s + (1-\theta )\,r,p}(\varOmega )). \end{aligned}$$

holds.

Proof

See Lemma 4.3 in [24].

Similarly as above we define the Sobolev spaces

$$\begin{aligned} W^{r,s}_p(\varSigma _T) :=L^p(0,T;W^{r,p}(\partial \varOmega )) \cap W^{s,p}(0,T;L^p(\partial \varOmega )) \end{aligned}$$

on the lateral boundary $\varSigma _T = (0,T) \times \partial \varOmega$ of the space-time cylinder $Q_T$.

The following trace theorem is proved in [22, Lem. 3.5]; see also Sect. 2 in [23]:

Lemma 2.3

For $p>1$, the trace operators

$$\begin{aligned} \gamma _0&:W^{2,1}_p(Q_T) \rightarrow W^{2-1/p,1-1/(2p)}_p(\varSigma _T) , \\ \gamma _1&:W^{2,1}_p(Q_T) \rightarrow W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T) \end{aligned}$$

defined by $\gamma _0 \rho = {\rho }{|}_{\varSigma _T}$ and $\gamma _1 \rho = \nabla \rho \cdot n_{\varSigma _T}$ are bounded and have a continuous right inverse.

Recall the differential equation (1.8), where an extension to $\mathbb {R}^2$ of the density function $\rho$ is used. For theoretical purposes we will use an extension operator fulfilling the following result from Lemma 6.37 in [32]:

Lemma 2.4

Let $\alpha \in (0,1)$ be a fixed number. There exists a linear, continuous extension operator

$$\begin{aligned} {E }:C^{1,\alpha }(\overline{\varOmega }) \rightarrow C^{1,\alpha }(\mathbb {R}^2) \end{aligned}$$

such that

$$\begin{aligned} \Vert {E }f\Vert _{C^{1,\alpha }(\mathbb {R}^2)} \le C_{E }\, \Vert f\Vert _{C^{1,\alpha }(\varOmega )} \quad \text {and}\quad \Vert {E }f\Vert _{L^\infty (\mathbb {R}^2)} \le C_{{E },\infty }\Vert f\Vert _{L^\infty (\varOmega )}, \end{aligned}$$

holds for all $f \in C^{1,\alpha }(\overline{\varOmega })$. For brevity we will write $\overline{f} :={E }f$.

3 Analysis of the Forward System

This section is devoted to showing the existence of strong solutions to the forward system (1.9) with boundary and initial conditions (1.4), (1.5), (1.10). We proceed in two steps. First we provide auxiliary results on equation (1.9b) as well as on linear parabolic equations. Then we prove existence of solutions to the complete forward system.

3.1 Preliminary Results

First, we study the regularized Eikonal equation (1.9b).

Lemma 3.1

For given $\widetilde{\rho } :Q_T \rightarrow \mathbb {R}$ consider the equation

$$\begin{aligned} -\delta _1\,\varDelta \phi + |\nabla \phi |^2 = \frac{1}{f(\widetilde{\rho })^2+\delta _2} \quad \text {in } Q_T , \end{aligned}$$

with boundary conditions (1.5). We have:

• (i) If $\widetilde{\rho } \in C([0,T];L^2(\varOmega ))\cap H^1(0,T;H^1(\varOmega )^*)$, then there exists a unique strong solution which satisfies $\phi \in L^\infty (0,T;W^{2,p}_\text {ND}(\varOmega ))\cap H^1(0,T;H^1(\varOmega ))$ for all $2 \le p < \infty$. Moreover, the a priori estimates

  $$\begin{aligned} \Vert \phi \Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} \le \widetilde{C}_\phi \quad \text {and} \quad \Vert \phi \Vert _{H^1(0,T;H^1_D (\varOmega ))} \le \widetilde{C}_\phi \,\Vert \widetilde{\rho }\Vert _{H^1(0,T;H^1(\varOmega )^*)} \end{aligned}$$

  (3.1)

  hold with a positive constant $\widetilde{C}_\phi$ depending on $p,\delta _1,\,\delta _2,\,\varOmega$ and T only.

• (ii) If $\widetilde{\rho } \in W^{2,1}_p(Q_T)$ holds, the strong solution $\phi$ additionally belongs to $W^{4,1}_p(Q_T)$ and satisfies the a priori estimate

  $$\begin{aligned} \Vert \phi \Vert _{W^{1,p}(0,T;W^{2,p}(\varOmega ))} \le \overline{C}_\phi \,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} , \end{aligned}$$

  (3.2)

  with a constant $\overline{C}_\phi$ depending on $p,\, \delta _1,\,\delta _2,\,\varOmega$ and T only.

• (iii) For any $\widetilde{\rho }_1, \widetilde{\rho }_2 \in C([0,T];L^2(\varOmega ))$, the corresponding solutions $\phi _1$ and $\phi _2$ satisfy the Lipschitz estimate

  $$\begin{aligned} \Vert \phi _1(\cdot ,t) -\phi _2(\cdot ,t)\Vert _{W^{2,2}(\varOmega )} \le \widehat{C}_\phi \, \Vert \widetilde{\rho }_1(\cdot , t) - \widetilde{\rho }_2(\cdot ,t)\Vert _{L^2(\varOmega )} \end{aligned}$$

  (3.3)

  for all $t\in [0,T]$, with $\widehat{C}_\phi$ depending on $p,\, \delta _1,\, \delta _2,\,\varOmega$ and T only.

From now on we shall use the definition $C_\phi :=\max \{{\widetilde{C_\phi }, \overline{C}_\phi , \widehat{C}_\phi }\}$.

Proof

We first show assertion (i). Note that due to the continuity of $\widetilde{\rho }$ in time it makes sense to define, for fixed $t \in [0,T]$, the function

$$\begin{aligned} q_t(x) :=\frac{1}{\delta _1^2} \frac{1}{f(\widetilde{\rho }(t,x))^2 + \delta _2} . \end{aligned}$$

(3.4)

By assumption (A4) on f we then have $q_t \in L^\infty (\varOmega )$.

Step 1: Existence. First note that an application of the transformation

$$\begin{aligned} \psi (t,\cdot ) = \textrm{e}^{-\phi (t,\cdot )/\delta _1} - 1 \end{aligned}$$

(3.5)

to (1.9b) yields for all $t\in [0,T]$ the linear problem

$$\begin{aligned} \begin{aligned} -\varDelta \psi (t,\cdot ) + q_t\,\psi (t,\cdot )&= -q_t(\cdot ){} & {} \text {in } \varOmega , \\ \psi (t,\cdot )&= 0{} & {} \text {on }{\partial \varOmega _O } , \\ \partial _n \psi (t,\cdot )&= 0{} & {} \text {on } \partial \varOmega _W . \end{aligned} \end{aligned}$$

(3.6)

As $q_t \in L^\infty (\varOmega )$ and $q_t > 0$ a.e. in $\varOmega$, the Lax-Milgram lemma yields for every $t\in [0,T]$ the existence of a unique weak solution $\psi (t,\cdot ) \in H^1_D (\varOmega )$ and the estimate

$$\begin{aligned} \Vert \psi (t,\cdot )\Vert _{H^1(\varOmega )} \le C \, \Vert q_t\Vert _{L^2(\varOmega )} . \end{aligned}$$

Elliptic regularity theory (see for instance Theorem 3.17 in [49]) then implies $\psi (t,\cdot ) \in W^{2,q}(\varOmega )$ for any $2 \le q < \infty$, and there holds

$$\begin{aligned} \Vert \psi (t,\cdot )\Vert _{W^{2,q}(\varOmega )} \le C \, \Vert q_t\Vert _{L^q(\varOmega )} =:C_\phi \end{aligned}$$

(3.7)

for $t\in [0,T]$. For given $x \in \varOmega$, the function $t \mapsto q_t(x)$ belongs to $L^\infty (0,T)$ again by (A4). Taking the supremum over $0 \le t \le T$ yields the regularity $\psi \in L^\infty (0,T;W^{2,q}(\varOmega ))$. Formally differentiating equation (3.6) with respect to time we see that $\partial _t \psi$ satisfies

$$\begin{aligned} -\varDelta \partial _t\psi (t,\cdot ) + q_t\partial _t\psi (t,\cdot ) = -(1+\psi (t,\cdot ))(\partial _t q_t) \quad \text {in } \varOmega . \end{aligned}$$

Using its definition and the fact that $\widetilde{\rho } \in H^1(0,T;H^1(\varOmega )^*)$, the derivative of $q_t(x)$ with respect to time is an element of $L^2(0,T;H^1(\varOmega )^*)$. Therefore, Sect. 2.2.2, Corollary, p. 99 in [49] (after replacing $\partial _t q_t(x)$ by a continuous in time approximation and passing to the limit) yields the existence of $\partial _t \psi \in L^2(0,T;H^1(\varOmega ))$ with

$$\begin{aligned} \Vert \partial _t \psi \Vert _{L^2(0,T;H^1(\varOmega ))} \le C \, \Vert \partial _t \widetilde{\rho }\Vert _{L^2(0,T;H^1(\varOmega )^*)} . \end{aligned}$$

(3.8)

Step 2: Strict lower bound. In order to invert (3.5) and obtain a solution to (1.9b), we need to ensure that $\psi > -1$ holds. This follows from Theorem 4 in [43], provided that we can show $\psi \ge -1$, $\Vert \psi \Vert _{L^\infty (Q_T)} < \infty$ and Hölder continuity of $\psi (t,\cdot )$ for a.a. $t\in (0,T)$. The last two assertions are a direct consequence of the embeddings $W^{2,p}(\varOmega ) \hookrightarrow L^\infty (\varOmega )$ and $W^{2,p}(\varOmega ) \hookrightarrow C^{0,1/2}(\varOmega )$ as $p \ge 2$ and we are in dimension 2, combined with (3.7). To show the first part we choose $\phi = (\psi +1)_-$, i. e., the negative part of $\psi +1$, as test function in the weak formulation of (3.6) and obtain

$$\begin{aligned} \int _\varOmega |\nabla (\psi (t,\cdot )+1)_-|^2 \mathop {}\!d x+ \int _\varOmega q_t(x)(\psi + 1)_-^2 \mathop {}\!d x= 0 . \end{aligned}$$

Since $q_t$ is strictly positive, this implies that $(\psi + 1)_- = 0$ holds a.e. in $\varOmega$. Thus there exists a positive constant $C_\psi$ s. t.

$$\begin{aligned} \psi \ge C_\psi > - 1 \text { a.e. in } \varOmega . \end{aligned}$$

We can therefore invert the transformation (3.5) and conclude that the solution of (1.9b), (1.5) fulfills the desired regularity, as the regularity is unaffected by the transformation (using $\psi \in L^\infty (0,T;L^\infty (\varOmega ))$ and $\psi > -1$). This ends the proof of assertion (i). To conclude (ii) we only have to show the additional regularity $\partial _t \psi \in L^p(0,T;W^{2,p}(\varOmega ))$ which follows under the assumption $\widetilde{\rho } \in W^{2,1}_p(Q_T)$ and standard elliptic regularity theory, see [49, Theorem 3.17].

Similarly, in case $\widetilde{\rho } \in W^{2,1}_p(Q_T)$ we obtain $\partial _t \psi \in L^p(0,T;W^{2,p}(\varOmega ))$, which gives the desired estimate.

Step 3: Lipschitz estimate. To show (iii), denote by $q_{t,1}$, $q_{t,2}$ the respective coefficients for $\widetilde{\rho }_1$ and $\widetilde{\rho }_2$ as in (3.4) and, analogously, let $\psi _1$ and $\psi _2$ be the respective solutions to (3.6). Then $\overline{\psi }= \psi _1 - \psi _2$ satisfies

$$\begin{aligned} \begin{aligned} -\varDelta \overline{\psi }(t,\cdot ) + q_{t,1}(x) \, \overline{\psi }(t,\cdot )&= (1 -\psi _2(t,\cdot )) \, \overline{q}_t(x){} & {} \text {in } \varOmega , \\ \overline{\psi }(t,\cdot )&= 0{} & {} \text {on }{\partial \varOmega _O } , \\ \partial _n \overline{\psi }(t,\cdot )&= 0{} & {} \text {on } \partial \varOmega _W , \end{aligned} \end{aligned}$$

(3.9)

for $t\in [0,T]$ with $\overline{q}_t = q_{t,1} - q_{t,2}$. Noting that $q_t$ is Lipschitz continuous as a function of $\widetilde{\rho }$ (due to (A4)) and applying the a priori estimate (3.7) to (3.9) yields

$$\begin{aligned} \Vert \overline{\psi }(t,\cdot )\Vert _{W^{2,2}(\varOmega )} \le C \, \Vert \overline{q}_t\Vert _{L^2(\varOmega )} \le C \, \Vert \widetilde{\rho }_1(t,\cdot ) - \widetilde{\rho }_2(t,\cdot )\Vert _{L^2(\varOmega )} , \end{aligned}$$

where we used the boundedness of $(1-\psi _2(t,\cdot ))$ in $L^\infty (\varOmega )$ and, again, $W^{2,p}(\varOmega ) \hookrightarrow L^\infty (\varOmega )$ and (3.7). This implies (iii) and completes the proof.

To obtain the desired $W^{2,1}_p(Q_T)$-regularity of the density function $\rho$ we will need the following lemma taken from [22] but adopted to our notation.

Lemma 3.2

Let assumptions (A1)–(A4) hold, $p\in [2,\infty )$ and $\varepsilon > 0$. Suppose that $c \in L^p(Q_T)$, $b \in L^p(0,T;L^\infty (\varOmega ))$, $r \in W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T)$ and $\rho _0\in W^{2-2/p,p}(\varOmega )$ are given. Then the problem

$$\begin{aligned} \begin{aligned} \partial _t \rho - \varepsilon \varDelta \rho + b\cdot \nabla \rho&= c{} & {} \text {in } Q_T , \\ \varepsilon \nabla \rho \cdot n&= r{} & {} \text {on }\varSigma _{W ,T} , \\ \rho&= 0{} & {} {\text {on }\varSigma _{O ,T}} , \\ \rho (0)&= \rho _0{} & {} \text {in } \varOmega \end{aligned} \end{aligned}$$

admits a unique strong solution $\rho \in W^{2,1}_p(Q_T)$ depending continuously on the input data c, r and $\rho _0$.

Proof

See Theorem 2.1 in [22].

Finally, we need the following regularity result for $p=2$ with flux boundary conditions.

Lemma 3.3

Given $h \in H^1(0,T;H^1(\varOmega )) \cap L^\infty (Q_T)$, $g \in L^\infty (\mathbb {R})$ and $\rho _0 \in H^1(\varOmega )$, the variational problem

$$\begin{aligned}{} & {} \int _{Q_T} \partial _t \rho \,\xi \mathop {}\!d x\mathop {}\!d t+ \varepsilon \int _{Q_T} \nabla \rho \cdot \nabla \xi \mathop {}\!d x\mathop {}\!d t- \int _{Q_T} g(\rho )\,h\cdot \nabla \xi \mathop {}\!d x\mathop {}\!d t\nonumber \\{} & {} \quad = -\eta \int _{\varSigma _T} {\chi _{\partial \varOmega _O }} \, \rho \, \xi \mathop {}\!d s_x \mathop {}\!d t\quad \text {for all } \xi \in L^2(0,T;H^1(\varOmega )) \end{aligned}$$

(3.10)

has a unique solution $\rho \in L^\infty (0,T;H^1(\varOmega ))\cap H^1(0,T;L^2(\varOmega ))$ with $\rho (0) = \rho _0$.

The proof mainly uses standard methods but since, to the best of the authors’ knowledge, a proof matching our boundary conditions is not available in the literature, we included it into Sect. A.

Next we define our notion of solution for the ODE (1.9c).

Definition 3.4

Fix $2< p < \infty$ and $\rho \in W^{2,1}_p(Q_T)$. Then for given f and u satisfying assumptions (A4) and (C1) and an initial value $x_0\in \mathbb {R}^2$ we say that $x :[0,T] \rightarrow \mathbb {R}^2$ is a solution to

$$\begin{aligned} \dot{x}(t) = f({\overline{\rho }(x(t),t)}) \, u(t), \end{aligned}$$

(3.11)

if it is absolutely continuous, satisfies (3.11) for a.a. $t \in [0,T]$ and $x(0) = x_0$.

We have the following result about the existence of a solution of (3.11).

Lemma 3.5

For given $2<p < 4$, $u\in L^\infty (0,T;\mathbb {R}^2)$ satisfying assumption (C1) and $\rho \in W^{2,1}_p(Q_T)$, there exists a unique, absolutely continuous solution $x :[0,T] \rightarrow \mathbb {R}^2$ to (3.11) satisfying $x(0) = x_0$. Furthermore, $x \in W^{1,\infty }(0,T)$ holds. For $\rho _1,\,\rho _2 \in W^{2,1}_p(Q_T)$, the corresponding solutions $x_1$ and $x_2$ satisfy

$$\begin{aligned} \Vert x_1-x_2\Vert _{L^\infty (0,T;\mathbb {R}^2)} \le C_s\,\Vert \rho _1 - \rho _2\Vert _{L^\infty (Q_T)} , \end{aligned}$$

(3.12)

where the constant $C_s$ depends on u, T, $C_{E }$, $C_{{E },\infty }$ and the Lipschitz constants of f and $\rho$.

Proof

First note that $\rho \in L^p(0,T;W^{2,p}(\varOmega )) \hookrightarrow L^p(0,T;C^{1,\alpha }(\varOmega ))$, for some $\alpha > 0$, so that the application of the extension operator from Lemma 2.4 is well-defined. Since also f is Lipschitz continuous with Lipschitz constant $L_f$ by assumption (A4), the function $f(\overline{\rho }(t,x))$ satisfies the Carathéodory conditions (see Definition A.1 in the Appendix) and thus there exists a solution in the sense of Definition 3.4; see Ch. I, Theorem 5.1 in [34]. Furthermore, as we also have $\overline{\rho }\in L^p(0,T;C^{0,1}(\mathbb {R}^2))$, we obtain with $C_u:=\Vert u\Vert _{L^\infty (0,T)}$, $C_\rho (t) = \Vert \rho (t)\Vert _{W^{2,p}(\varOmega )}$, and the property $C_\rho \in L^p(0,T)$ that

$$\begin{aligned}&|f(\overline{\rho }(t,x)) - f(\overline{\rho }(t,y))| \le L_f\,|\overline{\rho }(t,x)-\overline{\rho }(t,y)| \le L_f\,\Vert \overline{\rho }(t)\Vert _{C^{1,\alpha }(\mathbb {R}^2)}\,|x-y| \nonumber \\&\quad \le L_f\,C_E\,\Vert \rho (t)\Vert _{C^{1,\alpha }(\varOmega )}\,|x-y| \nonumber \\&\quad \le L_{f,{E },\rho }(t)\,|x-y| \quad \text { for} \quad \hbox {a.a.} \quad t \in [0,T] \text { and all } x, y \in \mathbb {R}^2, \end{aligned}$$

(3.13)

where $L_{f,{E },\rho }(t) :=L_f\,C_E\,C_\infty \,C_\rho (t)$ and with $C_\infty$ the embedding constant for $W^{2,p}(\varOmega ) \hookrightarrow C^{1,\alpha }(\overline{\varOmega })$. The estimate (3.13) implies uniqueness by Ch. I, Thm. 5.3 in [34]. The additional regularity $x \in W^{1,\infty }(0,T;\mathbb {R}^2)$ is a consequence of the boundedness in $L^\infty (0,T;\mathbb {R}^2)$ of the right-hand side of (3.11). To establish the stability estimate we show

$$\begin{aligned}&|x_1(t) - x_2(t)| \le \int _0^t |f(\overline{\rho }_1(s,x_1(s))) - f(\overline{\rho }_2(s,x_2(s)))|\,|u(s)|\mathop {}\!d s\\&\quad \le \int _0^t |f(\overline{\rho }_1(s,x_1(s))) - f(\overline{\rho }_1(s,x_2(s)))|\,{|u(s)|} \mathop {}\!d s\\&\qquad + \int _0^t |f(\overline{\rho }_1(s,x_2(s))) - f(\overline{\rho }_2(s,x_2(s)))|\,{|u(s)|} \mathop {}\!d s\\&\quad \le \int _0^t {C_u}\,L_{f,\rho ,{E }}(s)\,|x_1(s) - x_2(s)|\mathop {}\!d s+ L_f \, {C_u}\,C_{E }\, C_{{E },\infty }\, \,t\,\Vert \rho _1 - \rho _2\Vert _{L^\infty (Q_t)}, \end{aligned}$$

where we used that the extension ${E }$ is also continuous with respect to the $L^\infty$-norm. An application of Gronwall’s inequality in integral form then yields, for $t \in (0,T)$,

$$\begin{aligned} |x_1(t) - x_2(t)|&\le L_f \, C_u\,C_{E }\, C_{{E },\infty }\,t\,\Vert \rho _1 - \rho _2\Vert _{L^\infty (Q_T)} \, \exp \left( {\int _0^t C_u\,L_{f,\rho ,{E }}(r) \mathop {}\!d r}\right) \\&\le C_s ({u, L_f, L_{f,\rho ,{E }}, C_{E }, C_{{E },\infty }, T}) \, \Vert \rho _1 - \rho _2\Vert _{L^\infty (Q_T)} . \end{aligned}$$

Next, we state an existence and stability result for a regularized version of (3.11). Note that the following result requires less regularity for the density function $\rho$.

Lemma 3.6

Fix $2<p < \infty$ and $\rho \in C([0,T];L^2(\varOmega ))$. Then for given $u\in L^\infty (0,T)$, f satisfying assumption (A4) and (C1), and every $\gamma > 0$, there exists, a unique, absolutely continuous solution $x :[0,T]\rightarrow \mathbb {R}^2$ to

$$\begin{aligned} \dot{x}(t) = f({(\overline{{\varphi _\gamma } *\rho })(t,x(t))}) \, u(t) \end{aligned}$$

(3.14)

satisfying $x(0) = x_0 \in \mathbb {R}^2$ and $x \in W^{1,\infty }(0,T)$. Here, ${\varphi _\gamma }$ is a standard mollifier as in (2.2) and $*$ denotes the convolution w.r.t. to the x-variable. Furthermore, for $\rho _1,\,\rho _2 \in C([0,T];L^2(\varOmega ))$, the corresponding solutions $x_1$ and $x_2$ satisfy

$$\begin{aligned} |x_1(t)-x_2(t)| \le {C_1\,T\,e^{C_2\,t}} \,\Vert \rho _1 - \rho _2\Vert _{L^\infty (0,T;L^2(\varOmega ))} , \end{aligned}$$

(3.15)

where the constants $C_1, C_2$ depend on u, $\gamma$, $C_{E }$, $C_{{E },\infty }$ and the Lipschitz constants of f and $\rho _1$, $\rho _2$.

Proof

As f and $\overline{{\varphi _\gamma } *\rho }$ are Lipschitz continuous (see assumption (A4) and Lemma A.2), the function $f({(\overline{{\varphi _\gamma } *\rho })(t,x)})$ satisfies the Carathéodory conditions of Definition 3.4 and thus there exists a solution in the sense of Definition 3.4; see Ch. I, Thm. 5.1 in [34]. In particular there exists a positive function $L_{f,\gamma ,{E }} >0$ such that

$$\begin{aligned} |f((\overline{{\varphi _\gamma } *\rho })(t,x)) - f((\overline{{\varphi _\gamma } *\rho })(t,y))| \le {L_{f,\gamma ,{E }}\,\Vert \rho (t,\cdot )\Vert _{L^1(\varOmega )}}\,|x-y| \end{aligned}$$

(3.16)

for a.a. $t \in [0,T]$ and all $x, y \in \varOmega$. This implies uniqueness by Ch. I, Thm. 5.3 in [34].

To show the stability estimate (3.15) we proceed as in the proof of Lemma Lemma 3.5, employ the estimate (3.16), Young’s inequality for convolutions $\Vert \varphi _\gamma *\rho \Vert _{L^\infty (\varOmega )} \le \Vert \varphi _\gamma \Vert _{L^2(\varOmega )}\,\Vert \rho \Vert _{L^2(\varOmega )}$, and obtain

$$\begin{aligned} |x_1(t)-x_2(t)|&\le \int _0^t |f((\overline{\varphi _\gamma *\rho _1})(s,x_1(s))) - f((\overline{\varphi _\gamma *\rho _2})(s,x_1(s)))||u(s)|\mathop {}\!d s\\&\quad +\int _0^t |f(\overline{(\varphi _\gamma *\rho _2)}(s,x_1(s))) - f((\overline{\varphi _\gamma *\rho _2})(s,x_2(s)))||u(s)|\mathop {}\!d s\\&\le L_f\,C_{u, {E }}\,\Vert \varphi _\gamma \Vert _{L^2(\varOmega )}\,t\,\Vert \rho _1-\rho _2\Vert _{L^\infty (0,T;L^2(\varOmega ))} \\&\quad + L_{f,\gamma ,{E }}\,C_u\,\Vert \rho _2\Vert _{L^\infty (0,T;L^2(\varOmega ))} \int _0^t|x_1(s)-x_2(s)|\mathop {}\!d s. \end{aligned}$$

Then, (3.15) follows from the Gronwall Lemma.

We continue with the following Lipschitz estimate for the transport term in (1.9a), which is needed in the theorem right after the next.

Lemma 3.7

Given $\rho _1,\,\rho _2 \in C([0,T];H^1(\varOmega ))$, $\phi _1,\,\phi _2 \in L^\infty (0,T;W^{2,p}(\varOmega ))$ and ${\varvec{x}}_1,\,{\varvec{x}}_2 \in L^\infty (0,T;\mathbb {R}^2)$ the function $\beta (\rho ,\phi ,{\varvec{x}})$ satisfies the Lipschitz inequality

$$\begin{aligned}&\Vert \rho _1\,\beta (\rho _1,\phi _1,{\varvec{x}}_1) - \rho _2\,\beta (\rho _2,\phi _2,{\varvec{x}}_2)\Vert _{L^2(Q_T)} \\&\quad \le L_\beta \left( { \Vert \rho _1-\rho _2\Vert _{L^2(Q_T)} + \Vert \nabla \phi _1 - \nabla \phi _2\Vert _{L^2(Q_T)} + \Vert {\varvec{x}}_1 - {\varvec{x}}_2\Vert _{L^\infty (0,T;\mathbb {R}^2)^M}}\right) , \end{aligned}$$

where $L_\beta$ depends on the Lipschitz and boundedness constants of f, h and K.

Proof

We define $g(\rho ) = \rho \,f(\rho )$ and write $\rho _i\,\beta (\rho _i,\phi _i,{\varvec{x}}_i) = g(\rho _i) \, h(\nabla (\phi _i+\phi _K({\varvec{x}}_i;\cdot )))$, $i=1,2$. Using the Lipschitz continuity of g, h and $\phi _K$ we obtain

$$\begin{aligned}&\Vert \rho _1\,\beta (\rho _1,\phi _1,{\varvec{x}}_1) - \rho _2\,\beta (\rho _2,\phi _2,{\varvec{x}}_2)\Vert _{L^2(0,T;L^2(\varOmega ))} \\&\quad \le L_g\,\Vert \rho _1-\rho _2\Vert _{L^2(Q_T)}\,\Vert h ({\nabla (\phi _1+\phi _K({\varvec{x}}_1;\cdot ))})\Vert _{L^\infty (Q_T)} \\&\qquad + C_g \, L_h ({\Vert \nabla (\phi _1-\phi _2)\Vert _{L^2(Q_T)} + L_{\phi _K} \Vert {\varvec{x}}_1-{\varvec{x}}_2\Vert _{L^\infty (0,T;\mathbb {R}^2)^M}}). \end{aligned}$$

3.2 Existence for the Full Forward System

We are now in a position to show the following existence and uniqueness result.

Theorem 3.8

Let assumptions (A1), (A2), (A3) (A4), (K1), (K2) (C1) hold and fix $2<p < \infty$. Then for any given control ${\varvec{u}}= (u_1,\ldots , u_M)^T \in L^\infty (0,T;\mathbb {R}^2)^M$ and any $T>0$, there exists a unique solution $(\rho ,\phi ,{\varvec{x}})$ to (1.9) with initial and boundary conditions (1.4), (1.5), (1.10) which satisfies $\rho \in W^{2,1}_p(Q_T)$, $\phi \in L^\infty (0,T;W^{2,p}(\varOmega ))$ and ${\varvec{x}}$ is a solution to (1.9c) in the sense of Definition 3.4. Moreover, the a priori estimate

$$\begin{aligned} \Vert \rho \Vert _{W^{2,1}_p(Q_T)} + \Vert \phi \Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} \le C \Vert \rho _0\Vert _{W^{1,p}(\varOmega )}, \end{aligned}$$

holds with C depending only on the domain, the bounds for the coefficients and the respective kernel.

The structure of (1.9a)–(1.9b) is very similar to chemotaxis models with volume filling, see for instance [44], except for the additional nonlinearity of the squared gradient term in (1.9b), which can, however, be handled using Lemma 3.1. Therefore, the existence and uniqueness of solutions can be proved using Banach’s fixed point theorem, similar to, e. g., Thm. 3.1 in [27]. The main issue in our situation is the additional coupling to the system of ODEs (1.9c), which requires $\rho$ to be regular enough to allow point evaluations. Our strategy is to introduce an additional regularization in (1.9c) in order to be able to perform the fixed point argument in the relatively “large” space $C([0,T];L^2(\varOmega ))$. We then show additional regularity and pass to the limit to recover the original system.

Proof

The proof consists of two parts. First we show existence with (1.9c) replaced by the regularized version (3.14). Then, we perform the limit $\gamma \rightarrow 0$ to recover the original problem.

Step 1: Fixed point argument: We consider the fixed-point operator

$$\begin{aligned} S :C([0,T];L^2(\varOmega )) \rightarrow C([0,T];L^2(\varOmega )) , \quad \widetilde{\rho } \mapsto \rho _\gamma , \end{aligned}$$

(3.17)

where, for fixed $\gamma > 0$, $\rho _\gamma$ is the unique weak solution to the system

$$\begin{aligned}&\partial _t\rho _\gamma - \nabla \cdot ({\widetilde{\rho }\,\beta (\widetilde{\rho }, \phi _\gamma , {\varvec{x}}_\gamma )}) = \varepsilon \,\varDelta \rho _\gamma \quad \text {in } Q_T , \end{aligned}$$

(3.18)

$$\begin{aligned}&- \delta _1\,\varDelta \phi _\gamma + |\nabla \phi _\gamma |^2 = \frac{1}{f(\widetilde{\rho })^2+\delta _2} \quad \text {in } Q_T , \end{aligned}$$

(3.19)

$$\begin{aligned}&\dot{x}_{i,\gamma }(t) = f({(\overline{{\varphi _\gamma } *\widetilde{\rho }})(t,x_{i,\gamma }(t))})\,u_i(t) \quad \text {for } t \in [0,T], \quad i=1,\ldots , M \end{aligned}$$

(3.20)

subject to the boundary conditions

$$\begin{aligned} \varepsilon \,\nabla \rho _\gamma \cdot n = - \widetilde{\rho }\,\beta (\widetilde{\rho },\phi _\gamma ,{\varvec{x}}_\gamma )\cdot n - \eta \,\rho _\gamma \,{\chi _{\partial \varOmega _O }} \quad \text {on } \varSigma _T . \end{aligned}$$

The boundary conditions for $\phi$ and the initial conditions for all variables are as in (1.5) and (1.4), (1.10), respectively. The equation for $\rho _\gamma$ is understood in the weak sense, i. e.,

$$\begin{aligned} \int _\varOmega \partial _t \rho _\gamma \,\xi \mathop {}\!d x+ \varepsilon \int _\varOmega \nabla \rho _\gamma \cdot \nabla \xi \mathop {}\!d x= - \int _\varOmega \widetilde{\rho }\,\beta (\widetilde{\rho },\phi _\gamma ,{\varvec{x}}_\gamma )\cdot \nabla \xi \mathop {}\!d x- \eta \int _{{\partial \varOmega _O }} \rho _\gamma \, \xi \mathop {}\!d s_x \end{aligned}$$

(3.21)

for all $\xi \in H^1(\varOmega )$ and a.a. $t\in [0,T]$. Applying Lemma 3.3 with $g(\rho ) = \rho f(\rho )$ and $h = h({\nabla (\phi _\gamma + \phi _K({\varvec{x}}_\gamma ;\cdot ))})$ shows that there exists a unique weak solution

$$\begin{aligned} \rho _\gamma \in L^\infty (0,T;H^1(\varOmega ))\cap H^1(0,T;L^2(\varOmega )) \hookrightarrow C([0,T];L^2(\varOmega )) . \end{aligned}$$

Choosing $\rho$ itself as a test function in the weak formulation (3.21) yields, after an application of the weighted Young’s inequality with parameter $\varepsilon / 2$, the trace inequality for $H^1$-functions, and an integration in time, the a priori estimate

$$\begin{aligned} \Vert \rho _\gamma \Vert _{C([0,T];L^2(\varOmega ))}&\le C ({\Vert \widetilde{\rho }\,\beta (\widetilde{\rho },\phi _\gamma ,{\varvec{x}}_\gamma )\Vert _{L^2(0,T; L^2(\varOmega ))} + \Vert \rho _0\Vert _{L^2(\varOmega )}}) \\&\le C ({\sqrt{T\,|\varOmega |} + \Vert \rho _0\Vert _{L^2(\varOmega )}}) , \end{aligned}$$

where we used $|\tilde{\rho }\,\beta (\tilde{\rho },\phi _\gamma , {\varvec{x}}_\gamma )| \le 1$.

Next, we show that S is a contraction. Let $\widetilde{\rho }_1, \, \widetilde{\rho }_2 \in C([0,T];L^2(\varOmega ))$ be arbitrary. The corresponding densities, potentials and agent positions are denoted by $\rho _{\gamma ,1}$, $\rho _{\gamma ,2}$, $\phi _{\gamma ,1}$, $\phi _{\gamma ,2}$ and ${\varvec{x}}_{\gamma ,1}$, ${\varvec{x}}_{\gamma ,2}$, respectively. Using Lemma 3.7 we have

$$\begin{aligned}&\Vert \rho _{\gamma ,1}-\rho _{\gamma ,2}\Vert _{C([0,T];L^2(\varOmega ))} \nonumber \\&\quad \le \Vert \widetilde{\rho }_1\,\beta (\widetilde{\rho }_1,\phi _{\gamma ,1},{\varvec{x}}_{\gamma ,1}) - \widetilde{\rho }_2\,\beta (\widetilde{\rho }_2,\phi _{\gamma ,2},{\varvec{x}}_{\gamma ,2})\Vert _{L^2(0,T;L^2(\varOmega ))} \nonumber \\&\quad \le L_\beta ({\Vert \widetilde{\rho }_1-\widetilde{\rho }_2\Vert _{L^2(0,T;L^2(\varOmega ))} + \Vert \nabla (\phi _{\gamma ,1} - \phi _{\gamma ,2})\Vert _{L^2(0,T;L^2(\varOmega ))}} \nonumber \\&\quad + {\Vert {\varvec{x}}_{\gamma ,1} - {\varvec{x}}_{\gamma ,2}\Vert _{L^\infty (0,T)^M}}). \end{aligned}$$

(3.22)

Applying Lemmas 3.1 and 3.6 to the second and third terms on the right-hand side, respectively, yields

$$\begin{aligned} \Vert \rho _{\gamma ,1}-\rho _{\gamma ,2}\Vert _{C([0,T];L^2(\varOmega ))} \le {C\,\max \{\sqrt{T}, T\,e^{C_2\,T}\}}\, \Vert \widetilde{\rho }_1 - \widetilde{\rho }_2\Vert _{C([0,T];L^2(\varOmega ))} . \end{aligned}$$

Thus, we can again find T small enough so that S is a contraction and Banach’s fixed point theorem asserts the existence of a unique solution.

The box constraints $0 \le \rho _\gamma \le 1$ a.e. in $Q_T$ follow by applying Lem. 3 in [27], i. e., by testing with smoothed versions of the positive part of $\rho -1$ and $-\rho$, respectively. In view of these uniform estimates a standard continuation argument yields existence for arbitrary $T>0$.

So far we have shown

$$\begin{aligned} \rho _\gamma \in L^2(0,T;H^1(\varOmega ))\cap H^1(0,T;H^1(\varOmega )^*) \cap L^\infty (Q_T) , \end{aligned}$$

since the fixed-point satisfies the weak formulation (3.21). Note in particular that $\rho _\gamma$ is bounded in $L^2(0,T;H^1(\varOmega ))\cap H^1(0,T;H^1(\varOmega )^*) \cap L^\infty (Q_T)$ by a constant independent of $\gamma$. Furthermore, we denote the corresponding potential and agent trajectories by

$$\begin{aligned} \phi _\gamma \in L^\infty (0,T;W^{2,p}(\varOmega ))\cap H^1(0,T;H^1(\varOmega )) , \quad {\varvec{x}}_\gamma \in W^{1,\infty }(0,T)^M . \end{aligned}$$

Step 2: Additional regularity: To shorten the notation we write the nonlinear term in the form

$$\begin{aligned} \rho _\gamma \,\beta (\rho _\gamma ,\phi _\gamma ,{\varvec{x}}_\gamma ) = g(\rho _\gamma )\,h(\varPhi _\gamma ) \end{aligned}$$

with $g(\rho _\gamma ) :=\rho _\gamma \,f(\rho _\gamma )$ and $\varPhi _\gamma :=\nabla (\phi _\gamma +\phi _K({\varvec{x}}_\gamma ;\cdot ))$. From the product rule we obtain the following representation for the divergence,

$$\begin{aligned} \nabla \cdot ({\rho _\gamma \,\beta (\rho _\gamma ,\phi _\gamma ,{\varvec{x}}_\gamma )}) = g'(\rho _\gamma )\,\nabla \rho _\gamma \cdot h(\varPhi _\gamma ) + \rho _\gamma \,f(\rho _\gamma ) \, \nabla \cdot h(\varPhi _\gamma ) . \end{aligned}$$

Freezing the nonlinear terms allows us to understand (3.18) as a linear equation of the form

$$\begin{aligned} \begin{aligned} \partial _t \rho _\gamma - \varepsilon \,\varDelta \rho _\gamma + b(t,x)\cdot \nabla \rho _\gamma + c(t,x)\,\rho _\gamma&= 0{} & {} \text {in } Q_T , \\ \varepsilon \nabla \rho _\gamma \cdot n&= r{} & {} \text {on } \varSigma _T , \\ \rho (0)&= \rho _0{} & {} \text {in } \varOmega \end{aligned} \end{aligned}$$

(3.23)

with

$$\begin{aligned} \begin{aligned} b(t,x)&:=g'(\rho _\gamma )\,h(\varPhi _\gamma ) \in L^p(0,T;L^\infty (\varOmega )) , \\ c(t,x)&:=f(\rho _\gamma )\,\nabla \cdot h(\varPhi _\gamma ) \in L^p(Q_T) , \end{aligned} \end{aligned}$$

(3.24)

and

$$\begin{aligned} r(t,x) :=- g(\rho _\gamma ) \, h(\varPhi _\gamma ) \cdot n - \eta \, \rho _\gamma \, {\chi _{\partial \varOmega _O } } \in W^{\kappa ,\kappa /2}_p(\varSigma _T) \end{aligned}$$

(3.25)

with $\kappa = 1-1/p$. It remains to show the regularity claimed for b, c and r. Again, this follows from the Hölder inequality and the regularity already shown for $\rho _\gamma$, $\phi _\gamma$, ${\varvec{x}}_\gamma$. Together with the product and the chain rule this leads to

$$\begin{aligned} \Vert b\Vert _{L^\infty (Q_T)}&\le \Vert g'\Vert _{L^\infty (\mathbb {R})}\,\Vert h\Vert _{L^\infty (\mathbb {R}^2)}\le C_b , \end{aligned}$$

(3.26)

$$\begin{aligned} \Vert c\Vert _{L^p(Q_T)}&\le \Vert f\Vert _{W^{1,\infty }(\mathbb {R})}\,\Vert Dh\Vert _{W^{1,\infty }(\mathbb {R}^2)} ({\Vert \phi _\gamma \Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} + C_{\phi _K}}) \le C_c . \end{aligned}$$

(3.27)

Note that $C_b$ and $C_c$ are independent of $\gamma$ as, in particular, $\phi _\gamma$ can be bounded independently of $\rho _\gamma$ and thus of $\gamma$. To show the required regularity for r we proceed as follows. First, we show the estimates

$$\begin{aligned}&\Vert \partial _{x_i} ({g(\rho _\gamma )\,h(\varPhi _\gamma )})\Vert _{L^2(Q_T)} \nonumber \\&\quad \le C \, ({\Vert g\Vert _{L^\infty (\mathbb {R})}\,\Vert Dh\Vert _{L^\infty (\mathbb {R}^{2 \times 2})} \Vert \nabla ^2(\phi _\gamma + \phi _K({\varvec{x}}_\gamma ;\cdot ))\Vert _{L^2(\varOmega _T)}} \nonumber \\&\qquad + {\Vert g'\Vert _{L^\infty (\mathbb {R})}\Vert \nabla \rho _\gamma \Vert _{L^2(Q_T)}\Vert h\Vert _{L^\infty (\mathbb {R}^2)}}) \le C \end{aligned}$$

(3.28)

for $i=1,2$, as well as

$$\begin{aligned}&\Vert \partial _t ({g(\rho _\gamma )\,h(\varPhi _\gamma )})\Vert _{L^2(Q_T)} \nonumber \\&\quad \le C \, ({\Vert g\Vert _{L^\infty (\mathbb {R})}\Vert Dh\Vert _{L^\infty (\mathbb {R}^{2 \times 2})} \Vert \partial _t \nabla (\phi _\gamma + \phi _K({\varvec{x}}_\gamma ;\cdot ))\Vert _{L^2(Q_T)}} \nonumber \\&\qquad + {\Vert g'\Vert _{L^\infty (\mathbb {R})}\Vert \partial _t \rho _\gamma \Vert _{L^2(Q_T)}\Vert h\Vert _{L^\infty (\mathbb {R}^2)}}) \le C. \end{aligned}$$

(3.29)

This follows from the regularity already shown for $\rho _\gamma$, $\phi _\gamma$ (see Lemma 3.1) and ${\varvec{x}}_\gamma$ (see Lemma 3.6). With these considerations we conclude

$$\begin{aligned} -g(\rho _\gamma )\,h(\varPhi _\gamma ) - \eta \,\rho _\gamma \,{\chi _{\partial \varOmega _O }}\in W^{1,1}_2(Q_T) \hookrightarrow W^{1,1/2}_2(Q_T). \end{aligned}$$

This allows us to apply the trace Lemma 2.3, which provides

$$\begin{aligned} \Vert r\Vert _{W^{1/2,1/4}_2(\varSigma _T)} \le C_r . \end{aligned}$$

(3.30)

Collecting the properties (3.26)–(3.30) an application of Thm. 2.1 in [22] and Lemma 3.1 implies

$$\begin{aligned} \rho _\gamma \in W^{2,1}_2(Q_T) \quad \text {and} \quad \phi _\gamma \in L^\infty (0,T;H^3(\varOmega ))\cap H^1(0,T;H^1(\varOmega )) . \end{aligned}$$

We can even further improve the regularity of $\rho _\gamma$. Analogous to (3.28) and (3.29) we show

$$\begin{aligned} \Vert \partial _{x_i} ({g(\rho _\gamma )\,h(\varPhi _\gamma )})\Vert _{L^p(Q_T)} + \Vert \partial _t ({\partial _{x_i}g(\rho _\gamma )\,h(\varPhi _\gamma )})\Vert _{L^2(0,T;L^p(\varOmega ))} \le C \end{aligned}$$

and together with the trace Lemma 2.3 we deduce

$$\begin{aligned} r&\in L^p(0,T;W^{1,p}(\varOmega ))\cap W^{1,2}(0,T;L^p(\varOmega ))\\&\hookrightarrow W^{1,1/2}_p(Q_T) \hookrightarrow W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T) . \end{aligned}$$

This, (3.26) and (3.27) allow a further application of Thm. 2.1 in [22] with $p \le 4$, from which we infer the desired regularity

$$\begin{aligned} \rho _\gamma \in W^{2,1}_p(Q_T). \end{aligned}$$

The functions $\rho _\gamma$, $\phi _\gamma$ and ${\varvec{x}}_\gamma$ are thus a strong solution of the system

$$\begin{aligned}&\partial _t\rho _\gamma - \nabla \cdot ({\rho _\gamma \,\beta (\rho _\gamma ,\phi _\gamma ,{\varvec{x}}_\gamma )}) = \varepsilon \,\varDelta \rho _\gamma \quad \text {in } Q_T , \end{aligned}$$

(3.31)

$$\begin{aligned}&- \delta _1\,\varDelta \phi _\gamma + |\nabla \phi _\gamma |^2 = \frac{1}{f( \rho _\gamma )^2+\delta _2} \quad \text {in } Q_T , \end{aligned}$$

(3.32)

$$\begin{aligned}&\dot{x}_{i,\gamma } = f({(\overline{{\varphi _\gamma } *\rho _\gamma })(\cdot ,x_{i,\gamma }(\cdot ))}) \, u_i \quad \text {in } [0,T], \quad i = 1, \ldots , M \end{aligned}$$

(3.33)

that satisfy the boundary and initial conditions pointwise almost everywhere.

Step 3: Limit $\underline{\gamma \rightarrow 0}$: In order to recover a solution to (1.9), it remains to pass to the limit $\gamma \rightarrow 0$ in (3.31)–(3.33). As a first step in this direction, note that b(t, x), c(t, x) and r(t, x) defined in (3.24) and (3.25) are bounded independently of $\gamma$. Thus, understanding (3.31) as the linear equation (3.23) we obtain the following estimate, uniformly in $\gamma$,

$$\begin{aligned} \Vert \rho _\gamma \Vert _{W^{2,1}_p(Q_T)} \le C . \end{aligned}$$

This implies the existence of a sequence $\gamma _k$ with $\gamma _k\rightarrow 0$ that

$$\begin{aligned} \rho _{\gamma _k} \rightharpoonup \rho \text { in } W^{2,1}_p(Q_T) \quad \text {and} \quad \rho _{\gamma _k} \rightarrow \rho \text { in } L^p(0,T;W^{1,p}(\varOmega )), \end{aligned}$$

(3.34)

where the second convergence is a consequence of the Aubin-Lions lemma, Thm. II.5.16 in [10]. As, moreover, $\phi _\gamma$ is also uniformly bounded in $L^p(0,T;W^{2,p}(\varOmega ))$ by $C_\phi$, we also have

$$\begin{aligned} \begin{aligned} D^2 \phi _{\gamma _k}&\rightharpoonup D^2 \phi{} & {} \text {in } L^p(0,T;L^p(\varOmega ;\mathbb {R}^{d \times d})) , \\ \phi _{\gamma _k}&\rightarrow \phi{} & {} \text {in } L^p(0,T;W^{1,p}(\varOmega )) , \end{aligned} \end{aligned}$$

(3.35)

where the second convergence follows from Lemma A.3. We omit the index k in the following to shorten the notation. Passing to the limit $\gamma \rightarrow 0$, in the sense of distributions, in equation (1.9b) yields the validity of this equation also for the limit values $\phi$ and $\rho$. As for (1.9c), denote by $x_i(t)$ the solution of $\dot{x}_i(t) = f({\overline{\rho }(x_i(t),t)})\,u_i(t)$, $t \in [0,T]$, and $x_i(0) = x_{i,0}$, with $\rho$ denoting the limit from (3.34). Arguing similarly as in Lemma 3.6 and using $|u_i(t)| \le 1$ for a.a. $t\in (0,T)$ we obtain

$$\begin{aligned}&|x_i(t) - x_{i,\gamma }(t)| \\&\quad \le \int _0^t |f({\overline{\rho }(s,x_i(s))}) - f({(\overline{{\varphi _\gamma } *\rho _\gamma })(s,x_{i,\gamma }(s))})| \mathop {}\!d s\\&\quad \le \int _0^t |f({\overline{\rho }(s,x_i(s))}) - f({\overline{\rho }(s,x_{i,\gamma }(s))})| \mathop {}\!d s\\&\qquad + \int _0^t |f({\overline{\rho }(s,x_{i,\gamma }(s))}) - f({(\overline{{\varphi _\gamma } *\rho })(s,x_{i,\gamma }(s))})| \mathop {}\!d s\\&\qquad + \int _0^t |f({(\overline{{\varphi _\gamma } *\rho })(s,x_{i,\gamma }(s))}) - f({(\overline{{\varphi _\gamma } *\rho _\gamma })(s,x_{i,\gamma }(s))})| \mathop {}\!d s\\&\quad \le \int _0^t L_{f,\rho ,{E }}(s)\,|x_i(s)-x_{i,\gamma }(s)|\mathop {}\!d s\\&\qquad + L_f({t\,\Vert \rho - {\varphi _\gamma } *\rho \Vert _{C([0,T] \times \overline{\varOmega })} + \sqrt{t}\,\Vert {\varphi _\gamma }\Vert _{L^1(\varOmega )}\,\Vert \rho - \rho _\gamma \Vert _{L^2(0,T;L^\infty (\varOmega ))}}) . \end{aligned}$$

As $\Vert {\varphi _\gamma }\Vert _{L^1(\varOmega )}=1$, another application of Gronwall’s inequality yields

$$\begin{aligned}&|x_i(t) - x_{i,\gamma }(t)| \\&\quad \le C \, L_f ({t\,\Vert \rho - {\varphi _\gamma } *\rho \Vert _{C([0,T] \times \overline{\varOmega })} + \sqrt{t}\,\Vert \rho - \rho _\gamma \Vert _{L^2(0,T;L^\infty (\varOmega ))}})\\&\quad \cdot \exp \left( {\int _0^t L_{f,\rho ,{E }}(s) \mathop {}\!d s}\right) , \end{aligned}$$

and due to (2.3), (3.34) we see that as $\gamma \rightarrow 0$, the right-hand side converges to 0 and thus

$$\begin{aligned} x_{i,\gamma }(t) \rightarrow x_i(t)\;\text { for every } t \in (0,T), \quad i = 1, \ldots , M. \end{aligned}$$

(3.36)

Next, we pass to the limit in (3.31). The convergence of the linear terms is a direct consequence of (3.34). For the convection term we note that it can be written as

$$\begin{aligned} \nabla \cdot ({\rho _\gamma \,\beta (\rho _\gamma , \phi _\gamma , {\varvec{x}}_\gamma )})&= \nabla \cdot ({g(\rho _\gamma )\, h(\varPhi _\gamma )}) \\&= g'(\rho _\gamma )\,\nabla \rho _\gamma \cdot h(\varPhi _\gamma ) + g(\rho _\gamma )\, Dh(\varPhi _\gamma ) :\nabla ^2(\phi _\gamma + \phi _K({\varvec{x}}_\gamma ;\cdot )) . \end{aligned}$$

As both $\rho _\gamma$ and h are uniformly bounded, the convergences (3.34), (3.35) and (3.36) imply

$$\begin{aligned} \nabla \cdot (\rho _\gamma \,\beta (\rho _\gamma , \phi _\gamma , {\varvec{x}}_\gamma )) \rightharpoonup \nabla \cdot (\rho \,\beta (\rho , \phi , {\varvec{x}}))\quad \text { in } L^p(Q_T) . \end{aligned}$$

It remains to pass to the limit in the boundary conditions. The trace Lemma 2.3 implies that $\varepsilon \,\nabla \rho _\gamma \cdot n$ is bounded in $L^p(\varSigma _T)$ and thus converges weakly. The compact embedding $W^{2-1/p,1-1/(2p)}_p(\varSigma _T) \hookrightarrow L^p(\varSigma _T)$ implies the convergences

$$\begin{aligned} {\rho }{|}_{\varSigma _T} \rightarrow {\rho }{|}_{\varSigma _T} \quad \text {and} \quad \nabla {\phi _\gamma }{|}_{\varSigma _T} \rightarrow {\nabla \phi }{|}_{\varSigma _T}\;\text { in } L^p(\varSigma _T) , \end{aligned}$$

so that using the uniform boundedness of both $\rho _\gamma$ and h as well as (3.36) we conclude

$$\begin{aligned} \varepsilon \,\nabla \rho _\gamma \cdot n + \rho _\gamma \,\beta (\rho _\gamma ,\phi _\gamma ,{\varvec{x}}_\gamma )\cdot n \rightharpoonup \varepsilon \,\nabla \rho \cdot n + \rho \,\beta (\rho ,\phi ,{\varvec{x}})\cdot n\;\text { in } L^p(\varSigma _T) . \end{aligned}$$

Thus, the weak limit of $\rho _\gamma$, namely $\rho$, is the strong solution of (1.9) with boundary conditions (1.5). This completes the proof.

The previous theorem allows us to introduce the solution operator

$$\begin{aligned} S :\mathcal {U}\rightarrow \mathcal {Y}, \quad {\varvec{u}}\mapsto S({\varvec{u}}) :=(\rho ,\phi ,{\varvec{x}}) \end{aligned}$$

of (1.9) with boundary conditions (1.5) and initial conditions ${\varvec{x}}(0) = {\varvec{x}}^0$ and $\rho (\cdot ,0) = \rho _0$. The control and state spaces are

$$\begin{aligned} \mathcal {U}&:=L^\infty (0,T;\mathbb {R}^2)^M , \\ \mathcal {Y}&:=W^{2,1}_p(Q) \times ({L^\infty (0,T;W_\text {ND}^{2,p}(\varOmega ))\cap W^{1,p}(0,T;W^{1,p}(\varOmega ))}) \times W^{1,s}(0,T;\mathbb {R}^2)^M \end{aligned}$$

with $1/s=1/2+1/p$. Later, the operator S is referred to as control-to-state operator.

We conclude this section with an auxiliary result required later to show the existence of global solutions to an optimal control problem.

Lemma 3.9

The operator $S :\mathcal {U}\rightarrow \mathcal {Y}$ is weakly sequentially continuous.

Proof

Given a weakly convergent sequence $(u_n)_{n\in \mathbb {N}} \subset \mathcal {U}$ with $u_n\rightharpoonup u$ one has to show that the corresponding states $(\rho _n,\phi _n,{\varvec{x}}_n) = S(u_n)$ converge weakly in $\mathcal {Y}$ to $(\rho ,\phi ,{\varvec{x}}) = S(u)$. This follows from the same arguments as those in step 3 of the proof of Theorem 3.8, together with the uniqueness of solutions.

4 The Linearized System

In order to prove necessary optimality conditions for the optimal control problems introduced later, we investigate the differentiability of the control-to-state operator S. The desired results follow from the implicit function theorem applied to the equation $e({\varvec{y}},{\varvec{u}}) = 0$ with ${\varvec{y}}= S({\varvec{u}}) = (\rho ,\phi ,{\varvec{x}})$. Here, e corresponds to the strong formulation of our forward system (1.9a)–(1.9c), more precisely, there holds

$$\begin{aligned} e :\mathcal {Y}\times \mathcal {U}\rightarrow \mathcal {Z}\end{aligned}$$

with

$$\begin{aligned} \mathcal {Z}= & {} L^p(Q_T) \times W^{1-1/p,1/2-1/(2p)}_p(\varSigma _W ) \times W^{2(1-1/p),p}(\varOmega ) \\{} & {} \times ({L^\infty (0,T;L^p(\varOmega ))\cap W^{1,p}(0,T;W^{1,p}(\varOmega )^*)}) \times W^{1,s}(0,T;\mathbb {R}^2)^M \end{aligned}$$

defined by

$$\begin{aligned} e_1({\varvec{y}},{\varvec{u}})&:=\partial _t\rho - \varepsilon \,\varDelta \rho - \nabla \cdot (\rho \,\beta ({\varvec{y}})) \nonumber \\ e_2({\varvec{y}},{\varvec{u}})&:=(\varepsilon \nabla \rho + \rho \,\beta ({\varvec{y}}))\cdot n + {\chi _{\partial \varOmega _O }}\,\eta \,\rho \nonumber \\ e_3({\varvec{y}},{\varvec{u}})&:=\rho (0)-\rho _0\\ e_4({\varvec{y}},{\varvec{u}})&:=-\delta _1\,\varDelta \phi + |\nabla \phi |^2 - \frac{1}{f(\rho )^2+\delta _2} \nonumber \\ e_{5,i}({\varvec{y}},{\varvec{u}})(t)&:=x_i(t) - x_{i,0} - \int _0^t f(\overline{\rho }(s,x_i(s)))\,u_i(s) \mathop {}\!d s, \quad i=1,\ldots ,M . \nonumber \end{aligned}$$

(4.1)

Notice that from here on we write $\beta ({\varvec{y}})$ in place of $\beta (\rho ,\phi ,{\varvec{x}})$. Recall that $p\in (2,4)$ is a fixed number and the integrability index s of the space of agent trajectories is chosen such that $\frac{1}{s} = \frac{1}{2} + \frac{1}{p}$.

Lemma 4.1

Let ${\varvec{y}}= (\rho ,\phi ,{\varvec{x}}) \in \mathcal {Y}$ and define $\varPhi := \nabla (\phi +\phi _K({\varvec{x}};\cdot ))$. Then a constant $C_h>0$ exists such that

$$\begin{aligned} \Vert h(\varPhi )\Vert _{L^\infty (Q_T)} + \Vert h(\varPhi )\Vert _{L^\infty (0,T;W^{1,p}(\varOmega ))} +\Vert \partial _t h(\varPhi )\Vert _{L^p(Q_T)}&\le C_h \\ \Vert D^kh(\varPhi )\Vert _{L^\infty (Q_T)} + \Vert D^kh(\varPhi )\Vert _{L^\infty (0,T;W^{1,p}(\varOmega ))}&\le C_h, \quad k=1,2 . \end{aligned}$$

Note that

$$\begin{aligned} Dh({\varvec{y}}) = \chi _{\{|{\varvec{y}}| \le 1\},\varepsilon } \, id \end{aligned}$$

is a smoothed characteristic function. Moreover, the functions h and $D^k h$, $k=1,2$, are Lipschitz continuous. In particular, for any two vector fields $\varPhi _1, \varPhi _2 \in L^\infty (0,T;W^{1,p}(\varOmega ))$ there holds

$$\begin{aligned}&\Vert \nabla \cdot (Dh(\varPhi _1)-Dh(\varPhi _2))\Vert _{L^\infty (0,T;L^p(Q_T))} \nonumber \\&\quad \le C_h ({1+\Vert D\varPhi _1\Vert _{L^\infty (Q_T)}}) \Vert D\varPhi _1-D\varPhi _2\Vert _{L^\infty (0,T;L^p(\varOmega ))} . \end{aligned}$$

(4.2)

Proof

The result follows directly from the regularity properties of h.

To shorten the notation we also introduce the following constants,

$$\begin{aligned} C_\rho&:=\Vert \rho \Vert _{W^{2,1}_p(Q_T)} ,&C_\phi&:=\Vert \phi \Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} + \Vert \phi \Vert _{W^{1,p}(0,T;W^{1,p}(\varOmega ))} ,&\nonumber \\ C_{\varvec{x}}&:=\Vert {\varvec{x}}\Vert _{W^{1,s}(0,T)} , \end{aligned}$$

(4.3)

whose boundedness is guaranteed due to assumptions (A4), (K1), (K2) and (C1).

We confirm the assumptions of the implicit function theorem in the following lemmas.

Lemma 4.2

The operator $e :\mathcal {Y}\times \mathcal {U}\rightarrow \mathcal {Z}$ is continuously Fréchet differentiable.

Proof

We start with equation $e_1$. Again, we write the nonlinear term in the form

$$\begin{aligned} \rho \,\beta ({\varvec{y}}) = g(\rho )\,h(\varPhi ) \text { with } \varPhi = \nabla (\phi +\phi _K({\varvec{x}};\cdot )) . \end{aligned}$$

For given $(\rho ,\phi , {\varvec{x}}) \in \mathcal {Y}$ and $(\widetilde{\rho },\widetilde{\phi }, \widetilde{{\varvec{x}}}) \in \mathcal {Y}$, Taylor’s formula with integral remainder yields

$$\begin{aligned}&e_1(\rho +\widetilde{\rho },\phi ,{\varvec{x}};{\varvec{u}}) - e_1(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) \\&\quad = \partial _t\widetilde{\rho } - \varepsilon \varDelta \widetilde{\rho } - \nabla \cdot (g'(\rho )\,\widetilde{\rho }\,h(\varPhi )) - \,\int _0^1 \nabla \cdot [{(g'(\rho +s\,\widetilde{\rho }) - g'(\rho ))\,\widetilde{\rho }\,h(\varPhi )}]\mathop {}\!d s, \\&e_1(\rho ,\phi +\widetilde{\phi },{\varvec{x}};{\varvec{u}}) - e_1(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) \\&\quad = -\nabla \cdot ({g(\rho )\,Dh(\varPhi )\,\nabla \widetilde{\phi }}) - \int _0^1 \nabla \cdot [{g(\rho ) ({Dh(\widetilde{\varPhi }_\phi ) - Dh(\varPhi )}) \nabla \widetilde{\phi }}] \mathop {}\!d s, \\&e_1(\rho ,\phi ,{\varvec{x}}+\widetilde{{\varvec{x}}};{\varvec{u}}) - e_1(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) \\&\quad = \nabla \cdot ({g(\rho )\,Dh(\varPhi )\nabla ^2K(\cdot -x_i)\,\widetilde{x}_i}) \\&\qquad + \int _0^1 \nabla \cdot ({g(\rho ) [{Dh(\widetilde{\varPhi }_{\varvec{x}})\,\nabla ^2K(\cdot -x_i-s\,\widetilde{x}_i) - Dh(\varPhi )\,\nabla ^2K(\cdot -x_i)}] \widetilde{x}_i}) \mathop {}\!d s, \end{aligned}$$

for $\widetilde{{\varvec{x}}} = (0,\ldots 0, \widetilde{x}_i, 0, \ldots 0)^T$, $i=1,\ldots ,M$.

In the above equations we used the notation $\widetilde{\varPhi }_\phi = \nabla (\phi + s\,\widetilde{\phi } + \phi _K({\varvec{x}};\cdot ))$ and $\widetilde{\varPhi }_{\varvec{x}}= \nabla (\phi +\phi _K({\varvec{x}}+s\,\widetilde{{\varvec{x}}};\cdot ))$. In the following we derive bounds for the remainder terms (the terms depending nonlinearly on $\widetilde{\rho }$, $\widetilde{\phi }$ and $\widetilde{{\varvec{x}}}$ in the above equations), which we denote by $r_{1,\rho }(\widetilde{\rho })$, $r_{1,\phi }(\widetilde{\phi })$ and $r_{1,x_i}(\widetilde{{\varvec{x}}})$, $i=1,\ldots ,M$.

First, we apply the product rule, the Hölder inequality and the embedding $W^{2,1}_p(Q_T) \hookrightarrow L^\infty (Q_T)$, where we denote by $C_\infty$ the maximum of the embedding constant and 1. These arguments yield

$$\begin{aligned}&\Vert r_{1,\rho }(\widetilde{\rho })\Vert _{L^p(Q_T)} \nonumber \\&\quad \le \int _0^1 \Big ([{\Vert (g''(\rho +s\,\widetilde{\rho }) - g''(\rho ))\Vert _{L^\infty (Q_T)}\Vert \nabla \rho \Vert _{L^p(Q_T)}} \nonumber \\&\qquad \qquad + {s\,\Vert g''(\rho +s\,\widetilde{\rho })\Vert _{L^\infty (Q_T)}\Vert \nabla \widetilde{\rho }\Vert _{L^p(Q_T)}}] \cdot \Vert \widetilde{\rho }\Vert _{L^\infty (Q_T)}\Vert h(\varPhi )\Vert _{L^\infty (Q_T)} \nonumber \\&\quad + \Vert g'(\rho +s\,\widetilde{\rho })-g'(\rho )\Vert _{L^\infty (Q_T)} \nonumber \\&\qquad \cdot ({\Vert \nabla \widetilde{\rho }\Vert _{L^p(Q_T)}\,\Vert h(\varPhi )\Vert _{L^\infty (Q_T)} + \Vert \widetilde{\rho }\Vert _{L^\infty (Q_T)} \Vert \nabla \cdot h(\varPhi )\Vert _{L^p(Q_T)}})\Big )\mathop {}\!d s\nonumber \\&\quad \le C_g\,C_h({C_\rho \,+3\,C_\infty }) \,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}^2 = \textrm{o} ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}} ). \end{aligned}$$

(4.4)

Second, to show differentiability with respect to $\phi$ we confirm

$$\begin{aligned}&r_{1,\phi }(\widetilde{\phi }) = \int _0^1 \Big [g'(\rho )\nabla \rho ({Dh(\widetilde{\varPhi }_\phi ) - Dh(\varPhi )})\nabla \widetilde{\phi } \\&\quad + g(\rho ) ({s\,\nabla \cdot ({Dh(\widetilde{\varPhi }_\phi ) - Dh(\varPhi )}) \cdot \nabla \widetilde{\phi } + ({Dh(\widetilde{\varPhi }_\phi )-Dh(\varPhi )}) :\nabla ^2\widetilde{\phi }}) \Big ] \mathop {}\!d s. \end{aligned}$$

With the Hölder inequality, the Lipschitz properties of Dh and $D^2 h$, in particular (4.2), and the usual embeddings we can show

$$\begin{aligned} \begin{aligned} \Vert r_{1,\phi }(\widetilde{\phi })\Vert _{L^p(Q_T)}&\le C_g \, C_h \, C_\infty \, (C_\rho + 2) \, \Vert \widetilde{\phi }\Vert ^2_{L^\infty (0,T;W^{2,p}(\varOmega ))} \\&= \mathcal {O}({\Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))}}) . \end{aligned} \end{aligned}$$

Third, we derive an estimate for $r_{1,x_i}(\widetilde{{\varvec{x}}})$ with $\widetilde{{\varvec{x}}} = (0,\ldots ,0,\widetilde{x}_i,0,\ldots ,0)$. We use the notation $\phi _K = \phi _K({\varvec{x}};\cdot )$ and $\widetilde{\phi }_K = \phi _K({\varvec{x}}+s\widetilde{{\varvec{x}}};\cdot )$ as well as $K_i = K(\cdot - x_i)$ and $\widetilde{K}_i = K(\cdot -x_i-s\,\widetilde{x}_i)$, reformulate the remainder term by applying the product and chain rule and obtain

$$\begin{aligned}&r_{1,x_i}(\widetilde{x}_i) \\&\quad \le \int _0^1 \Big (|g'(\rho )\,\nabla \rho | \, ({|Dh(\widetilde{\varPhi }_{\varvec{x}})|\, |\nabla ^2 K_i - \nabla ^2 \widetilde{K}_i| + |Dh(\widetilde{\varPhi }_{\varvec{x}})- Dh(\varPhi )|\,|\nabla ^2 K_i|}) \\&\qquad + g(\rho ) \Big [|\nabla \cdot Dh(\widetilde{\varPhi }_{\varvec{x}})|\, |\nabla ^2\widetilde{K}_i - \nabla ^2 K_i| + |Dh(\widetilde{\varPhi })||\nabla ^3\widetilde{K}_i - \nabla ^3K_i| \\&\qquad + |\nabla ^3 K_i|\,|Dh(\widetilde{\varPhi }_{\varvec{x}})-Dh(\varPhi )| + |\nabla ^2 K_i|\,|\nabla \cdot (Dh(\widetilde{\varPhi }_{\varvec{x}})-Dh(\varPhi ))|\Big ] \Big ) \, |\widetilde{x}_i| \mathop {}\!d s. \end{aligned}$$

Here, $|\cdot |$ is an arbitrary vector, matrix or tensor norm, depending on the argument. With the Hölder inequality, the regularity of $(\rho ,\phi ,{\varvec{x}})$, in particular $\rho \in W^{2,1}_p(Q_T) \hookrightarrow L^\infty (0,T;W^{1,p}(\varOmega ))\cap L^p(0,T;W^{1,\infty }(\varOmega ))$, and the Lipschitz properties of h and K we deduce

$$\begin{aligned} \begin{aligned} \Vert r_{1,x_i}(\widetilde{{\varvec{x}}})\Vert _{L^p(Q_T)}&\le C_g\,(C_\rho +4)\,C_h\,C_{\phi _K}\,2\,(1+C_{\phi _K})\Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)}^2 \\&= \textrm{o} ({\Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)}}) . \end{aligned} \end{aligned}$$

The differentiability of $e_2$ can be shown with similar arguments. The Taylor formula yields

$$\begin{aligned}&e_2(\rho +\widetilde{\rho },\phi ,{\varvec{x}};{\varvec{u}}) - e_2(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) \\&\quad = (\varepsilon \nabla \widetilde{\rho } + \widetilde{\rho }\,\beta ({\varvec{y}})) \cdot n + {\chi _{\partial \varOmega _O }}\,\eta \,\widetilde{\rho } \\&\quad + g'(\rho )\,\widetilde{\rho }\,h(\varPhi )\cdot n + \int _0^1 (g'(\rho +s\,\widetilde{\rho })-g'(\rho ))\,\widetilde{\rho }\,h(\varPhi )\cdot n\mathop {}\!d s, \\&e_2(\rho ,\phi +\widetilde{\phi },{\varvec{x}};{\varvec{u}}) - e_2(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) \\&\quad = g(\rho )\,Dh(\varPhi )\,\nabla \widetilde{\phi }\,\cdot n + \int _0^1 g(\rho )({Dh(\widetilde{\varPhi }_\phi ) - Dh(\varPhi )})\nabla \widetilde{\phi }\cdot n \mathop {}\!d s, \\&e_2(\rho ,\phi ,{\varvec{x}}+\widetilde{{\varvec{x}}};{\varvec{u}}) - e_2(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) \\&\quad = -g(\rho )\,Dh(\varPhi )\,\nabla ^2K(\cdot -x_i)\,\widetilde{x}_i\cdot n \\&\qquad -\int _0^1 g(\rho )({Dh(\widetilde{\varPhi }_{{\varvec{x}}})\,\nabla ^2K(\cdot - x_i-s\,\widetilde{x}_i) - Dh(\varPhi )\,\nabla ^2 K(\cdot -x_i)})\,\widetilde{x}_i\cdot n \mathop {}\!d s, \end{aligned}$$

and it remains to estimate the remainder terms, i. e., the last terms on the right-hand sides of the previous equations. We abbreviate these terms by $r_{2,\rho }(\widetilde{\rho })\cdot n$, $r_{2,\phi }(\widetilde{\phi })\cdot n$ and $r_{2,{\varvec{x}}}(\widetilde{{\varvec{x}}})\cdot n$, respectively.

First, we show $\Vert r_{2,\rho }(\widetilde{\rho })\cdot n\Vert _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T)} = \textrm{o} ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}})$. We wish to apply the trace Lemma 2.3, for which we have to show the required time and space regularity of the extension onto $Q_T$. First, note that there holds

$$\begin{aligned} \Vert r_{2,\rho }(\widetilde{\rho })\Vert _{L^p(0,T;W^{1,p}(\varOmega ))} = \textrm{o} ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}}) , \end{aligned}$$

(4.5)

which can be concluded from the same arguments as in (4.4). Moreover, for the time derivative we show

$$\begin{aligned} \partial _t r_{2,\rho }(\widetilde{\rho })&= \int _0^1 \Big ( [{(g''(\rho +s\widetilde{\rho }) - g''(\rho ))\,\partial _t\rho + s\,g''(\rho +s\widetilde{\rho })\,\partial _t \widetilde{\rho }}] \widetilde{\rho }\,h(\varPhi ) \\&\qquad + ({g'(\rho +s\,\widetilde{\rho })- g'(\rho )}) \, ({\partial _t\widetilde{\rho }\,h(\varPhi ) + \widetilde{\rho }\,\partial _t h(\varPhi )}) \Big ) \mathop {}\!d s\end{aligned}$$

and with the usual arguments we obtain for $s^{-1} = 2^{-1}+p^{-1}$

$$\begin{aligned} \Vert r_{2,\rho }(\widetilde{\rho })\Vert _{W^{1,s}(0,T;L^p(\varOmega ))} \le C_\infty \,C_g\,C_h\,(3+C_\rho )\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}^2 = \textrm{o} ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}}) . \end{aligned}$$

(4.6)

With the estimates (4.5) and (4.6), the embedding

$$\begin{aligned} L^p(0,T;W^{1,p}(\varOmega ))\cap W^{1,s}(0,T;L^p(\varOmega )) \hookrightarrow W^{1,1/2}_p(Q_T) \end{aligned}$$

(4.7)

and the trace Lemma 2.3 we deduce

$$\begin{aligned} \Vert r_{2,\rho }(\widetilde{\rho })\cdot n\Vert _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T)} = \textrm{o} ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}}) , \end{aligned}$$

which implies the differentiability of $e_2$ with respect to $\rho$.

To show an estimate for $r_{2,\phi }(\widetilde{\phi })\cdot n$ we proceed in a similar fashion. With analogous arguments we deduce the estimates

$$\begin{aligned} \Vert \partial _{x_i} r_{2,\phi }(\widetilde{\phi })\Vert _{L^p(Q_T)}&\le C\,\Vert \nabla \widetilde{\phi }\Vert _{L^\infty (0,T;W^{1,p}(\varOmega ))}^2 , \\ \Vert \partial _t r_{2,\phi }(\widetilde{\phi })\Vert _{L^s(0,T;L^p(Q_T))}&\le C\left( C_\rho \,\Vert \nabla \widetilde{\phi }\Vert _{L^\infty (Q_T)}^2 + \Vert \partial _t \nabla \widetilde{\phi }\Vert _{L^p(Q_T)}\,\Vert \nabla \widetilde{\phi }\Vert _{L^\infty (Q_T)}\right) \end{aligned}$$

and exploit again (4.7) to arrive at

$$\begin{aligned} \Vert r_{2,\phi }(\widetilde{\phi })\Vert _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T)} = \mathcal {O}({\Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} + \Vert \widetilde{\phi }\Vert _{W^{1,p}(0,T;W^{1,p}(\varOmega ))}}) , \end{aligned}$$

which confirms the differentiability of $e_2$ w.r.t. $\phi$.

An analogous procedure is used to deduce an estimate for the remainder term $r_{2,{\varvec{x}}}(\widetilde{{\varvec{x}}})\cdot n$. For each direction $\widetilde{{\varvec{x}}} = (0,\ldots 0, \widetilde{x}_i, 0, \ldots 0)^T$, $i=1,\ldots ,M$, a direct calculation taking into account Lemma 4.1 and the Lipschitz continuity of the derivatives of K and h yields

$$\begin{aligned} \Vert \partial _{x_i} r_{2,{\varvec{x}}}(\widetilde{{\varvec{x}}})\Vert _{L^p(Q_T)}&\le C\,\Vert \widetilde{x}_i\Vert ^2_{L^\infty (0,T;\mathbb {R}^2)} , \\ \Vert \partial _t r_{2,{\varvec{x}}}(\widetilde{{\varvec{x}}})\Vert _{L^s(0,T;L^p(\varOmega ))}&\le C\,\Vert \widetilde{x}_i\Vert _{W^{1,s}(0,T;\mathbb {R}^2)}^2 . \end{aligned}$$

Using again (4.7) and the trace Lemma 2.3 leads to

$$\begin{aligned} \Vert r_{2,x_i}(\widetilde{{\varvec{x}}})\cdot n\Vert _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T)} = \mathcal {O}({\Vert \widetilde{x}_i\Vert _{W^{1,s}(0,T;\mathbb {R}^2)}}) , \end{aligned}$$

which proves the differentiability of $e_2$ w.r.t. ${\varvec{x}}$.

The component $e_3$ is trivially differentiable. For $e_4$ we show

$$\begin{aligned}&e_4(\rho +\widetilde{\rho },\phi ,{\varvec{x}}) - e_4(\rho ,\phi ,{\varvec{x}}) \\&\quad = \frac{2 f(\rho )\,f'(\rho )}{(f^2(\rho )+\delta _2)^2} \widetilde{\rho } + \int _0^1\left( {\frac{2\,f(\rho +s\widetilde{\rho })\,f'(\rho +s\widetilde{\rho })}{(f(\rho +s\widetilde{\rho })^2+\delta _2)^2}-\frac{2\,f(\rho )\,f'(\rho )}{(f(\rho )^2+\delta _2)^2}}\right) \widetilde{\rho } \mathop {}\!d s, \\&e_4(\rho ,\phi +\widetilde{\phi },{\varvec{x}}) - e_4(\rho ,\phi ,{\varvec{x}}) \\&\quad = \delta _1\varDelta \widetilde{\phi } + 2\nabla \phi ^T \nabla \widetilde{\phi } + |\nabla \widetilde{\phi }|^2 . \end{aligned}$$

Again, we denote the remainder terms (the terms which are nonlinear in $\widetilde{\rho }$ and $\widetilde{\phi }$) by $r_{4,\rho }(\widetilde{\rho })$ and $r_{4,\phi }(\widetilde{\phi })$. We show

$$\begin{aligned} \Vert r_{4,\phi }(\widetilde{\phi })\Vert _{L^\infty (0,T;L^p(\varOmega ))}&\le C\Vert \nabla \widetilde{\phi }\Vert _{L^\infty (Q_T)}^2 \le C\Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))}^2 , \\ \Vert \partial _t r_{4,\phi }(\widetilde{\phi })\Vert _{L^p(0,T;W^{1,p}(\varOmega )^*)}&\le C\Vert \nabla \widetilde{\phi }\Vert _{L^\infty (0,T;W^{1,p}(\varOmega ))}\,\Vert \partial _t\nabla \widetilde{\phi }\Vert _{L^p(Q_T)} . \end{aligned}$$

and both estimates together imply

$$\begin{aligned} \Vert r_{4,\phi }(\widetilde{\phi })\Vert _{\mathcal {Z}_4} = \mathcal {O}({\Vert \widetilde{\phi }\Vert _{\mathcal {Y}_2}}) . \end{aligned}$$

With similar arguments we can show the differentiability of $e_4$ with respect to $\rho$.

Finally we consider $e_5$, which is nonlinear in $\rho$ and $x_i$. We confirm by a simple computation that

$$\begin{aligned}&e_5(\rho +\widetilde{\rho },\phi ,{\varvec{x}};{\varvec{u}};t) - e_5(\rho ,\phi ,{\varvec{x}};{\varvec{u}};t) \\&\quad = \int _0^t f'(\overline{\rho }(\tau ,x_i(\tau )))\,\overline{\widetilde{\rho }}(\tau ,x_i(\tau ))\,u_i(\tau ) \mathop {}\!d \tau \\&\qquad +\int _0^1 \int _0^t ({f'({(\overline{\rho +s\,\widetilde{\rho }})(\tau ,x_i(\tau ))} )- f'({\overline{\rho }(\tau ,x_i(\tau ))})}) \overline{\widetilde{\rho }}(\tau ,x_i(\tau ))\,u_i(\tau ) \mathop {}\!d \tau \mathop {}\!d s, \\&e_5(\rho ,\phi ,{\varvec{x}}+\widetilde{{\varvec{x}}};{\varvec{u}};t) - e_5(\rho ,\phi ,{\varvec{x}};{\varvec{u}};t) \\&\quad = \widetilde{x}_i(t) - \int _0^t f'(\rho (\tau ,x_i(\tau )))\,\nabla \rho (\tau ,x_i(\tau )) \cdot \widetilde{x}_i(\tau )\,u_i(\tau ) \mathop {}\!d \tau \\&\qquad - \int _0^1 \int _0^t \Big (f'(\overline{\rho }(\tau ,(x_i+s\,\widetilde{x}_i)(\tau )))\,\nabla \overline{\rho }(\tau ,(x_i+s\,\widetilde{x}_i)(\tau )) \\&\qquad - f'(\overline{\rho }(\tau ,x_i(\tau )))\,\nabla \overline{\rho }(\tau ,x_i(\tau ))\Big ) \cdot \widetilde{x}_i(\tau ) \, u_i(\tau ) \mathop {}\!d \tau \mathop {}\!d s\end{aligned}$$

for $\widetilde{{\varvec{x}}}=(0 \ldots \widetilde{x}_i\ldots 0)^T$, $i=1,\ldots ,M$. Again, the remainder terms (the terms depending nonlinearly on $\widetilde{\rho }$ and $\widetilde{{\varvec{x}}}$ in the above equation) are denoted by $r_{5,\rho }(\widetilde{\rho })$ and $r_{5,{\varvec{x}}}(\widetilde{{\varvec{x}}})$.

Using the Lipschitz continuity of $f'$, assumption (C1) and $\Vert \overline{\rho }\Vert _{L^\infty (\mathbb {R}^2)} \le \Vert \rho \Vert _{L^\infty (\varOmega )}$ (see Lemma 2.4), we obtain

$$\begin{aligned} |r_{5,\rho }(\widetilde{\rho };t)| \le C_f\,C_{{E },\infty }^2\,\Vert \widetilde{\rho }\Vert _{L^2(0,T;L^\infty (\varOmega ))}^2 = \textrm{o} ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}}) . \end{aligned}$$

(4.8)

With similar arguments and the Hölder inequality with $1/s>2/p$ we deduce for the temporal derivative of $r_{5,\rho }$

$$\begin{aligned}&\Vert \partial _t r_{5,\rho }(\widetilde{\rho })\Vert _{L^s(0,T;\mathbb {R}^2)} \\&\quad \le \int _0^1 \Vert {({f'({(\overline{\rho +s\widetilde{\rho }})(\cdot ,x_i(\cdot ))}) - f'(\overline{\rho }(\cdot ,x_i(\cdot )))}) \, \overline{\widetilde{\rho }}({\cdot ,x_i(\cdot )}) \, u_i(\cdot )\Vert }_{L^s(0,T;\mathbb {R}^2)} \mathop {}\!d s\\&\quad \le C_f\,\Vert \overline{\widetilde{\rho }}\Vert _{L^p(0,T;L^\infty (\mathbb {R}^2))}^2 \le C_f\,C_{{E },\infty }^2\,\Vert \widetilde{\rho }\Vert _{L^p(0,T;L^\infty (\varOmega ))}^2 = \textrm{o} ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}}) . \end{aligned}$$

For the remainder term $r_{5,{\varvec{x}}}(\widetilde{{\varvec{x}}})$ with $\widetilde{{\varvec{x}}} = (0,\ldots 0, \widetilde{x}_i, 0, \ldots 0)^T$ we show

$$\begin{aligned}&|r_{5,{\varvec{x}}}(\widetilde{{\varvec{x}}};t)| \nonumber \\&\quad \le \Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)} \cdot \nonumber \\&\qquad \cdot \int _0^1\int _0^t \Big (|f'({\overline{\rho }(\tau ,(x_i+s\,\widetilde{x}_i)(t))}) - f'(\overline{\rho }(\tau ,x_i(t)))| |\nabla \overline{\rho }(\tau , (x_i+s\,\widetilde{x}_i)(\tau ))| \nonumber \\&\quad \qquad + |f'(\overline{\rho }(\tau ,x_i(\tau )))||\nabla \overline{\rho }(\tau ,(x_i+s\,\widetilde{x}_i)(\tau )) - \nabla \overline{\rho }(\tau ,x_i(\tau ))| \Big ) \mathop {}\!d \tau \mathop {}\!d s. \end{aligned}$$

(4.9)

Using Lipschitz estimates for $f'$ and the mean value theorem we conclude

$$\begin{aligned} |f'(\overline{\rho }(\tau ,(x_i+s\,\widetilde{x}_i)(\tau ))) - f'(\overline{\rho }(\tau ,x_i(\tau )))|&\le C_f\,\text {Lip}(\overline{\rho }(\tau ,\cdot ))\,|\widetilde{x}_i(\tau )| \nonumber \\&\le C_f \, C_{E }\, \Vert \rho (\tau ,\cdot )\Vert _{W^{2,p}(\varOmega )}\,|\widetilde{x}_i(\tau )|,\end{aligned}$$

(4.10)

$$\begin{aligned} |\nabla \overline{\rho }(\tau ,(x_i+s\,\widetilde{x}_i)(\tau )) - \nabla \overline{\rho }(\tau ,x_i(\tau ))|&\le \Vert \nabla \overline{\rho }(\tau ,\cdot )\Vert _{C^{0,\alpha }(\mathbb {R}^2)}\,|\widetilde{x}_i(\tau )|^\alpha \nonumber \\&\le C_{E }\, \Vert \rho (\tau ,\cdot )\Vert _{W^{2,p}(\varOmega )}\,|\widetilde{x}_i(\tau )|^\alpha . \end{aligned}$$

(4.11)

In the above estimates we used $\Vert \overline{\rho }\Vert _{C^{0,1}(\mathbb {R}^2)} \le \Vert \overline{\rho }\Vert _{C^{1,\alpha }(\varOmega )}$, the continuity of ${E }:C^{1,\alpha }(\overline{\varOmega })\rightarrow C^{1,\alpha }(\mathbb {R}^2)$, see Lemma 2.4, and the embedding $W^{2,p}(\varOmega ) \hookrightarrow C^{1,\alpha }(\overline{\varOmega })$, which is valid for $\alpha \in (0,1/2]$ due to $p>2$.

The insertion of (4.10) and (4.11) into (4.9) yields

$$\begin{aligned} |r_{5,{\varvec{x}}}(\widetilde{{\varvec{x}}};t)|&\le C_{E }\, C_f \, C_\rho \, \max \{{t^{1-2/p}, t^{1-1/p}}\} \left( {C_\rho \, \Vert \widetilde{x}_i\Vert ^2_{L^\infty (0,T;\mathbb {R}^2)} + \Vert \widetilde{x}_i\Vert ^{1+\alpha }_{L^\infty (0,T;\mathbb {R}^2)}}\right) \\&= \textrm{o} ({\Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)}}) . \end{aligned}$$

The above also uses an application of the Hölder inequality in time and $2/p < 1$.

Finally, the Hölder inequality in time with $1/s = 1/2+1/p$ yields the estimate

$$\begin{aligned}&\Vert \partial _t r_{5,{\varvec{x}}}(\widetilde{{\varvec{x}}})\Vert _{L^s(0,T;\mathbb {R}^2)} \\&\quad \le C_f\Big (\Vert \overline{\rho }(\cdot ,(x_i+s\,\widetilde{x}_i)(\cdot )) - \overline{\rho }(\cdot ,x_i(\cdot ))\Vert _{L^2(0,T)}\, \Vert \nabla \overline{\rho }(\cdot , (x_i+s\,\widetilde{x}_i)(\cdot ))\Vert _{L^p(0,T)} \\&\qquad + \Vert \overline{\rho }(\tau ,x_i(\tau ))\Vert _{L^2(0,T)} \Vert \nabla \overline{\rho }(\cdot ,(x_i+s\,\widetilde{x}_i)(\cdot )) \\&\qquad - \nabla \overline{\rho }(\cdot ,x_i(\cdot ))\Vert _{L^p(0,T)}\Big ) \Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)} \\&\quad = \textrm{o} ({\Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)}} ), \end{aligned}$$

where the last step follows again from (4.10) and (4.11). This confirms the partial differentiability of $e_5$.

Collecting the previous estimates proves the partial Fréchet differentiability of the operator e. Based on the properties of f, h and $\phi _K$ using again Lipschitz continuity and boundedness, the operators $\partial _\rho e$, $\partial _\phi e$ and $\partial _{x_i} e$ depend continuously on $(\rho ,\phi ,{\varvec{x}})$, which implies the continuous Fréchet differentiability. $\square$

Next, we show that the operator $\partial _{\varvec{y}}e({\varvec{y}},{\varvec{u}})$ is invertible.

Lemma 4.3

The operator $\partial _{\varvec{y}}e({\varvec{y}};{\varvec{u}})$ is bijective. That is, given ${\varvec{y}}= (\rho , \phi , {\varvec{x}})\in \mathcal {Y}$ and $F = (F_1,\ldots ,F_5)^T \in \mathcal {Z}$, the system $\partial _{\varvec{y}}e({\varvec{y}},{\varvec{u}}) \, \widetilde{{\varvec{y}}} = F$ given by

$$\begin{aligned} \begin{aligned} \partial _t \widetilde{\rho } - \varepsilon \,\varDelta \widetilde{\rho } - \nabla \cdot ({\widetilde{\rho }\,\beta ({\varvec{y}}) + \rho \,\frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}}\widetilde{{\varvec{y}}}} )= F_1&\quad \text {in } Q_T , \\ ({\varepsilon \,\nabla \widetilde{\rho } + \widetilde{\rho }\,\beta ({\varvec{y}}) + \rho \,\frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}} \widetilde{{\varvec{y}}}}) \cdot n + {\chi _{\partial \varOmega _O }}\,\eta \,\widetilde{\rho } = F_2&\quad \text {on } \varSigma _T , \\ \widetilde{\rho }(0) = F_3&\quad \text {in } \varOmega , \\ -\delta _1\,\varDelta \widetilde{\phi } + 2 \nabla \phi ^T \nabla \widetilde{\phi } + \frac{2 f(\rho )\,f'(\rho )}{(f^2(\rho )+\delta _2)^2} \widetilde{\rho } = F_4&\quad \text {in } Q_T , \\ \widetilde{x}_i(t) - \int _0^t f'(\overline{\rho }(s,x_i(s))) ({\nabla \overline{\rho }(s,x_i(s))^T \widetilde{x}_i(s) + \overline{\widetilde{\rho }}(s,x_i(s)) })u_i(s) \mathop {}\!d s&= F_{5,i}(t) , \end{aligned} \end{aligned}$$

(4.12)

for $t \in (0,T)$ and $i=1,\ldots ,M$, possesses a unique solution $\widetilde{{\varvec{y}}} = (\widetilde{\rho },\widetilde{\phi },\widetilde{{\varvec{x}}}) \in \mathcal {Y}$.

Proof

The strategy of the proof is to apply Banach’s fixed point theorem to the linear system (4.12) to avoid technicalities that arise from the fact that $\widetilde{{\varvec{x}}}(t)$ depends on $\widetilde{\rho }$ non-locally in time. To this end, we introduce three solution operators.

First, there is $F_\rho :\mathcal {Y}\rightarrow W_p^{2,1}(Q_T)$ which maps $\widetilde{{\varvec{y}}} = (\widetilde{\rho },\widetilde{\phi },\widetilde{{\varvec{x}}})$ to $\widehat{\rho }\in W^{2,1}_p(Q_T)$, which is defined as the solution to

$$\begin{aligned} \begin{aligned} \partial _t \widehat{\rho }- \varepsilon \varDelta \widehat{\rho } - \nabla \cdot ({\widetilde{\rho }\,\beta ({\varvec{y}}) + \rho \,\frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}}\widetilde{{\varvec{y}}}})&= F_1{} & {} \text {in } Q_T , \\ ({\varepsilon \nabla \widehat{\rho } + \widetilde{\rho }\,\beta ({\varvec{y}}) + \rho \,\frac{\partial \beta (\widetilde{{\varvec{y}}})}{\partial {\varvec{y}}}}) \cdot n + {\chi _{\partial \varOmega _O }}\,\eta \,\widetilde{\rho }&= F_2{} & {} \text {on } \varSigma _W , \\ \widehat{\rho }(0)&= F_3{} & {} \text {in } \varOmega . \end{aligned} \end{aligned}$$

(4.13)

Second, we define the operator

$$\begin{aligned} F_\phi :W^{2,1}_p(Q_T) \rightarrow L^\infty (0,T;W_\text {ND}^{2,p}(\varOmega )) \cap H^1(0,T;H^1(\varOmega )), \end{aligned}$$

which maps $\widetilde{\rho }$ to the solution $\widetilde{\phi }$ of

$$\begin{aligned} -\delta _1\varDelta \widetilde{\phi } + 2 \nabla \phi ^T \nabla \widetilde{\phi } + \frac{2 f(\rho )\,f'(\rho )}{(f^2(\rho )+\delta _2)^2} \widetilde{\rho } = F_4\quad \text {in}\ Q_T , \end{aligned}$$

(4.14)

together with the boundary conditions (1.5) (these are incorporated in the function space).

Third, we introduce $F_{\varvec{x}}:W^{2,1}_p(Q_T) \rightarrow W^{1,\infty }([0,T];\mathbb {R}^n)^M$ mapping $\widetilde{\rho }$ to the solution $\widetilde{{\varvec{x}}}=(\widetilde{x}_1,\ldots ,\widetilde{x}_M)^T$ of

$$\begin{aligned} \widetilde{x}_i(t) = \int _0^t f'({\overline{\rho }(s,x_i(s))}) ({\nabla \overline{\rho }(s,x_i(s))^T \widetilde{x}_i(s) + \overline{\widetilde{\rho }}(s,x_i(s))}) u_i(s) \mathop {}\!d s+ F_{5,i}(t) , \end{aligned}$$

for $i=1,\ldots ,M$.

The idea of the proof is to apply Banach’s fixed point theorem in the space $W^{2,1}_p(Q_T)$ to the operator

$$\begin{aligned} \widehat{\rho } = F(\widetilde{\rho }) :=F_\rho (\widetilde{\rho },F_\phi (\widetilde{\rho }),F_{{\varvec{x}}}(\widetilde{\rho })). \end{aligned}$$

(4.15)

As the operator F is affine, it suffices to show the boundedness of the linear part of F by a constant smaller than 1, which will imply that $F(\widetilde{\rho })$ is a contraction.

Step 1: Estimate for $\underline{F_\phi }$. Due to $\phi (t)\in W^{2,p}(\varOmega ) \hookrightarrow W^{1,\infty }(\varOmega )$ for $p>2$ and a.a. $t\in (0,T)$, we confirm that $\nabla \phi \in L^\infty (Q_T)$ holds. Moreover, taking into account $f\in W^{1,\infty }(\mathbb {R})$ and $\rho \in W^{2,p}(Q_T) \hookrightarrow L^\infty (Q_T)$, we get

$$\begin{aligned} \Vert {\frac{2 f(\rho )\,f'(\rho )}{(f^2(\rho )+\delta _2)^2}\,\widetilde{\rho }}\Vert _{L^\infty (0,T;L^p(\varOmega ))} \le 2\,C_f\,C_\rho \,\delta _2^{-2}\,\Vert \widetilde{\rho }\Vert _{L^\infty (0,T;L^p(\varOmega ))} . \end{aligned}$$

From Thm. 3.17 in [49], see also Prop. 2 in [41] for an application to the linearized Eikonal equation, we then deduce that the problem (4.14) possesses for a.a. $t\in (0,T)$ a unique solution $\widetilde{\phi }(t)\in W^{2,p}_\text {ND}(\varOmega )$, which fulfills the inequality

$$\begin{aligned} \Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} \le C({\Vert \widetilde{\rho }\Vert _{L^\infty (0,T;L^p(\varOmega ))} + \Vert F_4\Vert _{L^\infty (0,T;L^p(\varOmega ))}}) . \end{aligned}$$

(4.16)

Furthermore, formal differentiation of (4.14) with respect to time and exploiting $\partial _t F_4\in L^p(0,T;W^{1,p}(\varOmega )^*)$ yields

$$\begin{aligned}{} & {} \Vert \partial _t \widetilde{\phi }\Vert _{L^p(0,T;W^{1,p}(\varOmega ))} \nonumber \\{} & {} \quad \le C \, ({\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} + \Vert F_4\Vert _{L^\infty (0,T;L^p(\varOmega ))} + \Vert \partial _tF_4\Vert _{L^p(0,T;W^{1,p}(\varOmega )^*)}}) .\nonumber \\ \end{aligned}$$

(4.17)

In the second estimate C may depend on T, via norms of $\partial _t\phi$, $\partial _t \rho$, etc.. Since this implies that C is decreasing when T is decreasing, this does not affect the contraction argument.

Step 2: Estimate for $\underline{F_{{\varvec{x}}}}$. Next, we show the boundedness of the solution operator $F_{{\varvec{x}}}$. To this end, we define $C_\rho (t):=\Vert \rho (t,\cdot )\Vert _{W^{2,p}(\varOmega )}$ and deduce from the last equation in (4.12), using the Hölder inequality with $1/p+1/q=1$ and the continuity of ${E }:W^{2,p}(\varOmega ) \hookrightarrow C^{1,\alpha }(\overline{\varOmega })\rightarrow C^{1,\alpha }(\mathbb {R}^2)$, see Lemma 2.4:

$$\begin{aligned} |\widetilde{x}_i(t)|&\le \int _0^t|f'(\overline{\rho }(s,x_i(s)))\left( \nabla \overline{\rho }(s,x_i(s))^T \widetilde{x}_i(s) + \overline{\widetilde{\rho }}(s,x_i(s))\right) u_i(s)| \mathop {}\!d s+ |F_{5,i}(t)| \\&\le C_f\,C_{E }\,C_\infty ({\int _0^t C_\rho (s)\,|\widetilde{x}_i(s)|\mathop {}\!d s+ t^{1/q}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_t)}}) + \Vert F_{5,i}\Vert _{L^\infty (0,t)} . \end{aligned}$$

Due to the Gronwall inequality we obtain

$$\begin{aligned} |\widetilde{x}_i(t)| \le \big (C_f \, C_{E }\, C_\infty t^{1/q}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} + \Vert F_{5,i}\Vert _{L^\infty (0,t)}\big )\, \exp ({C_f\,C_{E }\,C_\infty \int _0^t C_\rho (s)\mathop {}\!d s}) . \end{aligned}$$

We can further estimate $\int _0^t C_\rho (s)\mathop {}\!d s\le t^{1/q}\,\Vert \rho \Vert _{W^{2,1}_p(Q_t)}$. The constants depending only on input data and the linearization point $\rho$ are shifted into the generic constants c, C, which depend on T in monotonically decreasing way only. This implies

$$\begin{aligned} \Vert \widetilde{x}_i\Vert _{L^\infty (0,T)} \le C\,e^{c\,T^{1/q}}({T^{1/q}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} + \Vert F_{5,i}\Vert _{L^\infty (0,T)}}) . \end{aligned}$$

(4.18)

From now on we will choose $T>0$ sufficiently small such that $C\,T^{1/q}\,e^{cT^{1/q}} < 1$.

Furthermore, we can show $L^s$-regularity for the time derivative of $\widetilde{x}_i$. From the classical formulation of the ODE and the usual Hölder arguments we deduce

$$\begin{aligned} |\dot{\widetilde{x}}_i(t)|&= |f'(\overline{\rho }(t,x_i(t))) ({\nabla \overline{\rho }(t,x_i(t))^T \widetilde{x}_i(t) + \overline{\widetilde{\rho }}(t,x_i(t))}) u_i(t) + \dot{F_5}(t)| \\&\le C_f\,C_{E }({C_{\rho }(t)\,|\widetilde{x}_i(t)| + C_{\widetilde{\rho }}(t)}) + |\dot{F}_{5,i}(t)| \end{aligned}$$

for a.a. $t\in [0,T]$. Taking the $L^s(0,T;\mathbb {R}^2)$-norm yields, after insertion of the estimate (4.18) and an application of the Hölder inequality with $1/s = 1/2 + 1/p$,

$$\begin{aligned} \Vert \widetilde{x}_i\Vert _{W^{1,s}(0,T;\mathbb {R}^2)}&\le C_f\,C_{E }\,T^{1/2}({\Vert \rho \Vert _{W^{2,1}_p(Q_T)}\,\Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)} + \Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)}}) \nonumber \\&\quad + \Vert F_{5,i}\Vert _{W^{1,s}(0,T;\mathbb {R}^2)} \nonumber \\&\le C\,T^{1/2}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} + ({1+C\,e^{c\,T^{1/q}}})\Vert F_{5,i}\Vert _{W^{1,s}(0,T;\mathbb {R}^2)}. \end{aligned}$$

(4.19)

Step 3: Estimate for $\underline{F_\rho }$: Next, we consider the solution operator of (4.13). By Lemma 3.2 we have the a priori estimate

$$\begin{aligned} \Vert \widehat{\rho }\Vert _{W^{2,1}_p(Q_T)}&\le C\Bigg [\Vert F_1\Vert _{L^p(Q_T)} + \Vert F_2\Vert _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _W )} + \Vert F_3\Vert _{W^{2(1-1/p),p}(\varOmega )} \nonumber \\&+ \Vert \nabla \cdot ({\widetilde{\rho }\,\beta ({\varvec{y}}) + \rho \frac{\partial \beta ({\varvec{y}})}{\partial y}\widetilde{y}})\Vert _{L^p(Q_T)} \nonumber \\&+ \Vert ({\widetilde{\rho }\,\beta ({\varvec{y}})+\rho \,\frac{\partial \beta ({\varvec{y}})}{\partial y}\widetilde{y}})\cdot n+{\chi _{\partial \varOmega _O }}\,\eta \,\widetilde{\rho }\Vert _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T)} \Bigg ] . \end{aligned}$$

(4.20)

It remains to discuss the last two terms on the right-hand side. First, we confirm with the product rule, the Hölder inequality and Lemma 4.1 that the estimate

$$\begin{aligned} \Vert \nabla \cdot ({\widetilde{\rho }\,\beta ({\varvec{y}})})\Vert _{L^p(Q_T)}&\le C_h\,T^{1/p} ({\Vert \nabla \widetilde{\rho }\,f(\rho ) + \widetilde{\rho }\,f'(\rho )\,\nabla \rho \Vert _{L^\infty (0,T;L^p(\varOmega ))} + {\Vert \widetilde{\rho }\, f(\rho )\Vert _{L^\infty (Q_T)}}}) \nonumber \\&\le C\,T^{1/p}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} \end{aligned}$$

(4.21)

holds. Note that we also exploited regularity results for the linearization point, more precisely, $f(\rho ),\rho \in W^{2,1}_p(Q_T) \hookrightarrow L^\infty (0,T;W^{1,p}(\varOmega )) \hookrightarrow L^\infty (Q_T)$ (see (2.4)), $\phi \in L^\infty (0,T;W^{2,p}(\varOmega ))$, and the assumptions on f, h and $\phi _K$.

Next, we show an estimate for the term involving $\rho \frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}}\widetilde{{\varvec{y}}}$ and we study all partial derivatives of $\beta$ separately. For the derivative with respect to $\rho$ we get

$$\begin{aligned}&\Vert \nabla \cdot ({\rho \frac{\partial \beta ({\varvec{y}})}{\partial \rho }\widetilde{\rho }})\Vert _{L^p(Q_T)} \nonumber \\&\quad = \Vert \nabla \cdot ({\rho \,f'(\rho )\,\widetilde{\rho }\,h(\varPhi )})\Vert _{L^p(Q_T)} \nonumber \\&\quad \le C_h ({\Vert \widetilde{\rho } \, \nabla \rho \, (f'(\rho ) + \rho \,f''(\rho )) + \rho \, \nabla \widetilde{\rho } \, f'(\rho )\Vert _{L^p(Q_T)}} \nonumber \\&\qquad + ){\Vert \rho \,f'(\rho )\,\widetilde{\rho }\Vert _{L^p(0,T;L^\infty (\varOmega ))}} \nonumber \\&\quad \le C \, T^{1/p} \Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} , \end{aligned}$$

(4.22)

where the last step follows from the Hölder inequality, $f\in W^{2,\infty }(\mathbb {R})$, Lemma 4.1 and the embedding $W^{2,1}_p(Q_T) \hookrightarrow L^\infty (0,T;W^{1,p}(\varOmega ))\hookrightarrow L^\infty (Q_T)$, taking into account the regularity of the linearization point $(\rho ,\phi ,{\varvec{x}})$.

For the terms involving $\frac{\partial \beta }{\partial \phi }$ we first confirm by a simple computation

$$\begin{aligned} \frac{\partial h(\varPhi )}{\partial \phi }\widetilde{\phi } = Dh(\varPhi ) \nabla \widetilde{\phi }. \end{aligned}$$

This implies with the product rule

$$\begin{aligned}{} & {} \Vert \nabla \cdot ({\rho \frac{\partial \beta ({\varvec{y}})}{\partial \phi }\widetilde{\phi }})\Vert _{L^p(Q_T)}\nonumber \\{} & {} \quad \le \Vert g'(\rho )\,\nabla \rho \cdot ({Dh(\varPhi )\,\nabla \widetilde{\phi }}) + g(\rho )\,\nabla \cdot ({Dh(\varPhi )\,\nabla \widetilde{\phi }})\Vert _{L^p(Q_T)} \end{aligned}$$

(4.23)

with $g(\rho ) :=\rho \,f(\rho )$ and $g'(\rho ) = \rho \,f'(\rho )+f(\rho )$. Due to $\widetilde{\phi }\in L^\infty (0,T;W^{2,p}(\varOmega ))$ and thus, $\nabla \widetilde{\phi }\in L^\infty (0,T;W^{1,p}(\varOmega )) \hookrightarrow L^\infty (Q_T)$, as well as $\nabla \rho \in L^\infty (0,T;L^p(\varOmega ))$ and Lemma 4.1, we get for the first term on the right-hand side of (4.23)

$$\begin{aligned}&\Vert g'(\rho ) \nabla \rho \cdot ({Dh(\varPhi )\,\nabla \widetilde{\phi }})\Vert _{L^p(Q_T)} \\&\quad \le T^{1/p}\,C_g\,C_h\,\Vert \nabla \rho \Vert _{L^\infty (0,T;L^p(\varOmega ))} \Vert \nabla \widetilde{\phi }\Vert _{L^\infty (Q_T)} \le C\,T^{1/p}\,\Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} . \end{aligned}$$

To bound the second term in (4.23), we insert the identity

$$\begin{aligned} \nabla \cdot ({Dh(\varPhi )\nabla \widetilde{\phi }}) = ({\nabla \cdot Dh(\varPhi )}) \cdot \nabla \widetilde{\phi } + Dh(\varPhi ):\nabla ^2\widetilde{\phi } \end{aligned}$$

and obtain, together with the bounds for Dh, from Lemma 4.1

$$\begin{aligned}&\Vert g(\rho )\,\nabla \cdot ({Dh(\varPhi )\,\nabla \widetilde{\phi }})\Vert _{L^p(Q_T)} \\&\quad \le C_g\,C_h\,T^{1/p}\,({\Vert \nabla \widetilde{\phi }\Vert _{L^\infty (Q_T)} + \Vert \nabla ^2 \widetilde{\phi }\Vert _{L^\infty (0,T;L^p(\varOmega ))}}) \\&\quad \le C\,T^{1/p}\,\Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} . \end{aligned}$$

Insertion of the two previous estimates into (4.23) then yields

$$\begin{aligned} \Vert \nabla \cdot ({\rho \frac{\partial \beta ({\varvec{y}})}{\partial \phi }\widetilde{\phi }})\Vert _{L^p(Q_T)} \le C\,T^{1/p}\,\Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} . \end{aligned}$$

(4.24)

The derivative with respect to $x_i$, $i=1,\ldots ,M$, is handled using the chain rule, Lemma 4.1 and assumption (K2), leading to

$$\begin{aligned} \Vert \nabla \cdot ({\rho \frac{\partial \beta ({\varvec{y}})}{\partial x_i}\widetilde{x}_i})\Vert _{L^p(Q_T)}&= \Vert \nabla \cdot ({g(\rho )\,Dh(\varPhi )\,\nabla ^2K(\cdot -x_i)\,\widetilde{x}_i})\Vert _{L^p(Q_T)} \nonumber \\&\le C\,T^{1/p}\,\Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)} . \end{aligned}$$

(4.25)

The estimates (4.21), (4.22), (4.24) and (4.25) give an estimate for the fourth term in (4.20), namely

$$\begin{aligned} \Vert \nabla \cdot ({\widetilde{\rho }\,\beta ({\varvec{y}}) + \rho \frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}}\widetilde{{\varvec{y}}}})\Vert _{L^p(Q_T)} \le C\,T^{1/p}\,\Vert \widetilde{{\varvec{y}}}\Vert _\mathcal {Y}. \end{aligned}$$

(4.26)

It remains to show the boundedness of the fifth term in (4.20) and we apply a trace theorem as in the proof of Lemma 4.2, see (4.7). This requires us to show appropriate bounds in the $L^p(0,T;W^{1,p}(\varOmega ))$- and $W^{1,s}(0,T;L^p(\varOmega ))$-norms.

First, when replacing $\nabla \cdot$ by an arbitrary partial derivative $\partial _{x_j}$ in the previous considerations, we get analogously to (4.26)

$$\begin{aligned} \Vert \widetilde{\rho }\,\beta ({\varvec{y}}) + \rho \,\frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}} \widetilde{{\varvec{y}}}\Vert _{L^p(0,T;W^{1,p}(\varOmega ))} \le C\,T^{1/p}\,\Vert \widetilde{{\varvec{y}}}\Vert _\mathcal {Y}. \end{aligned}$$

(4.27)

Furthermore, we study the temporal derivative of the expressions occurring in the fifth term of (4.20). The Hölder inequality and Lemma 4.1 yield

$$\begin{aligned}&\Vert \partial _t (\widetilde{\rho }\,\beta ({\varvec{y}}))\Vert _{L^s(0,T;L^p(\varOmega ))} \\&\quad \le T^{1/2}\,\Vert \partial _t ({\widetilde{\rho }\,f(\rho )\, h(\varPhi )})\Vert _{L^p(Q_T)} \\&\quad \le C_f\,C_h\,T^{1/2} ({\Vert \partial _t\widetilde{\rho }\Vert _{L^p(Q_T)}+\Vert \widetilde{\rho }\Vert _{L^\infty (Q_T)}\,(\Vert \partial _t\rho \Vert _{L^p(Q_T)} + 1)}) \\&\quad \le C_f\,C_h\,C_\infty \,(C_\rho +1)\,T^{1/2}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} . \end{aligned}$$

For the term involving $\frac{\partial \beta ({\varvec{y}})}{\partial \rho } \widetilde{\rho }$ we get

$$\begin{aligned}&\left\| \partial _t (\rho \,\frac{\partial \beta ({\varvec{y}})}{\partial \rho }\widetilde{\rho })\right\| _{L^s(0,T;L^p(\varOmega ))} \\&\quad \le T^{1/2}\,\Vert \partial _t (\rho \,f'(\rho )\,\widetilde{\rho }\,h(\varPhi ))\Vert _{L^p(Q_T)}\\&\quad \le C_f\,C_h\,T^{1/2}\,({2\,\Vert \partial _t \rho \Vert _{L^p(Q_T)}\,\Vert \widetilde{\rho }\Vert _{L^\infty (Q_T)}}\\&\qquad + ){\Vert \rho \Vert _{L^\infty (Q_T)}\,(\Vert \partial _t \widetilde{\rho }\Vert _{L^p(Q_T)} +\Vert \widetilde{\rho }\Vert _{L^\infty (Q_T)})} \\&\quad \le C_f\,C_h\,C_\rho \,C_\infty \,T^{1/2}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} . \end{aligned}$$

For the term depending on $\frac{\partial \beta ({\varvec{y}})}{\partial \phi }\widetilde{\phi }$ we obtain

$$\begin{aligned}&\left\| \partial _t ({\rho \,\frac{\partial \beta ({\varvec{y}})}{\partial \phi }\widetilde{\phi }})\right\| _{L^s(0,T;L^p(\varOmega ))} \le T^{1/2}\,\Vert \partial _t(g(\rho )\,Dh(\varPhi )\,\nabla \widetilde{\phi })\Vert _{L^p(Q_T)} \\&\quad \le C_g\,C_h\,T^{1/2}\,({({\Vert \partial _t \rho \Vert _{L^p(Q_T)}+1})\Vert \nabla \widetilde{\phi }\Vert _{L^\infty (Q_T)} + \Vert \partial _t\nabla \widetilde{\phi }\Vert _{L^p(Q_T)}}) \\&\quad \le C_g\,C_h(2+C_\rho )\,T^{1/2}\,({\Vert \widetilde{\phi }\Vert _{L^\infty (0,T;W^{2,p}(\varOmega ))} + \Vert \partial _t \widetilde{\phi }\Vert _{L^p(0,T;W^{1,p}(\varOmega ))}}) . \end{aligned}$$

Finally, for the derivative with respect to $x_i$ we derive

$$\begin{aligned}&\left\| \partial _t ({\rho \,\frac{\partial \beta ({\varvec{y}})}{\partial x_i}\widetilde{x}_i})\right\| _{L^s(0,T;L^p(\varOmega ))} \\&\quad = \Vert \partial _t(g(\rho )\,Dh(\varPhi )\,\nabla ^2K(\cdot -x_i)\,\widetilde{x}_i)\Vert _{L^s(0,T;L^p(\varOmega ))} \\&\quad \le C_g\,C_h\,C_K\,C_{\varvec{x}}({(C_\rho +1)\,T^{1/2}\,\Vert \widetilde{x}_i\Vert _{L^\infty (0,T;\mathbb {R}^2)} + 2\,\Vert \widetilde{x}_i\Vert _{W^{1,s}(0,T;\mathbb {R}^2)}}) . \end{aligned}$$

The previous four estimates and an embedding yield the desired property

$$\begin{aligned} \left\| \widetilde{\rho }\,\beta ({\varvec{y}})+\rho \,\frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}}\widetilde{{\varvec{y}}}\right\| _{W^{1/2,p}(0,T;L^p(\varOmega ))}&\le C\,\left\| \widetilde{\rho }\,\beta ({\varvec{y}})+\rho \,\frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}}\widetilde{{\varvec{y}}}\right\| _{W^{1,s}(0,T;L^p(\varOmega ))} \nonumber \\&\le C\,T^{1/2}({\Vert \widetilde{\rho }\Vert _{\mathcal {Y}_1} + \Vert \widetilde{\phi }\Vert _{\mathcal {Y}_2}}) + C\,\Vert \widetilde{{\varvec{x}}}\Vert _{\mathcal {Y}_3} . \end{aligned}$$

(4.28)

Next, we may combine the estimates (4.28) and (4.27) and apply the trace Lemma 2.3 to arrive at

$$\begin{aligned}&\left\| ({\widetilde{\rho }\,\beta ({\varvec{y}})+\rho \,\frac{\partial \beta ({\varvec{y}})}{\partial {\varvec{y}}}\widetilde{{\varvec{y}}}})\cdot n+{\chi _{\partial \varOmega _O }}\,\eta \,\widetilde{\rho }\right\| _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _T)} \nonumber \\&\quad \le C\,T^{1/p}({\Vert \widetilde{\rho }\Vert _{\mathcal {Y}_1} + \Vert \widetilde{\phi }\Vert _{\mathcal {Y}_2}}) + C\,\Vert \widetilde{{\varvec{x}}}\Vert _{\mathcal {Y}_3} . \end{aligned}$$

(4.29)

From the estimates (4.26) and (4.29) we deduce that there exists a unique solution $\widehat{\rho }\in W^{2,1}_p(Q_T)$ of (4.13). Moreover, together with (4.20) using $p>2$, we obtain for sufficiently small $T>0$ the a priori estimate

$$\begin{aligned}&\Vert \widehat{\rho }\Vert _{W^{2,1}_p(Q_T)} \nonumber \\&\quad \le C\,\Big (T^{1/p}({\Vert \widetilde{\rho }\Vert _{\mathcal {Y}_1} + \Vert \widetilde{\phi }\Vert _{\mathcal {Y}_2}}) + \Vert \widetilde{{\varvec{x}}}\Vert _{\mathcal {Y}_3} \nonumber \\&\qquad + \Vert F_1\Vert _{L^p(Q_T)} + \Vert F_2\Vert _{W^{1-1/p,1/2-1/(2p)}_p(\varSigma _W )} + \Vert F_3\Vert _{W^{2(1-1/p),p}(\varOmega )}\Big ) . \end{aligned}$$

(4.30)

Step 4: The fixed-point argument

Now, we can combine the estimates (4.16), (4.17), (4.19) and (4.30) to deduce the continuity of the operator $F(\widetilde{\rho })=F_\rho (\widetilde{\rho },F_\phi (\widetilde{\rho }),F_{\varvec{x}}(\widetilde{\rho }))$. This gives the estimate

$$\begin{aligned} \Vert F(\widetilde{\rho })\Vert _{W^{2,1}_p(Q_T)} \le C\,\Vert F\Vert _\mathcal {Z}+ C\,T^{1/p}\,\Vert \widetilde{\rho }\Vert _{W^{2,1}_p(Q_T)} \end{aligned}$$

for $T \le 1$.

As F is an affine linear operator we directly conclude the contraction property provided that T is sufficiently small. Thus, Banach’s fixed-point theorem provides the existence of a unique solution $\widetilde{\rho }\in W^{2,1}_p(Q_T)$ which can be extended to an arbitrary time-scale with a concatenation argument.

The regularities of the corresponding potential $\widetilde{\phi }$ and the agent trajectories $\widetilde{x}_i$, $i=1,\ldots ,M$, follow from (4.16), (4.17) and (4.19).

The previous two lemmas, together with the implicit function theorem, imply the differentiability of the control-to-state operator. The derivative in a direction $\widetilde{{\varvec{u}}}\in L^\infty (0,T;\mathbb {R}^2)^M$ can be computed by means of

$$\begin{aligned} e_{\varvec{y}}({\varvec{y}},{\varvec{u}}) \, \widetilde{{\varvec{y}}} = - e_{\varvec{u}}({\varvec{y}},{\varvec{u}}) \, \widetilde{{\varvec{u}}} . \end{aligned}$$

We summarize the final result in the following theorem.

Theorem 4.4

The control-to-state operator $S :\mathcal {U}\rightarrow \mathcal {Y}$ is Fréchet-differentiable and $\widetilde{{\varvec{y}}} = S'({\varvec{u}}) \, \widetilde{{\varvec{u}}}$ for given $\widetilde{{\varvec{u}}} \in \mathcal {U}$ is characterized by the unique solution of

$$\begin{aligned} \begin{aligned}&\partial _t \widetilde{\rho } - \varepsilon \varDelta \widetilde{\rho }\\&\quad - \nabla \cdot \left( \widetilde{\rho } \,\beta (\rho ,\phi ,{\varvec{x}}) + \rho \left( \frac{\partial \beta (\rho ,\phi ,{\varvec{x}})}{\partial \rho } \widetilde{\rho } + \frac{\partial \beta (\rho ,\phi ,{\varvec{x}})}{\partial \phi } \widetilde{\phi } + \frac{\partial \beta (\rho ,\phi ,{\varvec{x}})}{\partial {\varvec{x}}} \widetilde{{\varvec{x}}} \right) \right) = 0 ,\\&\quad - \delta _1 \varDelta \widetilde{\phi } + 2\nabla \phi \cdot \nabla \widetilde{\phi } + \frac{2 f(\rho )\,f'(\rho )}{(f^2(\rho )+\delta _2)^2} \widetilde{\rho } = 0 ,\\&\quad \dot{\widetilde{x}}_i - f'(\overline{\rho }(\cdot ,x_i)) \left( \nabla \overline{\rho }(\cdot ,x_i)^T \widetilde{x}_i + \overline{\widetilde{\rho }}(\cdot ,x_i) \right) u_i = f(\overline{\rho }(\cdot ,x_i)) \, \widetilde{u}_i, \end{aligned} \end{aligned}$$

(4.31)

for $i=1,\ldots ,M$, together with the boundary conditions (1.5) and homogeneous initial conditions

$$\begin{aligned} \widetilde{\rho }(\cdot ,0) = 0 \quad \text {and} \quad \widetilde{x}_i(0) = 0, \quad i=1,\ldots ,M. \end{aligned}$$

5 Optimal Control Problem

In this section we come back to the optimal control problem outlined in Sect. 1.1. We consider objective functionals of the form

$$\begin{aligned} J(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) :=\varPsi (\rho ,\phi ,{\varvec{x}}) + \frac{\alpha }{2\,T} \Vert {\varvec{u}}\Vert _{H^1(0,T;\mathbb {R}^2)^M}^2, \end{aligned}$$

(5.1)

where $\alpha >0$ is a regularization parameter and the functional $\varPsi :\mathcal {Y}\rightarrow \mathbb {R}$ fulfills the following assumption:

• (J1) The functional $\varPsi :\mathcal {Y}\rightarrow \mathbb {R}$ is weakly lower semi-continuous and bounded from below on $\{(\rho ,\phi ,{\varvec{x}})\in \mathcal {Y}:\rho \ge 0\ \text {a.\,e.\ in}\ Q_T\}$.

The assumed weak lower semi-continuity is fulfilled, e. g., when $\varPsi$ is convex and continuous. As the density part $\rho$ of each solution of the forward model (1.9) is non-negative it suffices to assume boundedness of $\varPsi$ on a subset only.

We consider the following optimal control problem:

$$\begin{aligned} \left\{ \begin{aligned} \text {Minimize} \quad&J(\rho ,\phi ,{\varvec{x}};{\varvec{u}}) \quad \text {where } (\rho ,\phi ,{\varvec{x}}; {\varvec{u}}) \in \mathcal {Y}\times \mathcal {U}\\ \text {subject to} \quad&(\rho ,\phi ,{\varvec{x}}) = S({\varvec{u}}) \\ \text {and} \quad&{\varvec{u}}\in \mathcal {U}_{\text {ad}}. \end{aligned} \right. \end{aligned}$$

(5.2)

The set of admissible controls is defined by

$$\begin{aligned} \mathcal {U}_{\text {ad}}:={\{}{{\varvec{u}}\in \mathcal {U}}{|}{\Vert u_i\Vert _{L^\infty (0,T)} \le 1 \text { for } i=1,\ldots ,M}{\}} . \end{aligned}$$

This general framework covers in particular the two applications mentioned in Sect. 1

With standard arguments, see, e. g., Ch. 4.4 in [50], (J1) and Lemma 3.9 imply the following existence result for (5.2).

Theorem 5.1

Let assumptions (A1)–(A4), (K1)–(K2), (C1) and (J1) hold. Then problem (5.2) possesses at least one global solution ${\varvec{u}}^*\in \mathcal {U}_{\text {ad}}$.

Moreover, with the properties shown for S in Sects. 3 and 4, in particular the Fréchet differentiability, we may deduce the following first-order necessary optimality condition, which is a direct consequence of the chain rule applied to the reduced optimization problem with objective $j({\varvec{u}}) :=J(S({\varvec{u}});{\varvec{u}})$ and subject to ${\varvec{u}}\in \mathcal {U}_{\text {ad}}$:

Theorem 5.2

Let assumptions (A1)–(A4), (K1)–(K2), (C1) and (J1) hold. Then each local minimizer $({\varvec{y}}^*,{\varvec{u}}^*)\in \mathcal {Y}\times \mathcal {U}_{\text {ad}}$ of (5.2) fulfills

$$\begin{aligned} \langle {\varPsi '({\varvec{y}}^*)\,,\delta {\varvec{y}}}\rangle _{\mathcal {Y}^* \times \mathcal {Y}} + \frac{\alpha }{T} \, ({{\varvec{u}}^*\,,\,\delta {\varvec{u}}})_{H^1(0,T)} \ge 0 \quad \text {for all } \delta {\varvec{u}}\in \mathcal {T}_{\mathcal {U}_{\text {ad}}}{({{\varvec{u}}^*})}, \end{aligned}$$

where ${\varvec{y}}^* = S({\varvec{u}}^*)$ and $\delta {\varvec{y}}= S'({\varvec{u}}^*) \, \delta {\varvec{u}}$ and $S'$ is characterized by the system from Theorem 4.4. Moreover, $\mathcal {T}_{\mathcal {U}_{\text {ad}}}{({{\varvec{u}}^*})}$ denotes the tangent cone to $\mathcal {U}_{\text {ad}}$ at ${\varvec{u}}^*$.

The optimality conditions from Theorem 5.2 can be rewritten by expressing $\langle {\varPsi '({\varvec{y}}^*)\,,\delta {\varvec{y}}}\rangle _{\mathcal {Y}^* \times \mathcal {Y}}$ in terms of a suitably defined adjoint state. This then serves as the starting point for gradient-based optimization methods. Such methods, as well as appropriate discretization of the forward and adjoint systems and numerical results will be discussed in a forthcoming publication.

Auxiliary Results

First we present the proof of Lemma 3.3. It is based on choosing derivatives of the solution as test functions and a Galerkin approximation as well as the uniqueness of weak solutions, see also Ch. III, Sect. 6 in [42].

Proof of Lemma 3.3

Existence in $L^2(0,T;H^1(\varOmega )) \cap H^1(0,T;H^1(\varOmega )^*)$ follows by a standard Galerkin approximation. To this end we introduce $\rho ^N(t,x) = \sum _{i=0}^N d_i(t)\varphi _i(x)$, where $(\varphi _i)_{i\in \mathbb {N}}$ is an orthonormal basis in $H^1(\varOmega )$ which is also orthogonal w.r.t. $L^2(\varOmega )$. The coefficients satisfy $d_i(0) = ({\rho _0\,,\,\varphi _i})_{H^1(\varOmega )}$. Choosing $\rho ^N$ as a test function results in

$$\begin{aligned}&\int _0^T\frac{\mathop {}\!d }{\mathop {}\!d t} \int _\varOmega (\rho ^N)^2\mathop {}\!d t\mathop {}\!d x+ \varepsilon \int _{Q_T} |\nabla \rho ^N|^2\mathop {}\!d t\mathop {}\!d x\nonumber \\&\quad \le \int _{Q_T} |\nabla \rho ^N||g(\rho ^N)\,h|\mathop {}\!d t\mathop {}\!d x-\eta \int _{{\varSigma _O }} (\rho ^N)^2 \mathop {}\!d t\mathop {}\!d s_x \nonumber \\&\quad \le \frac{\varepsilon }{2} \int _{Q_T} |\nabla \rho ^N|^2\mathop {}\!d t\mathop {}\!d x+ \frac{1}{2\varepsilon }\int _{Q_T} |g(\rho ^N) \, h|^2 \mathop {}\!d t\mathop {}\!d x, \end{aligned}$$

(A.1)

i. e.

$$\begin{aligned} \Vert \rho ^N\Vert _{L^\infty (0,T;L^2(\varOmega ))} + \Vert \nabla \rho ^N\Vert _{L^2(Q_T)} \le C ({\Vert h\Vert _{L^2(0,T;L^2(\varOmega ))} + \Vert \rho _0\Vert _{L^2(\varOmega )}}), \end{aligned}$$

(A.2)

where we used $\Vert \rho ^N(\cdot ,0)\Vert _{L^2(\varOmega )} \le C \, \Vert \rho _0\Vert _{L^2(\varOmega )}$. This readily implies that $\rho ^N$ is bounded, uniformly in N, in

$$\begin{aligned} L^2(0,T;H^1(\varOmega )) \cap H^1(0,T;H^1(\varOmega )^*) \cap C([0,T];L^2(\varOmega )) , \end{aligned}$$

see, e. g., Ch. 7.1.2, Thm. 2 in [29]. A subsequent application of the Aubin-Lions Lemma (Thm. II.5.16 in [10]), ensures the convergence of $\rho ^N$ strongly in $L^2(Q_T)$ to $\rho$ and thus the convergence of the nonlinear term

$$\begin{aligned} \int _{Q_T} g(\rho ^N)\,h\cdot \nabla \xi \mathop {}\!d t\mathop {}\!d x\rightarrow \int _{Q_T} g(\rho )\,h\cdot \nabla \xi \mathop {}\!d t\mathop {}\!d x\end{aligned}$$

as $N\rightarrow \infty$ which follows by dominated convergence, noting that $g(\rho ^N)\,h\cdot \nabla \xi$ is uniformly bounded in $L^1(0,T;L^1(\varOmega ))$. Hence the limit satisfies the weak formulation (3.10).

To show uniqueness, assume that there are two solutions $\rho _1$ and $\rho _2$. Testing each equation with $\rho _1 - \rho _2$, integration w.r.t. to t, and using the Lipschitz continuity of g together with an application of the weighted Cauchy inequality yields

$$\begin{aligned}&\frac{1}{2} \int _\varOmega (\rho _1(t,\cdot )-\rho _2(t,\cdot ))^2\mathop {}\!d x+ \frac{\varepsilon }{2} \int _{Q_t} |\nabla (\rho _1-\rho _2)|^2 \mathop {}\!d t\mathop {}\!d x\\&\quad \le \frac{1}{2\varepsilon }\Vert h\Vert _{L^\infty (0,t;L^\infty (\varOmega ))} \int _{Q_t} (\rho _1- \rho _2)^2\mathop {}\!d t\mathop {}\!d x- \eta \int _{\varSigma _{D,t}} (\rho _1-\rho _2)^2 \mathop {}\!d t\mathop {}\!d s_x \end{aligned}$$

for a.a. $t \in (0,T)$. As the last term on the right-hand side is non-positive, the uniqueness follows from Gronwall’s inequality as both $\rho _1$ and $\rho _2$ have the same initial datum.

To show the additional regularity, we choose, for finite N, the test function $\partial _t \rho ^N$ and integrate w.r.t. time. We obtain

$$\begin{aligned} \int _{Q_t} (\partial _t\rho ^N)^2 \mathop {}\!d t\mathop {}\!d x&+ \varepsilon \int _{Q_t} \nabla \rho ^N\cdot \nabla (\partial _t\rho ^N)\mathop {}\!d t\mathop {}\!d x- \int _{Q_t} g(\rho ^N)\,h\cdot \nabla (\partial _t \rho ^N)\mathop {}\!d t\mathop {}\!d x\end{aligned}$$

(A.3)

$$\begin{aligned}&\qquad = - \eta \int _{\varSigma _{D,t}} \rho _N\,\partial _t\rho _N \mathop {}\!d t\mathop {}\!d s_x . \end{aligned}$$

(A.4)

For the second term on the left-hand side we have

$$\begin{aligned} \varepsilon \int _{Q_t} \nabla \rho ^N\cdot \nabla (\partial _t\rho ^N)\mathop {}\!d t\mathop {}\!d x&= \frac{\varepsilon }{2}\int _0^t\frac{\mathop {}\!d }{\mathop {}\!d t}\int _\varOmega |\nabla \rho ^N|^2\mathop {}\!d t\mathop {}\!d x\\&= \frac{\varepsilon }{2}\int _\varOmega |\nabla \rho ^N(t,\cdot )|^2\mathop {}\!d x- \frac{\varepsilon }{2}\int _\varOmega |\nabla \rho ^N(0,x)|^2\mathop {}\!d x, \end{aligned}$$

while for the third one, integration by parts in time together with the (weighted) Cauchy’s inequality gives

$$\begin{aligned}&-\int _{Q_t} g(\rho ^N)\,h\cdot \nabla (\partial _t \rho ^N)\mathop {}\!d t\mathop {}\!d x\\&\quad = \int _{Q_t} \partial _t \rho ^N g'(\rho ^N)\, h\cdot \nabla \rho ^N\mathop {}\!d t\mathop {}\!d x\\ {}&\quad + \int _{Q_t} g(\rho ^N)\,\partial _t h\cdot \nabla \rho ^N\mathop {}\!d t\mathop {}\!d x- \int _{\varOmega } |{g(\rho ^N)\,h \cdot \nabla \rho ^N\mathop {}\!d x}^t_0 \\&\quad \le \frac{1}{2}\Vert \partial _t\rho ^N\Vert _{L^2(0,T;L^2(\varOmega ))}^2 + \frac{1}{2}({\Vert g'(\rho ^N)h\Vert _{L^\infty (0,T;L^\infty (\varOmega ))}^2 + 1}) \Vert \nabla \rho ^N\Vert _{L^2(0,T;L^2(\varOmega ))}^2 \\&\qquad + \Vert g(\rho ^N)\,\partial _t h\Vert _{L^2(0,T;L^2(\varOmega ))}^2 + \frac{\varepsilon }{4}\Vert \nabla \rho ^N(t,\cdot )\Vert _{L^2(\varOmega )}^2+ \frac{1}{4\varepsilon }\Vert g(\rho ^N)\,h\Vert _{L^\infty (0,T;L^2(\varOmega ))}^2 \\&\qquad + \Vert g(\rho (0,\cdot )^N)\,h(0,\cdot ) \, \nabla \rho ^N(0,\cdot )\Vert _{L^2(\varOmega )}^2 . \end{aligned}$$

The term on the right-hand side of (A.3) becomes

$$\begin{aligned} - \eta \int _{\varSigma _{D,t}} \rho _N \partial _t\rho _N \mathop {}\!d t\mathop {}\!d s_x = -\frac{\eta }{2} \int _\varOmega (\rho ^N(t,x))^2 \mathop {}\!d x+ \frac{\eta }{2} \int _\varOmega (\rho ^N(x,0))^2 \mathop {}\!d x. \end{aligned}$$

Combining all these estimates, an application of the trace theorem applied to $\Vert \rho _0\Vert _{L^2({\partial \varOmega _O })}$ and using

$$\begin{aligned} \Vert \rho ^N(\cdot ,0)\Vert _{H^1(\varOmega )} \le \Vert \rho _0\Vert _{H^1(\varOmega )} \end{aligned}$$

yields

$$\begin{aligned}{} & {} \frac{1}{2}\Vert \partial _t \rho ^N\Vert _{L^2(0,T;L^2(\varOmega ))}^2 + \frac{\varepsilon }{4}\Vert \nabla \rho ^N\Vert _{L^\infty (0,T;L^2(\varOmega ))}^2\nonumber \\{} & {} \le C\Vert \rho _0^N\Vert _{H^1(\varOmega )}^2 + c \, \Vert \nabla \rho ^N\Vert _{L^2(0,T;L^2(\varOmega ))}^2 + d , \end{aligned}$$

(A.5)

where

$$\begin{aligned} c&= \frac{1}{2} ({\Vert g'(\rho ^N)\,h\Vert _{L^\infty (0,T;L^\infty (\varOmega ))}^2 + 1}) , \\ C&= \Vert g(\rho (0,\cdot )^N)\,h(0,\cdot )\Vert _{L^\infty (\varOmega )}^2, \\ d&= \frac{1}{2}\Vert g(\rho ^N)\,\partial _t h\Vert _{L^2(0,T;L^2(\varOmega ))}^2 + \frac{1}{4\varepsilon }\Vert g(\rho ^N)\,h\Vert _{L^\infty (0,T;L^2(\varOmega ))}^2. \end{aligned}$$

Using (A.2), we can estimate the second term on the right-hand side which yields the desired regularity for $\rho ^N$. We also infer the existence of a weakly converging subsequence in $L^2(0,T;L^2(\varOmega ))\cap L^\infty (0,T;H^1(\varOmega ))$ which, by the weak lower semi-continuity of the norms, yields the bound also for the limit $\rho$. Reinserting this into (A.5) yields the assertion.

We also state the definition of Carathéodory conditions.

Definition A.1

We say that a function $r:Q_T \rightarrow \mathbb {R}^2$ satisfies the Carathéodory conditions whenever

• for each fixed $x \in \mathbb {R}^d$, the function $t \mapsto r(t,x)$ is measurable,

• for a.a. $t \in [0,T]$, the function $x \mapsto r(t,x)$ is continuous,

• there exists an integrable function m s. t.

  $$\begin{aligned} |r(t,x)| \le m(t) \text { for all } x \in \varOmega \text { and a.a.} t \in [0,T] . \end{aligned}$$

Lemma A.2

Recall that $\varphi _\gamma$ is a standard mollifier. For given $\rho \in L^\infty (0,T;L^1(\varOmega ))$, the function $\rho _\gamma :={\varphi _\gamma } *\rho$ is Lipschitz continuous in x with the Lipschitz constant depending on $\gamma$ and $\Vert \rho \Vert _{L^1(\varOmega )} \in L^\infty (0,T)$, only.

Proof

By definition, $\varphi _\gamma$ is Lipschitz continuous and we denote its Lipschitz constant by $L_\gamma$. Then, for a.a. $t \in (0,T)$, we have

$$\begin{aligned} |\rho _\gamma (t,x) - \rho _\gamma (t,y)|&\le \int _\varOmega |{\varphi _\gamma }(x-z) - {\varphi _\gamma }(y-z)||\rho (t,z)| \mathop {}\!d z\\&\le L_\gamma \int _\varOmega |x-y||\rho (t,z)| \mathop {}\!d z\le L_\gamma \Vert \rho (t,\cdot )\Vert _{L^1(\varOmega )}|x-y| . \end{aligned}$$

Lemma A.3

Fix $2 \le p < \infty$ and denote by $(\rho _{\gamma _k})_k \in W^{2,1}_p(Q_T)$ a bounded family of functions, i. e. there exists a constant $C_\rho$, independent of $\gamma$, s. t.

$$\begin{aligned} \Vert \rho _{\gamma _k}\Vert _{W^{2,1}_p(Q_T)} \le C_\rho . \end{aligned}$$

Denoting by $\phi _{\gamma _k}$ the sequence of solutions to (1.9b) with $\rho _{\gamma _k}$ on the right-hand side, there exists an element $\phi \in L^p(0,T;W^{1,p}(\varOmega ))$ and a subsequence, again denoted by $\phi _{\gamma _k}$, s. t.

$$\begin{aligned} \phi _{\gamma _k} \rightarrow \phi \text { in } L^p(0,T;W^{1,p}(\varOmega )). \end{aligned}$$

(A.6)

Proof

Instead of (1.9b), we work with the transformed version (3.6) and denote its solution by $\psi _\gamma$. We seek to apply a variant of the Aubin-Lions Lemma, introduced in [48], where instead of time derivatives we have to control finite differences in time. Indeed, as $W^{2,p}(\varOmega ) \hookrightarrow \hookrightarrow W^{1,p}(\varOmega )$, it suffices to show

$$\begin{aligned} \int _0^{T-h} \Vert \psi _\gamma (t+h,\cdot ) - \psi _\gamma (t,\cdot )\Vert _{W^{1,p}(\varOmega )}^p \mathop {}\!d t\le \mathcal {O}(h) . \end{aligned}$$

To this end, note that the difference $\psi _\gamma (t+h,\cdot ) - \psi _\gamma (t,\cdot )$ satisfies the equation

$$\begin{aligned} \varDelta (\psi _\gamma (t+h,\cdot ) - \psi _\gamma (t,\cdot )) + q_t (\psi _\gamma (t+h,\cdot ) - \psi _\gamma (t,\cdot )) = (1-\psi (t,\cdot ))(q_{t+h} - q_t), \end{aligned}$$

with boundary conditions (3.6) and the corresponding a priori estimate

$$\begin{aligned} \Vert \psi _\gamma (t+h,\cdot ) - \psi _\gamma (t,\cdot )\Vert _{W^{1,p}(\varOmega )} \le C(C_\phi )\Vert q_{t+h} - q_t\Vert _{L^p(\varOmega )} . \end{aligned}$$

This yields

$$\begin{aligned} \int _0^{T-h} \Vert \psi _\gamma (t+h,\cdot ) - \psi _\gamma (t,\cdot )\Vert _{W^{1,p}(\varOmega )}^p \mathop {}\!d t\le C\int _0^{T-h} \Vert q_{t+h}(\cdot )- q_t(\cdot )\Vert _{W^{1,p}(\varOmega )}^p \mathop {}\!d t. \end{aligned}$$

(A.7)

To estimate the right-hand side, note that $\rho _\gamma \in W^{1,p}(0,T;L^p(\varOmega )) \hookrightarrow C^{0,1/2}([0,T]; L^p(\varOmega ))$, and thus there exists a constant $L_\rho$, depending only on $C_\phi$ but not on $\gamma$, s. t. for all $t \ge 0$ and $h > 0$ satisfying $t+h \le T$, the inequality

$$\begin{aligned} \Vert \rho _\gamma (t+h,\cdot ) - \rho _\gamma (t,\cdot )\Vert _{L^p(\varOmega )} \le L_\rho \sqrt{h} \end{aligned}$$

holds. The definition of $q_t$ in (3.4) then directly implies

$$\begin{aligned} \Vert q_{t+h}(\cdot ) - q_t(\cdot )\Vert _{L^p(\varOmega )} \le CL_\rho \sqrt{h} , \end{aligned}$$

which, when inserted into the right-hand side of (A.7), yields

$$\begin{aligned} \int _0^{T-h} \Vert \psi _\gamma (t+h,\cdot ) - \psi _\gamma (t,\cdot )\Vert _{W^{1,p}(\varOmega )}^p \mathop {}\!d t\le (CL_\rho )^p \, (T-h) \, h^{p/2} . \end{aligned}$$

Since $2 \le p$ holds, this is the desired estimate. Thus, Lemma on p. 1011 in [48] ensures the existence of a subsequence that strongly converges in $L^p(0,T;W^{1,p}(\varOmega ))$. By inverting the transformation (3.5), the result for $\phi _\gamma$ is obtained. $\square$