$\int_{\Omega}u^{\varepsilon}(t)m_{1}\leq C(T-t)$

we also deduce, letting first $\varepsilon\to 0$ and then $t\to T$, that
$\int_{\Omega}u(T)m_{1}\leq 0$. As for $m$, if $w$ denotes the weak limit
(e.g. in $L^{2}(Q)$) of $m^{\varepsilon}H_{p}(Du^{\varepsilon})$, then $m$ is
a bounded weak solution of the continuity equation

$\begin{cases}m_{t}-\operatorname{div}(w)=0&\\\ m\in
C([0,T],{\mathcal{P}}(\Omega))\,,\quad
m(0)=m_{0}\,,\,m(T)=m_{1}\,.&\end{cases}$ (6.24)

The duality between (6.23) and (6.24) was exploited in [31, Section 5.3];
roughly speaking, it is possible to cross multiply (6.23) and (6.24). Using
that $\int_{\Omega}u(T)m_{1}\leq 0$, and reasoning as in [9, Thm 1.15], then
one shows

$\limsup_{\varepsilon\to
0}\int_{0}^{T}\int_{\Omega}f(m^{\varepsilon})m^{\varepsilon}\leq\int_{\Omega}\bar{f}\,m\,dxdt\,\,.$

A standard Minty-type argument implies that $\bar{f}=f(m)$, and since $f$ is
increasing, one also deduces that $m^{\varepsilon}\to m$ a.e. in $Q$, and
therefore in $L^{p}(Q)$ for every $p<\infty$. Moreover, the identification of
$\bar{f}$ goes together with the identification of $w=m\,H_{p}(Du)$, using the
convexity of $(m,w)\mapsto L(w/m)m$. Hence one follows arguments in [31] or
[9]: $(u,m)$ satisfy the equality

$\int_{0}^{T}\int_{\Omega}m\left(H_{p}(Du)Du-H(Du)\right)+(f(m)+V)m=\int_{\Omega}m_{0}\,u(0)\,dx$

and $\int_{\Omega}u(T)m_{1}\,dx$ is proved to be zero. It follows from the
above equality that $m$ is a minimum of the limit functional and that $(u,m)$
is a weak solution of (6.20) in the sense of Definition 6.7. The uniqueness of
$m$ and of $u$ up to $m$-negligible sets can be proved as in [9, Thm 1.16].
$\square$

#### The compact case.

All results proved here holds in more generality if the state space is the
flat torus ${\mathbb{T}}^{d}=\frac{\mathbb{R}^{d}}{\mathbb{Z}^{d}}$. In fact,
requiring the Hamiltonian to be radial was only needed to handle the
regularity near the boundary $\partial\Omega$. Therefore, all results stated
remain true in the torus removing the assumption (2.11). In the same spirit,
the results could be extended to a compact Riemannian manifold without
boundary, although this requires to use Bochner’s formula to handle the
gradient estimates in the Riemannian setting. We will exploit this case in
future work.

##  7 Appendix: existence of solutions to the elliptic problem

Here we show the existence of solutions to problem (3.4) where, for
simplicity, we fix $\varepsilon=1$, $T=1$. In order to construct a solution,
we follow the classical approach and we build a family of one-parameter
problems. For $\tau\in[0,1]$, we set

$m^{\tau}:=\varphi^{\tau}(-u_{t}+H(Du)-\tau
V(x))\,,\qquad\varphi^{\tau}:=(\tau\,f(\cdot)+\log(\cdot))^{-1}$ (7.1)

and we consider the problem

$\begin{cases}-{\rm
tr}({\mathcal{A}}^{\tau}(x,{\mathcal{D}}u){\mathcal{D}}^{2}u)+\tau DV(x)\cdot
H_{p}(Du)+\rho u=0&\hbox{in $(0,1)\times\Omega$,}\\\
-u_{t}+H(Du)-\tau(f(m_{1})+\tilde{V}(x))-\delta u=\psi_{1}(x)&\hbox{at $t=1$,
$x\in\Omega$,}\\\ -u_{t}+H(Du)-\tau(f(m_{0})+\tilde{V}(x))+\delta
u=\psi_{0}(x)&\hbox{at $t=0$, $x\in\Omega$,}\\\
H_{p}(Du)\cdot\vec{\nu}=0&\text{on $(0,1)\times\partial\Omega$,}\end{cases}$
(7.2)

where

${\mathcal{A}}^{\tau}=\left(a_{ij}^{\tau}((t,x),{\mathcal{D}}u)\right):=\begin{pmatrix}1&-H_{p}(Du)\\\
-H_{p}(Du)&H_{p}(Du)\otimes
H_{p}(Du)\end{pmatrix}+(1+\tau\,m^{\tau}\,f^{\prime}(m^{\tau}))\begin{pmatrix}0&0\\\
0&H_{pp}\end{pmatrix}\,$

and where $\tilde{V},\psi_{0},\psi_{1}\in C^{1,\alpha}(\overline{\Omega})$.
Problem (3.4) actually corresponds to $\tau=1$, $\tilde{V}=V$,
$\psi_{0}=\log(m_{0}),\psi_{1}=\log(m_{1})$.

To this problem we are going to apply Theorem 1.2 (Chapter X, page 462) in
[19]. Unfortunately, we can not apply directly this result because our domain
does not look regular in the time-space environment; otherwise said, our
nonlinear Neumann boundary condition is only piecewisely defined at the
$(t,x)$ boundary. To overcome this technical issue, we will use the radial
structure of the Hamiltonian $H$ (see (2.11)). This assumption simplifies the
Neumann condition on $\partial\Omega$ and allows us to follow the classical
steps up to using an extra reflection argument, in order to infer the second
order estimates. A similar result for general Hamiltonians (without requiring
(2.11)) would need to develop different technical tools to handle the $C^{1}$
estimates for piecewisely defined nonlinear boundary value problems; this is
beyond our present scopes. We recall that this kind of quasilinear problems
(with nonlinear boundary conditions) can be solved provided one shows
(uniformly in $\tau$):

  * A global Lipschitz bound for (smooth) solutions of (7.2)

  * A global $C^{1,\alpha}$ estimate for (smooth) solutions of (7.2)

  * Solvability and $C^{2,\alpha}$ estimates for the linearized version of (7.2)

We split those three tasks in the next propositions.

###### Proposition 7.1

Assume that $V\in C^{2,\alpha}(\overline{\Omega})$, that $H\in
C^{3,1}(\mathbb{R}^{d})$ satisfies (2.11) and (2.5)-(2.6) for some $q>1$ and
$\varpi>0$, and that $f\in C^{2,1}(0,\infty)$ satisfies

$\exists\,\,\beta,s_{0}>0\,:\,|f^{\prime\prime}(s)|\,s^{2}\leq\beta\,(1+f^{\prime}(s)s)\qquad\hbox{
for all $s\geq s_{0}$. }$ (7.3)

Let $m_{0},m_{1},\tilde{V},\psi_{0},\psi_{1}\in
C^{1,\alpha}(\overline{\Omega})$. Then there exists $M$, independent of
$\tau\in[0,1]$, such that any $u\in C^{2,\alpha}(\overline{Q})$ which is a
solution of (7.2) satisfies

$\|u\|_{\infty}+\|{\mathcal{D}}u\|_{\infty}\leq M$ (7.4)

where $M$ depends on
$\delta,\|m_{0}\|_{W^{1,\infty}(\Omega)},\|m_{1}\|_{W^{1,\infty}(\Omega)},\|\psi_{0}\|_{W^{1,\infty}(\Omega)},\|\psi_{1}\|_{W^{1,\infty}(\Omega)},\|\tilde{V}\|_{W^{1,\infty}(\Omega)},\|V\|_{W^{2,\infty}(\Omega)}$
and on $H,f$.

  * _Proof_.

This is the result of Theorem 3.4 (together with Lemma 3.2), up to minor
remarks. The vector field $V$ is here replaced by $\tau V(x)$ in the interior
and by $\tau\tilde{V}$ on the boundary, but this makes no difference in the
proof. The condition (7.3) is slightly more restrictive than (3.8) in Theorem
3.4; this way we can assert that $\tau f(\cdot)$ satisfies the same condition
as $f$, with the same constant. Hence, the gradient estimates are independent
of $\tau$. We point out that the regularity required on $V,H,f$ guarantee that
any $C^{2}$ solution $u$ actually belongs to $C^{3}(Q)$, so that the gradient
bounds can be classically derived. Finally, we notice that the bound on $u$
depends on $\delta$ from maximum principle, as in Lemma 3.2.  $\square$

###### Proposition 7.2

Under the same assumptions of Proposition 7.1, there exists $C$, independent
of $\tau\in[0,1]$, such that any $u\in C^{2,\alpha}(\overline{Q})$ which is a
solution of (7.2) satisfies

$\|u\|_{C^{1,\alpha}(\overline{Q})}\leq C$ (7.5)

where $C$ depends on the constant $M$ in (7.4) (as well as on the bounds of
$H_{p},H_{pp},H_{ppp}$ for $|p|\leq M$) and again on
$\|m_{0}\|_{W^{1,\infty}(\Omega)},\|m_{1}\|_{W^{1,\infty}(\Omega)},\|\psi_{0}\|_{W^{1,\infty}(\Omega)},\|\psi_{1}\|_{W^{1,\infty}(\Omega)},\|\tilde{V}\|_{W^{1,\infty}(\Omega)}$.

  * _Proof_.

The elliptic problem in (7.2) can be rephrased as

$\begin{cases}L_{\tau}(u):=-a_{ij}^{\tau}(\eta,{\mathcal{D}}u)u_{\eta_{i}\eta_{j}}+\tau
DV(x)\cdot H_{p}(Du)+\rho u=0&\hbox{ in $Q$,}\\\
N^{\tau}(u)=(\psi_{1},\psi_{0})&\hbox{on $\Sigma$,}\\\
H_{p}(Du)\cdot\vec{\nu}=0&\text{on $(0,1)\times\partial\Omega$,}\end{cases}$
(7.6)

where $\eta=(t,x)$, $\Sigma=(\\{1\\}\times\Omega)\cup(\\{0\\}\times\Omega)$
and the boundary operator $N^{\tau}$ is defined as

$N_{\tau}(u):=\begin{cases}-u_{t}+H(Du)-\tau(f(m_{1})+\tilde{V}(x))-\delta
u&\hbox{at $t=1$, $x\in\Omega$,}\\\
-u_{t}+H(Du)-\tau(f(m_{0})+\tilde{V}(x))+\delta u&\hbox{at $t=0$,
$x\in\Omega$.}\end{cases}$ (7.7)

Due to the estimate (7.4) in Proposition 7.1, we have that the coefficients
$a_{ij}^{\tau}$ need only be considered in the compact set
$\mathcal{M}:=\\{\eta\in\overline{Q},|u|\leq M,|{\mathcal{D}}u|\leq M\\}$.
Assuming $V$ of class $C^{2,\alpha}$, $f$ of class $C^{2,1}$ and $H$ of class
$C^{3,1}$, we get that the coefficients $a_{ij}^{\tau}$, as well as the first
order coefficients, are $C^{1,\alpha}$ in their arguments, with
$\frac{\partial a_{ij}^{\tau}(y,p)}{\partial p}$ being Lipschitz in $p$. In
particular, for the interior regularity we can apply the classical results
(see e.g. [19, Chapter 6, Thm 1.1]) and we deduce the $C^{1,\alpha}$ interior
estimate for $u$.

We are only concerned now with the boundary regularity. Without loss of
generality (because $H_{pp}$ is bounded below), we assume here that

$H(p)=h(|p|^{2})\,,$ (7.8)

for some $C^{2}$ function $h$ such that $h^{\prime}>0$. With this notation we
have $H_{p}=2h^{\prime}\,p$ and
$H_{pp}=2h^{\prime}I_{d}+4h^{\prime\prime}(p\otimes p)$. Now we use a (local)
reflection argument in the $x$ variable. To this purpose, recall that for
$x\in\Omega$, we denote by $d(x)$ the distance of $x$ to the boundary. It is
well known (see e.g. [15]) that if $\Omega$ is of class $C^{k}$, $d(x)$ is a
$C^{k}$ function in a neighborhood of the boundary. More precisely, there
exists a positive number $d_{0}>0$ such that, if we denote
$\Omega_{d_{0}}:=\\{x\in\Omega\,:\,d(x)<d_{0}\\}$, then any point
$x\in\Omega_{d_{0}}$ can be represented as

$x=-\delta\,\vec{n}(\bar{x})+\bar{x}$ (7.9)

where $\delta=d(x)$ and $\bar{x}$ is the projection of $x$ onto
$\partial\Omega$ ($\vec{n}(\bar{x})$ is the outward unit vector at $\bar{x}$).
To ease notations, we denote hereafter by $N$ the dimension of the euclidean
space (i.e. $\Omega\subset\mathbb{R}^{N}$). Therefore, if we represent
$\bar{x}\in\partial\Omega$ in a system of coordinates of some local chart, we
may assume that $\bar{x}=(y^{\prime},g(y^{\prime}))$ for some $C^{k}$ function
$g$ defined on some open subset $U\subset\mathbb{R}^{N-1}$. Hence (7.9) yields
a local diffeomorphism between $\Omega_{d_{0}}$ and $U\times(0,d_{0})$ defined
as

$T(x):=(y^{\prime},\delta)\quad\iff\quad
x=-\delta\,\vec{n}(\bar{x})+\bar{x}\,,\quad\bar{x}=(y^{\prime},g(y^{\prime}))\in\partial\Omega.$

Let us now set $y_{N}:=\delta$ so that $y=(y^{\prime},y_{N})$ and $y=T(x)$
will shortly denote the above change of variables. We notice that assuming
$\Omega$ of class $C^{3}$ guarantees that the mapping $T$ is of class $C^{2}$.
The advantage of this change of coordinates lies in the property that

$\nabla T_{k}(x)\cdot\nabla T_{N}(x)=\nabla T_{k}(x)\cdot\nabla
d(x)=0\quad\forall k=1,\ldots,N-1$ (7.10)

because $T_{N}(x)=d(x)$ (the distance function) and the projection coordinates
remain constant if we move along the normal direction (i.e. $\nabla
T_{k}(x)\cdot\nabla d(x)=0$ for all $N-1$-coordinates $y_{k}=T_{k}(x)$).

Now we set

$\tilde{u}(t,y):=u(t,T^{-1}y)\,\qquad\rightsquigarrow\qquad
u(t,x)=\tilde{u}(t,T(x))\,$

and we compute:

$u_{x_{i}}=\sum_{k=1}^{N}\tilde{u}_{y_{k}}\frac{\partial T_{k}}{\partial
x_{i}}\quad;\quad
u_{x_{i}x_{j}}=\sum_{k,\ell=1}^{N}\tilde{u}_{y_{k}y_{\ell}}\frac{\partial
T_{k}}{\partial x_{i}}\frac{\partial T_{\ell}}{\partial
x_{j}}+\sum_{k=1}^{N}\tilde{u}_{y_{k}}\frac{\partial^{2}T_{k}}{\partial
x_{i}\partial x_{j}}$ (7.11)

so the equation (7.6) is transformed into

$\tilde{u}_{tt}-2\sum_{k}\left(DT_{k}\cdot
H_{p}(Du)\right)\tilde{u}_{ty_{k}}+\sum_{k,\ell=1}^{N}\tilde{a}^{\tau}_{kl}\tilde{u}_{y_{k}y_{\ell}}+\sum_{k=1}^{N}\tilde{b}_{k}^{\tau}\,\tilde{u}_{y_{k}}=\rho\tilde{u}$
(7.12)

where

$\tilde{a}^{\tau}_{kl}:=\sum_{i,j}a_{ij}^{\tau}\,\frac{\partial
T_{k}}{\partial x_{i}}\frac{\partial T_{\ell}}{\partial
x_{j}}\,;\qquad\tilde{b}_{k}^{\tau}:=\sum_{i,j}a_{ij}^{\tau}\,\frac{\partial^{2}T_{k}}{\partial
x_{i}\partial x_{j}}-\tau 2h^{\prime}DV\cdot DT_{k}$ (7.13)

and where the coefficients $a_{ij}^{\tau}$ are computed on $Du$; recall that
$Du=MD\tilde{u}$ for $M=(m_{ij})=\left(\frac{\partial T_{j}}{\partial
x_{i}}\right)$.

We extend now $\tilde{u}$ for $y_{N}<0$ by setting

$u^{*}(t,y):=\begin{cases}\tilde{u}(t,y^{\prime},y_{N})&\hbox{if
$y^{\prime}\in U,y_{N}\in(0,d_{0})$}\\\
\tilde{u}(t,y^{\prime},-y_{N})&\hbox{if $y^{\prime}\in
U,y_{N}\in(-d_{0},0)$}\end{cases}$

Notice that (7.10) implies

$Du(t,x)\cdot Dd(x)=\sum_{k=1}^{N}\tilde{u}_{y_{k}}(DT_{k}\cdot
Dd(x))=\tilde{u}_{y_{N}}\,|Dd(x)|^{2}=\tilde{u}_{y_{N}}\,\,,\qquad
x\in\Omega_{d_{0}}$

so the Neumann condition $H_{p}(Du)\cdot\vec{\nu}=0$ together with (7.8)
translate into $\tilde{u}_{y_{N}}=0$; this implies that $u^{*}$ is a $C^{1}$
function. Now we look at the equation satisfied by $u^{*}$; this involves the
even reflection (across $y_{N}=0$) of coefficients $\tilde{a}^{\tau}_{kl}$ for
$k,\ell\neq N$ or $k=\ell=N$ and the odd reflection of coefficients
$\tilde{a}^{\tau}_{kl}$ for $k\neq\ell$ and $k$ or $\ell=N$. Let us check that
the structure condition (7.8) and the property (7.10) yield a nice equation
for $u^{*}$. First of all we observe that, due to (7.10),

$\displaystyle|Du|^{2}$
$\displaystyle=\sum_{i}u_{x_{i}}^{2}=\sum_{i}\left(\tilde{u}_{y_{N}}\frac{\partial
T_{N}}{\partial x_{i}}+\sum_{p=1}^{N-1}\tilde{u}_{y_{p}}\frac{\partial
T_{p}}{\partial x_{i}}\right)^{2}$
$\displaystyle=\tilde{u}_{y_{N}}^{2}+\sum_{i}\left(\sum_{p=1}^{N-1}\tilde{u}_{y_{p}}\frac{\partial
T_{p}}{\partial
x_{i}}\right)^{2}=|u^{*}_{y_{N}}|^{2}+|M^{\prime}D^{\prime}u^{*}|^{2}$
$\displaystyle\qquad\hbox{for $M^{\prime}_{ip}=(\frac{\partial T_{p}}{\partial
x_{i}})\,,\,\,D^{\prime}u^{*}=(u^{*}_{y_{p}})_{p=1,\ldots,N-1}$}\quad\rightsquigarrow\quad(M^{\prime}D^{\prime}u^{*})_{i}:=\sum_{p=1}^{N-1}u^{*}_{y_{p}}\frac{\partial
T_{p}}{\partial x_{i}}$

In particular, by definition of $m^{\tau}$ in (7.1), we have

$(1+\tau\,m^{\tau}\,f^{\prime}(m^{\tau}))={\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})$
(7.14)

for some smooth function ${\mathcal{\kappa}}^{\tau}$, where $V^{*}(y)$ is the
even reflection of $V(T^{-1}(y))$.

Then we compute from (7.13), using the definition of $a_{ij}^{\tau}$ and the
radial structure of $H$, and due to condition (7.10),

$\displaystyle\forall k\neq N,\qquad\tilde{a}^{\tau}_{kN}$
$\displaystyle=\sum_{i,j}\left[H_{p_{i}}\otimes
H_{p_{j}}+H_{p_{i}p_{j}}(1+\tau\,m^{\tau}\,f^{\prime}(m^{\tau}))\right]\frac{\partial
T_{k}}{\partial x_{i}}\frac{\partial T_{N}}{\partial x_{j}}$
$\displaystyle=\sum_{i,j}\left(4(h^{\prime})^{2}+4h^{\prime\prime}{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\right)u_{x_{i}}u_{x_{j}}\frac{\partial
T_{k}}{\partial x_{i}}\frac{\partial T_{N}}{\partial x_{j}}$
$\displaystyle=\left(4(h^{\prime})^{2}+4h^{\prime\prime}{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\right)\sum_{p=1}^{N-1}DT_{p}\cdot
DT_{k}\tilde{u}_{y_{p}}\tilde{u}_{y_{N}}$

where we used (7.11) (and once more (7.10)) in the last step. In particular,
we notice that $\tilde{a}^{\tau}_{kN}$ is an odd coefficient with respect to
reflection; hence, using the above notation,

$\forall k\neq
N,\qquad\tilde{a}^{\tau}_{kN}\tilde{u}_{y_{k}y_{N}}=\left(4(h^{\prime})^{2}+4h^{\prime\prime}{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\right)(DT_{k}\cdot
MD^{\prime}u^{*})u^{*}_{y_{N}}u^{*}_{y_{k}y_{N}}\,.$

With similar computations we find

$\displaystyle\tilde{a}^{\tau}_{NN}$
$\displaystyle=\left(4(h^{\prime})^{2}+4h^{\prime\prime}{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\right)(u^{*}_{y_{N}})^{2}+2h^{\prime}\,{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\,.$

Still using (7.8), (7.11) and (7.10), we observe that

$\displaystyle\sum_{k}\left(DT_{k}\cdot H_{p}(Du)\right)\tilde{u}_{ty_{k}}$
$\displaystyle=2h^{\prime}(\tilde{u}_{y_{N}}\tilde{u}_{ty_{N}}+\sum_{k,\ell=1}^{N-1}\left(DT_{k}\cdot
DT_{\ell})\tilde{u}_{y_{\ell}}\tilde{u}_{ty_{k}}\right)$
$\displaystyle=2h^{\prime}(u^{*}_{y_{N}}u^{*}_{ty_{N}}+\sum_{k,\ell=1}^{N-1}\left(DT_{k}\cdot
DT_{\ell}\right)u^{*}_{y_{\ell}}u^{*}_{ty_{k}})$

Therefore, $u^{*}(t,y)$ solves the following elliptic equation for
$t\in(0,1),y\in B$, where $B\subset\mathbb{R}^{N}$ is any open set which is
contained in $U\times(-d_{0},d_{0})$:

$u_{tt}^{*}-4h^{\prime}(u^{*}_{y_{N}}u^{*}_{ty_{N}}+\sum_{k,\ell=1}^{N-1}\left(DT_{k}\cdot
DT_{\ell}\right)u^{*}_{y_{\ell}}u^{*}_{ty_{k}})+\sum_{k,\ell=1}^{N}a^{*}_{kl}u^{*}_{y_{k}y_{\ell}}+b^{*}(t,y)=\rho
u^{*}$ (7.15)

where

$\displaystyle{\small
a^{*}_{kl}=\left(4(h^{\prime})^{2}+4h^{\prime\prime}{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\right)(DT_{k}\cdot
M^{\prime}D^{\prime}u^{*})(DT_{\ell}\cdot
M^{\prime}D^{\prime}u^{*})+2h^{\prime}\,{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})(DT_{k}\cdot
DT_{\ell})}\quad\hbox{if $k,\ell\neq N$,}$ $\displaystyle
a^{*}_{kN}=\left(4(h^{\prime})^{2}+4h^{\prime\prime}{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\right)\left(DT_{k}\cdot
M^{\prime}D^{\prime}u^{*}\right)u^{*}_{y_{N}}\quad\hbox{if $k\neq N$,}$
$\displaystyle
a^{*}_{NN}=\left(4(h^{\prime})^{2}+4h^{\prime\prime}{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\right)|u^{*}_{y_{N}}|^{2}+2h^{\prime}\,{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})$

and where $b^{*}(t,y)$ is the even extension (through $y_{N}=0$) of the term
$\tilde{b}^{\tau}\cdot D\tilde{u}$. Now we apply to equation (7.15) the
$C^{1,\alpha}$-estimates (up to the boundary) for elliptic equations with
nonlinear first order condition at the boundary, see e.g. [24, Lemma 2.3]222We
warn the reader that, literally, the result in [24, Lemma 2.3] assumes the
coefficients $a_{ij}(y,{\mathcal{D}}u)$ to be $C^{1}$ with respect to $y$;
however this is not necessary, as can be readily checked by inspecting the
proof as well as by noticing that the estimate only requires $\frac{\partial
a_{ij}}{\partial y}$ to be bounded. In fact, the Lipschitz character of
$a_{ij}$ with respect to $y$ is enough, which is here preserved by the
reflection argument.. Let us point out that the boundary here is represented
by the hypersurfaces $(\Sigma_{0}\cap B)\cup(\Sigma_{1}\cap B)$; the solution
$u^{*}$ satisfies a uniform bound $|u^{*}|+|{\mathcal{D}}u^{*}|\leq K$ (by
Proposition 7.1) and the coefficients $a_{ij}^{*}((t,y),{\mathcal{D}}u^{*})$
are $C^{1}$ with respect to ${\mathcal{D}}u^{*}$ and $\frac{\partial
a_{ij}^{*}}{\partial q}$ as well as $\frac{\partial a_{ij}^{*}}{\partial y}$
are uniformly bounded (because $H$ is $C^{3,1}$ in $p$ and $V$ is Lipschitz).
Note that $a_{ij}^{*}$ depends on $y$ only through (the even reflections of)
$V(y)$ and $T^{-1}(y)$. As for the coefficient $b^{*}(t,y)$, it is actually
continuous (thanks to the fact that $u^{*}_{y_{N}}=0$) and, what only matters,
is uniformly bounded by some constant depending on $K,\|DV\|_{\infty}$ and
$\|D^{2}T\|_{\infty}$. Finally, the boundary operator $N^{\tau}$ is Lipschitz
with respect to $y$ and $C^{2}$ with respect to the gradient variable.
According to [24, Lemma 2.3], ${\mathcal{D}}u^{*}$ satisfies a bound in the
Hölder norm, which is uniform with respect to $\tau$. With a standard
localization argument, based on a partition of unity, we conclude the bound
(7.5).  $\square$

###### Remark 7.3

We stress that a $C^{1,\alpha}$ estimate (up to the boundary) can also be
obtained if we look at the elliptic equation in the divergence form, which is
clearly inherited from the continuity equation in (2.1). In the notations used
in the previous proof, we observe that the function $u^{*}$, obtained after
local change of coordinates and reflection, satisfies the equation

$\displaystyle\partial_{t}m^{\tau}-\sum_{k=1}^{N}\frac{\partial}{\partial
y_{k}}\left(m^{\tau}2h^{\prime}(|MDu^{*}|^{2})DT_{k}\cdot(MDu^{*})\right)+\rho\,u^{*}\frac{m^{\tau}}{1+\tau
f^{\prime}(m^{\tau})m^{\tau}}=0$ where
$m^{\tau}=\varphi^{\tau}\left(-u_{t}+h(|MDu^{*}|^{2})-\tau V\right)$,

where we recall that $\varphi^{\tau}$ is defined in (7.1) and
$M=(m_{ij})=(\frac{\partial T_{j}}{\partial x_{i}})$,
$|MDu^{*}|^{2}=|u_{y_{N}}^{*}|^{2}+\sum_{i=1}^{N}\left(\sum_{\ell=1}^{N-1}\frac{\partial
T_{\ell}}{\partial x_{i}}u^{*}_{y_{\ell}}\right)^{2}$. Hence the equation can
be rewritten as a nonlinear divergence form equation on $u^{*}$:

$\partial_{t}\left(A_{0}(y,{\mathcal{D}}u^{*})\right)+\sum_{k=1}^{N}\frac{\partial}{\partial
y_{k}}\left(A_{k}(y,{\mathcal{D}}u^{*})\right)=\beta(y,u^{*},{\mathcal{D}}u^{*})\qquad(t,y)\in(0,1)\times
B$

complemented with a co-normal boundary condition at $t=0,t=1$:

$\displaystyle{\bf A}\cdot{\bf n}$
$\displaystyle=\left(A_{0}(y,{\mathcal{D}}u^{*})\mathop{|}_{t=1},-A_{0}(y,{\mathcal{D}}u^{*})\mathop{|}_{t=0}\right)$
$\displaystyle=\left(\varphi^{\tau}(\psi_{1}^{*}+\tau(f(m_{1})^{*}+\tilde{V}^{*}-V^{*})+\delta
u^{*}),-\varphi^{\tau}(\psi_{0}^{*}+\tau(f(m_{0})^{*}+\tilde{V}^{*}-V^{*})-\delta
u^{*})\right)$

where $(\cdot)^{*}$ denotes the even extension of the various functions.
Notice that when $\tilde{V}=V$, $\delta=0$ and $\psi_{i}=\log(m_{i})$, then we
have ${\bf A}\cdot{\bf n}=(m_{1}^{*},-m_{0}^{*})$.

Since $|\beta(y,u^{*},{\mathcal{D}}u^{*})|\leq C$ (because
${\mathcal{D}}u^{*}$ is bounded), and due to the regularity of the boundary
terms, if $H_{pp}$ is nondegenerate the $C^{1,\alpha}$ estimate for $u^{*}$
can then be deduced by [25, Thm 2].

We mention this alternative approach because it could be exploited if one aims
at generalizing the $C^{1,\alpha}$ estimate to possibly degenerate
Hamiltonians, e.g. satisfying (2.5) with $\varpi=0$.

Last ingredient is the well-posedness of the linearized problem around one
solution $u$ of (7.2) and, correspondingly, a compactness property for
sequences of solutions of (7.2). In order to use the $C^{2,\alpha}$
regularity, here we need to use a compatibility condition at the boundary
$\partial\Omega$. To this purpose, we denote
$C^{1,\alpha}_{N}(\overline{\Omega})$ the space of functions $\phi\in
C^{1,\alpha}(\overline{\Omega})$ such that $D\phi\cdot\vec{\nu}=0$ on
$\partial\Omega$. We also denote $a_{ij}^{\tau}=a_{ij}^{\tau}(\eta,q)$.

###### Proposition 7.4

In addition to the conditions of Proposition 7.1, assume that $\partial\Omega$
is of class $C^{4}$, and that the functions $f(m_{0})$, $f(m_{1})$,
$\tilde{V}(x)$, $\psi_{0},\psi_{1}$ belong to
$C^{1,\alpha}_{N}(\overline{\Omega})$. Then the solution $u$ of (7.6) belongs
to $C^{2,\alpha}(\overline{Q})$, and for every $\vartheta\in
C^{0,\alpha}(\overline{Q}),\zeta=(\zeta_{1},\zeta_{0})\in
C^{1,\alpha}_{N}(\overline{\Omega})$ there exists a unique $\phi\in
C^{2,\alpha}(\overline{Q})$ which is a solution of the linear problem

$\begin{cases}-a_{ij}^{\tau}((t,x),{\mathcal{D}}u)\phi_{ij}-\frac{\partial
a^{\tau}_{ij}}{\partial q_{k}}((t,x),{\mathcal{D}}u)u_{ij}\phi_{k}+\tau
H_{pp}(Du)DV\cdot D\phi+\rho\phi=\vartheta&\hbox{in $Q$ }\\\
-\phi_{t}+H_{p}(Du)\cdot D\phi=\delta\phi+\zeta_{1}&\hbox{on $\Sigma_{T}$,}\\\
-\phi_{t}+H_{p}(Du)\cdot D\phi+\delta\phi=\zeta_{0}&\hbox{on $\Sigma_{0}$,}\\\
D\phi\cdot\vec{\nu}=0&\text{on $(0,T)\times\partial\Omega$.}\end{cases}$
(7.16)

Moreover, if $\psi_{1m},\psi_{0m}$ are sequences which converge in
$C^{1,\alpha}_{N}(\overline{\Omega})$, then the corresponding solutions
$u_{m}$ of (7.6) are relatively compact in $C^{2,\alpha}(\overline{Q})$.

  * _Proof_.

The first assertion follows from Proposition 7.2 if we come back to the
equation (7.15) complemented with the boundary condition satisfied by $u^{*}$,
which reads as

$\begin{cases}-u_{t}^{*}+h(|MDu^{*}|^{2})-\delta
u^{*}=\tau(f(m_{1})^{*}+\tilde{V}^{*})+\psi_{1}^{*}&\quad\hbox{at $t=1$,}\\\
-u_{t}^{*}+h(|MDu^{*}|^{2})+\delta
u^{*}=\tau(f(m_{0})^{*}+\tilde{V}^{*})+\psi_{0}^{*}&\quad\hbox{at
$t=0$.}\end{cases}$

Since the boundary data are assumed to belong to
$C^{1,\alpha}_{N}(\overline{\Omega})$, they are reflected into $C^{1,\alpha}$
functions. This implies that $u^{*}_{t}-h(|MDu^{*}|^{2})\in C^{1,\alpha}(U)$.
One can therefore apply the results of [23] to invoke that $u^{*}\in
C^{2,\alpha}$, which leads to the $C^{2,\alpha}$ regularity of $u$ up to the
boundary.

Let us now consider the linearized problem (7.16). Since $u\in
C^{2,\alpha}(\overline{Q})$, (7.16) is a linear problem with Hölder
coefficients, and Schauder’s theory applies to get, at least, the interior
estimates of $\phi$. Let us only check the boundary regularity. With the
notations of Proposition 7.2, we localize near the boundary with the change of
variable $y=T(x)$, and then we extend the problem by reflection, setting
$\phi^{*}(t,y)$ the (even) reflection of
$\tilde{\phi}(t,y):=\phi(t,T^{-1}(y))$ through the hyperplane $\\{y_{N}=0\\}$.
From the computations of Proposition 7.2, we know that the matrix
${\mathcal{A}}^{\tau}$ is transformed into a new matrix ${\mathcal{A}^{*}}$
which depends on $(y,{\mathcal{D}}u^{*})$ as follows:

${\mathcal{A}^{*}}=\left(a_{ij}^{*}\right)=\begin{pmatrix}1&-2h^{\prime}\Lambda_{u}\\\
-2h^{\prime}\,\Lambda_{u}&[4(h^{\prime})^{2}]\Lambda_{u}\otimes\Lambda_{u}\end{pmatrix}+{\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})\begin{pmatrix}0&0\\\
0&(2h^{\prime})\left(DT_{i}\cdot
DT_{j}\right)+(4h^{\prime\prime})\Lambda_{u}\otimes\Lambda_{u}\end{pmatrix}$

with $\Lambda_{u}:=\left((DT_{k}\cdot M^{\prime}D^{\prime}u^{*})_{k=1,\dots
N-1},u^{*}_{y_{N}}\right)$ and
${\mathcal{\kappa}}^{\tau}(V^{*}(y),{\mathcal{D}}u^{*})$ defined as in (7.14)
.

We look now at the first order terms. By definition of $a_{ij}^{\tau}$, we
have

$\begin{split}\frac{\partial a^{\tau}_{ij}}{\partial
q_{k}}((t,x),{\mathcal{D}}u)u_{ij}\phi_{k}&=-2H_{pp}Du_{t}\cdot
D\phi+2D^{2}uH_{p}\cdot
H_{pp}D\phi+\tau\,g^{\prime}(m^{\tau})(\varphi^{\tau})^{\prime}\,[-\phi_{t}+H_{p}D\phi]\,{\rm
tr}(H_{pp}D^{2}u)\\\ &\qquad\quad+(1+\tau\,g(m^{\tau})){\rm tr}(H_{ppp}D\phi
D^{2}u)\,\end{split}$ (7.17)

where $g(s)=sf^{\prime}(s)$. Using the radial structure of $H$ ((7.8)) and the
chosen reference frame where (7.10) holds true, one can check (with long but
routine computations, similarly as in Proposition 7.2) that all terms in
(7.17) can be transformed into functions of $\phi^{*}$ and $u^{*}$ except for
a few terms which are only first order in $u$ and terms which involve second
derivatives of the map $T$. Then, after the change of coordinates and the
reflection, $\phi^{*}$ satisfies a linear elliptic problem of the following
form in the open set $(0,T)\times B$:

$a_{ij}^{*}\phi^{*}_{ij}+B^{*}_{i}\phi^{*}_{i}+c^{*}+\rho\phi^{*}=\vartheta^{*}\qquad\hbox{$(t,y)\in(0,T)\times
B$}\,.$ (7.18)

Now we wish to apply to $\phi^{*}$ the Hölder estimates up to the boundary in
order to conclude that $\phi^{*}\in C^{2,\alpha}([0,T]\times
U\times(-d_{0},d_{0}))$ for some open set $U\subset\mathbb{R}^{N-1}$. The
regularity of $\phi^{*}$ will yield the $C^{2,\alpha}$ regularity of $\phi$ in
$[0,T]\times({\mathcal{B}}\cap\overline{\Omega})$ for any ball $\mathcal{B}$
such that $\partial\Omega\cap{\mathcal{B}}$ is the portion of a smooth graph.

In (7.18), the ($d+1$-dimensional) vector field $B^{*}$ depends on
${\mathcal{D}}^{2}u^{*}$ (which is Hölder continuous) and on the (even
reflection) of $DT,DV,D\tilde{u}$, which are Lipschitz continuous. Therefore,
overall $B^{*}$ belongs to $C^{0,\alpha}$ for some $\alpha\in(0,1)$.
Conversely, the term $c^{*}$ in (7.18) is the even reflection of a function
depending on $D\phi$, $D\tilde{u}$ and $D^{2}T$; in this case we cannot say
that $c^{*}$ is reflected into a Hölder continuous function until we establish
the Hölder regularity for $D\phi$.

Therefore, since the coefficient $c^{*}$ does not belong to $C^{0,\alpha}$ a
priori, we cannot apply the Schauder’s estimates in one shot to infer the
$C^{2,\alpha}$ regularity of $\phi^{*}$. However, the coefficients
$a_{ij}^{*}$ belong to $C^{0,1}$ (because $u^{*}$ is $C^{2,\alpha}$, while
$V,DT$ are $C^{1}$ and so they are reflected into a Lipschitz function);
hence, by rewriting the equation in divergence form, $\phi^{*}$ is, in
particular, a bounded weak solution to the equation

$(a_{ij}^{*}\phi^{*}_{j})_{i}+H(y,\phi^{*},\phi^{*}_{i})=0\qquad(t,y)\in(0,T)\times
U\times(-d_{0},d_{0})$ (7.19)

for some function $H(y,\phi^{*},\phi^{*}_{i})$ which satisfies
$|H(y,\phi^{*},\phi^{*}_{i})|\leq C_{0}(1+|D\phi^{*}|)$. In addition,
$\phi^{*}$ satisfies a co-normal boundary condition at
$\Sigma_{0},\Sigma_{T}$. By regularity for divergence form equations (which is
true for general nonlinear problems, see e.g. [25, Thm 2]), we deduce that
$\phi^{*}\in C^{1,\alpha}$ and

$\|\phi^{*}\|_{C^{1,\alpha}}\leq K\qquad\hbox{for some
$K=K(\|\phi\|_{\infty},\|u\|_{C^{2,\alpha}},\|T\|_{2},\|DV\|_{\infty},\|\vartheta\|_{\infty},\|\zeta\|_{\alpha},Q)$.}$

We can now use this information, which implies that the term $c^{*}$ in (7.18)
belongs to $C^{0,\alpha}$. Since $a^{*}_{ij}$, $B^{*}$, $\vartheta^{*}$ also
belong to $C^{0,\alpha}$, and since the boundary data $\zeta^{*}\in
C^{1,\alpha}$ (because $\zeta$ satisfies a Neumann condition), the classical
Schauder’s estimates (see e.g. [23, Lemma 1]) imply that $\phi^{*}\in
C^{2,\alpha}$ and

$\|\phi^{*}\|_{C^{2,\alpha}}\leq K\qquad\hbox{for some
$K=K(\|\phi\|_{\infty},\|u\|_{C^{2,\alpha}},\|T\|_{C^{2,\alpha}},\|V\|_{C^{1,\alpha}},\|\vartheta\|_{\alpha},\|\zeta\|_{C^{1,\alpha}},Q)$.}$

We notice that $\|\phi\|_{\infty}$ is estimated uniformly by maximum
principle. We also recall that, requiring $\partial\Omega$ of class $C^{4}$,
we have that $T$ is of class $C^{3}$. Finally, we proved that all solutions to
the linear problem (7.16) satisfy the estimate

$\|\phi\|_{C^{2,\alpha}}\leq
C\left(\|\vartheta\|_{\alpha}+\|\zeta\|_{C^{1,\alpha}}\right)$

for some $C=C(\|u\|_{C^{2,\alpha}},\partial\Omega,\|V\|_{C^{1,\alpha}})$. By
linear theory, this implies that problem (7.16) is uniquely solvable for every
$\tau\in[0,1]$.

Now we prove the last assertion of the Proposition. To this goal, we first
observe that, since $\psi_{m}=(\psi_{0m},\psi_{1m})$ is convergent in
$C^{1,\alpha}_{N}(\overline{\Omega})$, by Proposition 7.2 the sequence $u_{m}$
is bounded in $C^{1,\alpha}(\overline{Q})$, and therefore it is relatively
compact in $C^{1}(\overline{Q})$ by Ascoli-Arzelà theorem. We observe that
$w^{m,n}:=u_{m}-u_{n}$ solves

$\left\\{\begin{aligned}
&-a_{ij}^{\tau}(x,{\mathcal{D}}u_{m})w^{m,n}_{ij}+B_{k}^{n,m}w^{m,n}_{k}+\tau\,\Gamma^{n,m}DV(x)\cdot
Dw^{m,n}+\rho w^{m,n}=0&&\text{in }(0,1)\times\Omega,\\\
&-\partial_{t}w^{m,n}+\beta^{n,m}\cdot Dw^{m,n}-\delta
w^{m,n}=\psi_{1n}-\psi_{1m}&&\hbox{at $t=1$, $x\in\Omega$,}\\\
&-\partial_{t}w^{m,n}+\beta^{n,m}\cdot Dw^{m,n}+\delta
w^{m,n}=\psi_{0n}-\psi_{0m}&&\hbox{at $t=0$, $x\in\Omega$,}\\\
&Dw^{m,n}\cdot\vec{\nu}=0&&\text{on
$(0,1)\times\partial\Omega$,}\end{aligned}\right.$

where

$\displaystyle
B_{k}^{n,m}=\sum_{ij}(u_{n})_{ij}\int_{0}^{1}\frac{a_{ij}^{\tau}}{\partial
q_{k}}(x,\lambda Du_{m}+(1-\lambda)Du_{n})d\lambda$
$\displaystyle\beta^{n,m}=\int_{0}^{1}H_{p}(\lambda
Du_{m}+(1-\lambda)Du_{n})d\lambda\,,\qquad\Gamma^{n,m}=\int_{0}^{1}H_{pp}(\lambda
Du_{m}+(1-\lambda)Du_{n})d\lambda\,.$

This is a linear problem for $w^{m,n}$ which is of the same kind as (7.16).
Therefore, we proceed as we did before, by localizing near the boundary and
reflecting the solution through (the straightened part of) $\partial\Omega$.
Then we apply to the reflected function $(w^{n,m})^{*}$ the $C^{2,\alpha}$
estimates as before ([23, Lemma 1]). We only stress that, differently than it
was for the function $\phi^{*}$ above, now we already know that $w^{n,m}$ is
bounded in $C^{1,\alpha}(\overline{Q})$; this allows us to apply directly the
Schauder estimates because (in the notation used before)
$a_{ij}^{*}(w^{n,m})^{*}_{ij}$ belongs to (and can be estimated in)
$C^{0,\alpha}(\overline{Q})$. Applying the precise form of [23, Lemma 1], we
conclude that

$\|w^{n,m}\|_{C^{2,\alpha}}\leq
C_{0}\|w^{n,m}\|_{C^{1}}(1+\|u_{n}\|_{C^{2,\alpha}}+\|u_{m}\|_{C^{2,\alpha}})+C_{1}\|w^{n,m}\|_{C^{1,\alpha}}(\|D^{2}u_{n}\|_{\infty}+\|D^{2}u_{m}\|_{\infty})\,.$
(7.20)

From this estimate we conclude as in [23, Lemma 2]; one shows first that
$\|u_{n}\|_{C^{2,\alpha}}$ is bounded (using the smallness of
$\|w^{n,m}\|_{C^{1}}$ in (7.20)) and this implies the compactness of $u_{n}$
in $C^{C^{1,\alpha}}$. Hence (7.20) yields the convergence of $u_{n}$ in
$C^{2,\alpha}(\overline{Q})$.  $\square$

We are ready to prove the existence of solutions to the elliptic problem
(3.4).

Proof of Theorem 3.8. We start by assuming that the data are more regular,
namely that $\partial\Omega$ is of class $C^{4}$, $V\in
C^{2,\alpha}(\overline{\Omega})$, $f\in C^{2,1}(0,\infty)$ satisfies (7.3) and
$H\in C^{3,1}(\mathbb{R}^{d})$ satisfies (7.8). We also assume that the
boundary data $f(m_{0})$, $f(m_{1})$, $\tilde{V}(x)$, $\psi_{0},\psi_{1}$
belong to $C^{1,\alpha}_{N}(\overline{\Omega})$. Under those conditions, we
can use Propositions 7.1–7.4 and we establish the existence of a solution with
a continuity method, following the same proof as in [19, Chapter X, Thm 1.1 &
Thm 1.2]. Let us sketch some detail: we set

$X:=\\{v\in
C^{2,\alpha}(\overline{Q}):Dv\cdot\vec{\nu}=0\,\,\hbox{on}\,(0,T)\times\partial\Omega\\}$

and $X^{\prime}:=C^{0,\alpha}(\overline{Q})\times
C^{1,\alpha}(\overline{\Omega})^{2}$, and we define a family of mappings from
$X$ into $X^{\prime}$ by setting
$\Phi(u,\tau):=(L^{\tau}(u),N^{\tau}_{0}(u),N^{\tau}_{1}(u))$, where
$N^{\tau}_{0},N^{\tau}_{1}$ are the restrictions to $t=0$, $t=1$ respectively,
of the boundary operator $N^{\tau}$ in (7.7). We also set $X^{\prime}_{0}$ the
closed subset of $X^{\prime}$ made of elements $(0,\psi_{0},\psi_{1})$, with
$(\psi_{0},\psi_{1})\in C^{1,\alpha}_{N}(\overline{\Omega})^{2}$. Thanks to
Proposition 7.4, the mapping $\Phi$ is locally invertible in a neighborhood of
points $u\in X$ such that $\Phi(u,\tau)\in X_{0}^{\prime}$, and in addition a
sequence $u_{k}$ is compact in $X$ if $\Phi(u_{k},\tau)$ is converging in
$X^{\prime}_{0}$. Finally, for $\tau=0$ and $\psi_{0}=\psi_{1}=0$, the problem
(7.6) has the unique solution $u=0$. Therefore, as in [19, Chapter X, Thm 1.1]
one concludes that the problem with $\tau=1$ is uniquely solvable.

In a second step, we approximate the data with smoother sequences. In
particular, we approximate the domain $\Omega$ with a sequence of $C^{4}$
domains $\Omega_{n}$, and the functions $m_{0},m_{1}$ with sequences
$m_{0n},m_{1n}\in C^{1,\alpha}(\overline{\Omega})$ such that
$\|m_{0n}\|_{W^{1,\infty}(\Omega)},\|m_{1n}\|_{W^{1,\infty}(\Omega)}$ are
bounded, as well as $\|\log(m_{0n})\|_{W^{1,\infty}(\Omega)}$ and
$\|\log(m_{1n})\|_{W^{1,\infty}(\Omega)}$. It is also possible to build
$m_{0n},m_{1n}$ in a way that those functions be constant along the normal in
a small neighborhood of $\partial\Omega$, so that they belong to
$C^{1,\alpha}_{N}(\overline{\Omega})$. Similarly we take a sequence
$\tilde{V}_{n}\in C^{1,\alpha}_{N}(\overline{\Omega})$ which converges
uniformly to $V$ and such that $\|\tilde{V}_{n}\|_{W^{1,\infty}(\Omega)}$ is
bounded. Finally, we take another sequence $V_{n}\in
C^{2,\alpha}(\overline{\Omega})$ converging to $V$ and such that
$\|V_{n}\|_{W^{2,\infty}(\Omega)}$ is bounded. As for the function $f$, we
approximate it with $f_{n}(m)=f(\frac{m}{1+\frac{1}{n}m})$ which satisfies the
stronger condition (7.3) (for some constant $\beta$ possibly dependent on $n$)
but also satisfies the weaker condition (3.8) for some constant which is
uniform with respect to $n$. As for the function $H$, by assumption it is
$C^{3}$ and satisfies (2.11) and (2.5)-(2.6) for some $\varpi>0$; hence, up to
replacing $h$ with $h(\sqrt{n^{-1}+t^{2}})$, and up to a further
regularization if needed, we can assume that $H$ is approximated by some
$H^{n}$ which is $C^{3,1}(\mathbb{R}^{d})$ and satisfies (7.8), and in a way
that (2.5)-(2.6) hold uniformly for $H^{n}$. We build so far an approximating
problem

$\left\\{\begin{aligned} &-{\rm
tr}\left({\mathcal{A}}_{n}(x,{\mathcal{D}}u_{n})\,{\mathcal{D}}^{2}u_{n}\right)+\rho
u_{n}+DV_{n}(x)\cdot H_{p}^{n}(Du_{n})=0&&\text{in }Q,\\\
&-\partial_{t}u_{n}+H^{n}(Du_{n})=\delta
u_{n}+f^{\varepsilon}(m_{1n})+\tilde{V}_{n}(x)&&\hbox{at $t=1$,
$x\in\Omega$,}\\\ &-\partial_{t}u_{n}+H^{n}(Du_{n})+\delta
u_{n}=f^{\varepsilon}(m_{0n})+\tilde{V}_{n}(x)&&\hbox{at $t=0$,
$x\in\Omega$,}\\\ &\,\,H_{p}^{n}(Du_{n})\cdot\vec{\nu}=0&&\text{on
$(0,1)\times\partial\Omega$.}\end{aligned}\right.$

which admits a solution $u_{n}\in C^{2,\alpha}(\overline{Q})$ by what we
proved in the first part. This solution $u_{n}$ also belongs to $C^{3}(Q)$ and
we can apply Lemma 3.2 and Theorem 3.4 to infer a uniform estimate for
$\|u_{n}\|_{W^{1,\infty}(Q)}$. By Proposition 7.2, we also deduce that $u_{n}$
is uniformly bounded in $C^{1,\alpha}(\overline{Q})$, because the
$C^{1,\alpha}$ estimate only depends on the bound for
$\|u_{n}\|_{W^{1,\infty}(Q)}$ and on the Lipschitz bounds of the boundary
terms $f^{\varepsilon}(m_{in})$, $\tilde{V}_{n}$. Hence $u_{n}$ is relatively
compact in $C^{1}(\overline{Q})$. By interior regularity, we also have that
$u_{n}$ is bounded in $C^{2,\alpha}(Q)$, so it is relatively compact in
$C^{2}$ if restricted to compact subsets in the interior. This is enough to
pass to the limit and conclude that, up to subsequences, $u_{n}$ converges
towards some $u\in C^{2,\alpha}(Q)\cap C^{1,\alpha}(\overline{Q})$ which
solves problem (3.4) and satisfies estimate (3.18).  $\square$

Acknowledgement. I warmly thank Giuseppe Savaré for stimulating my interest in
this problem and for sharing with me several opinions and hints on the topics
of the paper.

I also thank Filippo Santambrogio for pointing me out reference [33].

## References

  * [1] L. Ambrosio, N. Gigli, G. Savaré, Gradient flows in metric spaces and in the space of probability measures. Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, 2nd ed. (2008).
  * [2] T. Bakaryan, R. Ferreira, D. Gomes, Some estimates for the planning problem with potential, Nonlinear Differ. Equ. Appl. NoDEA (2021).
  * [3] J.-D. Benamou, Y. Brenier, A computational fluid mechanics solution to the Monge- Kantorovich mass transfer problem, Numer. Math., 84 (2000), 375-393.
  * [4] J.-D. Benamou, , G. Carlier, F. Santambrogio, Variational mean field games. Active Particles, Volume 1. Birkhäuser, Cham, (2017), 141-171.
  * [5] P. Cardaliaguet, Weak solutions for first order mean field games with local coupling, in Analysis and geometry in control theory and its applications, Springer 2015, pp. 111–158.
  * [6] P. Cardaliaguet, P. J. Graber, Mean field games systems of first order, ESAIM Control Optim. Calc. Var. 21(2015), 690-722.
  * [7] P. Cardaliaguet, P. J. Graber, A. Porretta, D. Tonon Second order mean field games with degenerate diffusion and local coupling, NoDEA Nonlinear Differential Equations Appl., 22 (2015), 1287-1317.
  * [8] P. Cardaliaguet, A. R. Mészáros, F. Santambrogio, First order mean field games with density constraints: pressure equals price, SIAM J. Control Optim., 54 (2016), 2672-2709.
  * [9] P. Cardaliaguet, A. Porretta, An introduction to Mean Field Game theory, in: Mean field games, 1-158. Lecture Notes in Math., 2281 (CIME Found. Subser.), Springer, Cham, (2020).
  * [10] C. Clason, D.A. Lorenz, H. Mahler, B. Wirth, Entropic regularization of continuous optimal transport problems, arXiv:1906.01333 (2020).
  * [11] E. DiBenedetto: Degenerate parabolic equations, Springer-Verlag, 1993.
  * [12] L. Chizat, G. Peyre, B. Schmitzer, F.-X. Vialard, Unbalanced optimal transport: dynamic and Kantorovich formulations, J. Funct. Anal. 274 (2018), 3090-3123.
  * [13] I. Gentil, C. Léonard, L. Ripani, About the analogy between optimal transport and minimal entropy, Ann. Fac. Sci. Toulouse Math. Série 6 (2017), 569-600.
  * [14] N. Gigli, L. Tamanini, Benamou-Brenier and duality formulas for the entropic cost on $RCD^{*}(K,N)$ spaces, Probability Theory and Related Fields 176 (2020), 1-34.
  * [15] D. Gilbarg, N. Trudinger. Elliptic Partial Differential Equations of Second Order, 2nd ed. Springer, Berlin (1983).
  * [16] D. Gomes, T. Seneci, Displacement convexity for first-order mean-field games, Minimax Theory and Appl. 3 (2018), 261-284.
  * [17] P. J. Graber, A. R. Mészáros, F. Silva, D. Tonon, The planning problem in mean field games as regularized mass transport, Calc. Var. Partial Differential Equations 58 (2019).
  * [18] O. A. Ladyženskaja, V.A. Solonnikov, N.N. Ural’ceva: Linear and quasi-linear equations of parabolic type. Translations of Mathematical Monographs, Vol. 23 American Mathematical Society, Providence, R.I. 1967.
  * [19] O. A. Ladyzhenskaya, N. N. Ural’tseva: Linear and quasilinear elliptic equations. Academic Press, New York (1968).
  * [20] J.-M. Lasry and P.-L. Lions, Jeux à champ moyen. II-Horizon fini et contr$\hat{o}$le optimal. Comptes Rendus Mathématique, 343 (2006), 679-684.
  * [21] J.-M. Lasry and P.-L. Lions, Mean field games, Japanese journal of mathematics 2 (2007), 229-260.
  * [22] H. Lavenant, F. Santambrogio, Optimal density evolution with congestion: $L^{\infty}$ bounds via flow interchange techniques and applications to variational Mean Field Games, Comm. P.D.E. 43 (2018), 1761-1802.
  * [23] G. Lieberman, Solvability of quasilinear elliptic equations with nonlinear boundary conditions, Trans. Amer. Math. Soc. 273 (1982), 753-765.
  * [24] G. Lieberman, The nonlinear oblique derivative problem for quasilinear elliptic equations, Nonlinear Anal. T.M.A. 8 (1984), 49-65.
  * [25] G. Lieberman, Boundary regularity for solutions of degenerate elliptic equations, Nonlinear Anal. 12 (1988), 1203-1219.
  * [26] M. Liero, A. Mielke, G. Savaré, Optimal entropy-transport problems and a new Hellinger-Kantorovich distance between positive measures, Invent. Math. 211 (2018), 969-1117.
  * [27] P.-L. Lions, Cours at Collège de France (2009-2010), https://www.college-de-france.fr/site/pierre-louis-lions/course-2009-2010.htm.
  * [28] R. J. McCann, A convexity principle for interacting gases, Adv. Math. 128 (1997), 153-179.
  * [29] N. Mimikos-Stamatopoulos, S. Munoz, Regularity of one-dimensional first-order mean field games and the planning problem, preprint arXiv:2204.06474 (2022).
  * [30] S. Munoz, Classical and weak solutions to local first order mean field games through elliptic regularity, Ann. I. H. Poincaré Anal. Nonlinéaire 39 (2022), 1-39.
  * [31] C. Orrieri, A. Porretta, G. Savaré, A variational approach to the mean field planning problem, J. Funct. Anal. 277 (2019), 1868-1957.
  * [32] F. Santambrogio: Optimal transport for applied mathematicians. Progress in Nonlinear Differential Equations and their applications 87, Birkhäuser, 2015.
  * [33] F. Santambrogio, X.-J. Wang, Convexity of the support of the displacement interpolation: Counterexamples, Applied Math. Letters 58 (2015), 152-158.
  * [34] G. Stampacchia, Le problème de Dirichlet pour les équations elliptiques du seconde ordre à coefficients discontinus Ann. Inst. Fourier (Grenoble) 15 (1965), 189–258.
  * [35] C. Villani: Topics in Optimal transportation. Number 58 in Graduate Studies in Mathematics. American Mathematical Soc., 2003.


11institutetext: Max-Planck-Institute for Plasma Physics, Garching 85748,
Germany 22institutetext: CS Department, Technical University of Munich,
Garching 85748, Germany 22email: {yehor.yudin, david.coster, udo.v.toussaint,
<EMAIL_ADDRESS>
# Epistemic and Aleatoric Uncertainty Quantification and Surrogate Modelling
in High-Performance Multiscale Plasma Physics Simulations††thanks: Supported
by organisation “Munich School of Data Science”

Yehor Yudin 1122 0000-0001-8867-32073 David Coster 11 0000-0002-2470-9706 Udo
von Toussaint 11 0000-0002-8867-1014 Frank Jenko 1122 0000-0001-6686-1469

###### Abstract

This work suggests several methods of uncertainty treatment in multiscale
modelling and describes their application to a system of coupled turbulent
transport simulations of a tokamak plasma. We propose a method to quantify the
usually aleatoric uncertainty of a system in a quasistationary state,
estimating the mean values and their errors for quantities of interest, which
is average heat fluxes in the case of turbulence simulations. The method
defines the stationarity of the system and suggests a way to balance the
computational cost of simulation and the accuracy of estimation. This allows,
contrary to many approaches, to incorporate of aleatoric uncertainties in the
analysis of the model and to have a quantifiable decision for simulation
runtime. Furthermore, the paper describes methods for quantifying the
epistemic uncertainty of a model and the results of such a procedure for
turbulence simulations, identifying the model’s sensitivity to particular
input parameters and sensitivity to uncertainties in total. Finally, we
introduce a surrogate model approach based on Gaussian Process Regression and
present a preliminary result of training and analysing the performance of such
a model based on turbulence simulation data. Such an approach shows a
potential to significantly decrease the computational cost of the uncertainty
propagation for the given model, making it feasible on current HPC systems.

###### Keywords:

Uncertainty Quantification Sensitivity Analysis Plasma Physics Multiscale
Modelling High Performance Computing Surrogate Modelling Gaussian Process
Regression

##  1 Introduction

Thermonuclear fusion is a prospective source of clean and abundant renewable
energy with a potential for worldwide deployment, independent of geography and
politics. To achieve fusion conditions, one needs to bring plasma temperature
and density to sufficiently high values and isolate the plasma from the
external matter for a sufficiently long time. One of the most studied methods
to confine a plasma is to do it in a toroidal magnetic trap. The most
developed type of device for confinement is a tokamak, an axisymmetric
toroidal chamber with magnetic coils. One of the critical phenomena
interesting for designing and analysing such confinement devices is the
distribution of heat and density of the plasma over the minor radius of the
torus. In the latest decades, it was discovered that transport properties of
plasmas are dominated by microscopic dynamics of turbulent nature[7].

One of the challenges of such a study is the multiscale nature of the
occurring phenomena, with properties of the processes happening across the
size of an entire device, of the order of meters, during the confinement
relevant times, of the order of seconds, are influenced and driven by
microscopic processes, of the order of millimetres and microseconds.

One of the techniques in computational fields applied to tackle such
challenges is coupling single-scale models, each solving a subset of equations
describing processes happening on a particular scale and being implemented as
a separate computer code. Such codes can, in principle, be run as a standalone
program, accepting inputs and producing outputs, adhering to defined data
structures.

In this work, we study the simulation of temperature profiles in a tokamak
plasma, where plasma properties on the smallest scales define heat transport,
and the equations for transport, equilibrium and microscopic turbulence are
solved via separate codes, similar to the approach described by Falchetto et
al. [7].

The following practical step interesting for scientific purposes is estimating
uncertainties in the quantities of interest (QoIs) studied on the largest
scale, considering that those uncertainties come from different sources.

To handle uncertainties, a quantity of interest can be treated as a random
function of its parameters, also understood as random variables. In order to
estimate uncertainties of function values using only a finite number of
evaluations, assumptions can be made on the function, and an intermediate
model or metamodel can be created based on a given sample of evaluations. This
surrogate model describes the function’s dependency on a selected subset of
parameters and serves as a replacement for evaluations of the original
function for practical purposes.

In this work, we apply surrogate modelling to the turbulent transport plasma
simulation, replacing the solution of the microscopic turbulence equation with
a fast model able to infer mean values and uncertainties of turbulent heat
flux.

##  2 Methodology

###  2.1 Epistemic Uncertainty

One of the most common concerns in science is the uncertainty of inference or
predictions due to a lack of knowledge of actual values of specific parameters
or lack of knowledge incorporated into the model, in other words, epistemic
uncertainty.

In theoretical fields, particularly in computational areas, mathematical
models often involve unknown parameters that can only be estimated
probabilistically based on prior knowledge, observations or experiments. The
uncertainty propagation problem means finding an approximation of the
probability density function of quantities of interest in the model, given a
set of uncertain parameters that are understood as random variables with known
probability density functions.

As it is often difficult to retrieve a particular parametric form of QoI PDFs,
the useful statistics describing the properties of output random variables are
the statistical moments of the density.

####  2.1.1 Polynomial Chaos Expansion

One of the most commonly used methods for uncertainty quantification is
Polynomial Chaos Expansion (PCE)[18]. The core idea of the method is to expand
the function of output $Y$ of inputs $X$ using the finite basis of orthonormal
polynomials $P_{i}(x)$ up to the order of $p$, each depending on a subset of
input components and being orthogonal to their PDFs:

$Y\approx\hat{Y}(X)=\sum_{i=0}^{N}c_{i}P_{i}(X)$ (1)

For the function $Y=f_{X}(x)$, the expansion coefficient $\omega_{k}$ and the
ordinate values $x_{k}$ could be determined via a quadrature scheme based on a
Spectral Projection method exploiting the orthogonality of the basis
polynomials and using normalisation factors $H_{i}$ for polynomials $P_{i}$:

$c_{i}=\frac{1}{H_{i}}\int_{\Omega}Y(x)P_{i}(x)f_{X}(x)dx\approx\frac{1}{H_{i}}\sum_{k=1}^{N}Y(x_{k})P_{i}(x_{k})\omega_{k}$
(2)

The benefit of such an orthonormal approach is that the expansion coefficients
are sufficient to calculate the estimates of the moments of the output PDF.

$\hat{\mathbb{E}}[Y]=c_{1},\quad\hat{\mathbb{V}}[Y]=\sum_{i=2}^{N}\mathbb{E}[P_{i}^{2}(X)]c_{i}$
(3)

However, this method suffers from the curse of dimensionality as the number of
required samples grows exponentially, like a binomial coefficient, with the
number of dimensions $d:N=(d+p)!/p!\cdot d!$

####  2.1.2 Sensitivity Analysis

A particular type of statistics over a sample of code’s QoI responses aims to
define the contribution of a single parameter, or a subset of parameters, to
the uncertainty of the output. In this work, we study Sobol indices[17], which
is a global variance-based sensitivity metric. The quantity we estimate is the
variance of the mathematical expectations of the function output, conditioned
on a subset of input components $\\{i\\}$ and normalised by the total variance
of the output. Here we use the first-order Sobol indices, identifying the
influence of a single parameter $i$; higher-order Sobol indices, determining
the influence of a subset of interacting parameters $\\{i\\}$; and total Sobol
indices, identifying the total effect of each parameter, using a set of all
components $-i$, excluding a single $i$ component:

$S_{\\{i\\}}=\frac{\mathbb{V}[\mathbb{E}[Y|X_{\\{i\\}}]]}{\mathbb{V}[Y]}\,,\quad
S_{i}^{tot}=1-\frac{\mathbb{V}[\mathbb{E}[Y|X_{-i}]]}{\mathbb{V}[Y]}$ (4)

Using such an indicator, one could judge which parameters are most important
for the behaviour of the computational model and decide whether a particular
parameter’s value can be fixed in subsequent studies or whether more values
should be tried out.

###  2.2 Aleatoric Uncertainty

One of the types of uncertainties, different from the epistemic, is the
aleatoric uncertainty coming from the inherent properties of the model and its
solver. In the first place, it is related to the dynamic stochastic behaviour
of the quantities of interests, for which the model is solved.

Several methods for analysing aleatoric uncertainties commonly involve two
principles: separating epistemic uncertainties and performing statistical
analyses of output quantities. One approach uses replica sampling, in which a
certain number of model replicas are solved for each set of input parameters,
followed by analysing the distributions of quantities of interest values[19].
Other methods are suitable for analysing the quasistationary nonlinear
behaviour of the model in the vicinity of a certain attractor over a
considerable time of solution[2].

This work utilises a practical measure for an aleatoric uncertainty based on
analysis of the standard error of the mean values of quantities of interest
$E=\mathrm{SEM}[y]=\sigma[y]/\sqrt{n_{\mathrm{eff}}}$ for a sample of
observations $\mathbf{Y}^{\mathrm{eff}}$, each taken in an autocorrelation
time window of length $t_{\mathrm{A}}$ apart from the previous one. Such an
approach allows judging the requirements of a single solution for a given set
of parameters and controlling the level of aleatoric uncertainty. However,
such a measure only provides information on the second moments of the quantity
of interest distribution. Furthermore, it is not always well incorporable into
a single framework of analysis of both epistemic and aleatoric uncertainties.

###  2.3 Surrogate Modelling

This work primarily focuses on applying Gaussian Process Regression (GPR), a
non-parametric probabilistic machine learning model fully defined by the
covariance matrix of observed data samples and capable of providing a
posterior PDF of output given any input values from the considered
support[15].

When the Matern-$3/2$ kernel $k$ describes covariance between any two
parametric points $x_{i},x_{j}$, determining a PDF $p(f(x_{\ast}))$ of a
posterior for new data $x_{\ast}$ using a standard GPR based on a sample of
observed inputs $X$ and outputs $y$ with a Normal likelihood looks like
following:

$\displaystyle
p(f(x_{\ast})|\mathbf{X},\mathbf{y},x_{\ast})\sim\mathcal{N}(\mu(x_{\ast}),\sigma^{2}(x_{\ast}))$
(5)
$\displaystyle\sigma^{2}(x_{\ast})=K(x_{\ast},x_{\ast})-K(x_{\ast},\mathbf{X})^{\top}\left(K(\mathbf{X},\mathbf{X})+\sigma_{e}^{2}\mathbf{(}I)\right)^{-1}K(x_{\ast},\mathbf{X})$
(6) $\displaystyle K=\left(k(x_{i},x_{j})\right)_{i,j=1..N}$ (7)
$\displaystyle
k(x_{i},x_{j})=\sigma^{2}\left(1+\sqrt{3}\frac{|x_{i}-x_{j}|}{l}\right)\exp\left(-\sqrt{3}\frac{|x_{i}-x_{j}|}{l}\right)$
(8)

The method allows for capturing uncertainties in the underlying data that come
from the lack of information on the parametric dependencies and internal
noise. Furthermore, this method has a high prediction uncertainty in the
region of the parametric space for which there was not enough training data,
as the posterior would have a higher variance for points decorrelated with the
observed sample. Such a high uncertainty indicates the untrustworthiness of
the regression results and the model’s limitations.

One of the ideas behind the surrogate modelling approach is to create flexible
and adaptive surrogates based on new data. Online Learning allows adding new
data into the training sample to improve the surrogate gradually used as a
proxy solver[10]. Active Learning uses the regression model’s prediction
uncertainty and information about sought dependency to determine new points in
the parameter space for the evaluation of the function of interest[14].

In this work, we apply Gaussian Process Regression trained on a data set
describing performed turbulence simulations. We trained such a surrogate model
for the dependency of heat fluxes on the temperature values and its gradient
for a particular configurational location in a tokamak plasma.

##  3 Numerical Results

###  3.1 Computational Model

The object of the study is the effects of turbulent processes happening at
microscopic time and space scales on the temperature profiles of a nuclear
fusion device across the confinement scales. For that, we use the Multiscale
Fusion Workflow (MFW)[11] consisting of three main models, each implemented as
a separate computer code, serving as an independent workflow component capable
of exchanging its solution with others in a black-box fashion. The models are
the equilibrium code that describes plasma geometry, the transport code that
evolves the plasma profiles and the turbulence code used to compute effective
transport coefficients (Fig 1). The code that solves equilibrium Grad-
Shafranov equations for plasma density and the magnetic field is CHEASE[12], a
high-resolution fixed-boundary solver. Evolution of temperature and density is
performed with ETS[3], a 1D code for energy and particle transport.

Figure 1: Main components of the Multiscale Fusion Workflow. The turbulence
model is the one that requires the most computational resources to solve, and
that exposes the highest sensitivity to its inputs.

In this work, we analyse the turbulence component of the workflow, for which
we use the 3D nonlinear gyrofluid turbulence code GEM[16] in its local flux
tube approximation. Such turbulence models are very sensitive to input
parameters and computationally expensive to solve. Thus, it is a prime point
for performing a study of uncertainties in the joint multiscale model and the
component for which a fast and easy-to-evaluate surrogate model would be of
interest to create. As discussed in [4], performing UQ on the MFW with GEM is
currently computationally unreasonable since 100 000 core hours are required
for a single run, whereas a surrogate model should be at least 10 000 times
cheaper (1 core instead of 1024, and a factor of 10 per cycle around the MFW
loop).

Codes store their respective solution in standardised data structures of
Consistent Physical Objects (CPO)[8], and the turbulence quantities of
interests are stored in a structure describing core transport properties. In
principle, the workflow implemented to analyse the flux time traces and
prepare a dataset for surrogate training and utilisation can support any code
capable of interfacing via CPOs.

The resulting data structures are passed from a code to a code instance in a
point-to-point fashion using the MUSCLE2 coupling library[1]. Currently, a
version supporting various turbulence models is being implemented using the
MUSCLE3 library[20].

###  3.2 Simulations and resulting time traces

A number of simulations of turbulence code for different input parameter
values were performed to study the behaviour of quantities of interest in
turbulence. Furthermore, this data was used to analyse what information can be
included in a surrogate model and as an initial data set for the model
training.

Transport equations describe temperature and density behaviour using a 1D
approach, where quantities are functions of radial coordinate and time but
must be solved globally for the entire domain. Turbulence occurs on a smaller
scale, and behaviour around a single magnetic flux tube or a small range of
the radial domain typically defines it very well. In order to interpolate the
turbulence solution for the entire device, the tokamak turbulence model for
the ASDEX Upgrade[7] must be solved for at least eight different radial
positions[4].

The turbulent model solution’s quantity of interest is the heat flux averaged
across a magnetic field surface at a specific radial position in the toroidal
device. The computed result is nonlinear on microsecond timescales, with
several features and stages typical for all relevant parameter space.
Initially non-turbulent, the effective QoI value grows exponentially as
turbulence develops until nonlinear effects dominate, causing turbulence to
saturate and QoI to behave quasistationary. Although mean values remain around
a particular mean value, the QoI is also locally chaotic and fluctuates
substantially over time.

Defining whether the system is saturated, accurately representing the mean
output values for a particular case, and decoupling fluctuating behaviour are
crucial to creating an effective surrogate model parametrising the transport
properties.

Macroscopic kinetic profiles behave similarly to micro-level turbulent fluxes
over confinement-relevant times. Starting from initial conditions, they
converge towards a stationary level self-consistent with the coupled model.
The profiles dynamically and nonlinearly interact with local turbulent
behaviour, resulting in a quasistationary solution. In order to capture the
dependency of turbulent fluxes in relevant parameter space, a surrogate model
must be based on a large sample of code solutions covering the relevant
regions.

Creating a sample for surrogate training involves running the same code for
multiple input parameter values and managing the simulation and training data,
which pose a significant challenge. Using the resulting surrogate as a
microscopic code replacement requires ensuring it captures QoI dependency well
for all practically possible parameter values. This requires a training
dataset covering enough parameter space information; the surrogate should
indicate its applicability limits and epistemic uncertainty due to a lack of
training data.

Given the practical challenges, this work suggests another novel computational
workflow besides one of the multiscale simulations to solve the turbulent
transport problem. The new workflow analyses turbulent fluxes time traces on
the microscopic temporal scales and defines the stationarity of the QoI as
well as its mean values and estimation error.

###  3.3 Aleatoric uncertainty

One of the work’s challenges was to separate fundamental irreducible
uncertainty from epistemic parametric uncertainty in a fluctuating system. It
was important to establish the stationary mean values for analysing key
parametric dependencies and to use fluctuations only for statistical error
estimation. However, solving the model for a long enough time to gather
sufficient data is computationally demanding and gathering a large sample
would be prohibitively expensive. Estimating error levels during the course of
a simulation run can help decide when to stop the simulation and significantly
save computational costs.

The procedure for analysing the QoI time traces $y(t)$, representing a model
solution for a single parametric point, constitutes of the following steps:

  1. 1.

For the time-traces $y(t)$ of length $t_{n}$ representing a model solution in
a scalar quantity of interest at time steps $\\{t\\}$, we select a part in the
saturated phase:

    1. (a)

Define the $y(t)$ ramp-up phase duration $t_{\mathrm{r.u.}}$: here, in
practice, for the long term, we chose an initial $15\%$ of readings

    2. (b)

Discard the readings from the ramp-up phase

  2. 2.

Downsample the readings:

    1. (a)

Calculate the Auto-Correlation Time:
$t_{\mathrm{A}}=\mathrm{ACT}[y(t)]=\mathrm{min}\,t^{\ast}:\\\
\frac{1}{(t_{n}-t^{\ast})}\sum_{t=t^{\ast}}^{t_{n}}\left(y(t)-\bar{y}\right)\cdot\left(y(t-t^{\ast}\right)-\bar{y})<\frac{1}{t_{n}}\sum_{t=t_{1}}^{t_{n}}y^{2}(t)\cdot
e^{-1}$

    2. (b)

Split time series in saturated phase into
$n_{\mathrm{eff}}=\lfloor\frac{t_{n}}{t_{\mathrm{A}}}\rfloor$ windows

    3. (c)

For every autocorrelation time window of length $t_{\mathrm{A}}$ choose a
downsampled reading as a mean value
$y_{i}^{\mathrm{eff}}=\frac{1}{t_{\mathrm{A}}}\sum_{j=i\cdot
t_{\mathrm{A}}}^{(i+1)\cdot t_{\mathrm{A}}}y_{j}$ for an effective time step
$t_{i}^{\mathrm{eff}}=\frac{1}{t_{\mathrm{A}}}\sum_{j=i\cdot
t_{\mathrm{A}}}^{(i+1)\cdot t_{\mathrm{A}}}j$

    4. (d)

Collect downsampled readings into a new set
$\mathbf{Y}^{\mathrm{eff}}=\\{y_{i}^{\mathrm{eff}}\\}$

  3. 3.

Test the stationarity of the resulting time series:

    1. (a)

Here: compute an ordinary least-squares linear multivariate regression model
of downsampled QoI readings over the time

    2. (b)

Apply Normal Equations to find coefficients:
$\hat{\mathbf{\theta}}=\left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{Y}$
where $\mathbf{X}$ consists of effective time steps $t_{i}^{\mathrm{eff}}$

    3. (c)

Test if the linear regression coefficients are below a chosen relative
tolerance: $\hat{\mathbf{\theta}}<\mathbf{\epsilon}_{\mathrm{tol}}$

      1. i.

If the coefficients are too large, continue running the simulation for another
$t_{\mathrm{run}}$ time steps

      2. ii.

If the coefficients are small enough, stop the simulation and proceed to the
statistics calculation

  4. 4.

Compute the essential statistics estimates:

    1. (a)

Mean:
$\mu[y]=\frac{1}{n_{\mathrm{eff}}^{\mathrm{}}}\,\sum_{y\in\mathbf{Y}^{\mathrm{eff}}}y$

    2. (b)

Standard deviation:
$\sigma[y]=\left(\frac{1}{n_{\mathrm{eff}}^{\mathrm{}}-1}\,\sum_{y\in\mathbf{Y}^{\mathrm{eff}}}(y-\mu[y])^{2}\right)^{-\frac{1}{2}}$

    3. (c)

Standard Error of the Mean:
$\mathrm{SEM}[y]=\sigma[y]/\sqrt{n_{\mathrm{eff}}^{\mathrm{}}}$

We applied this procedure for the time traces produced by turbulence code
iteratively, in particular to ion heat flux values over time, sampling after a
number of solution time steps, and analysed the behaviour of the mean estimate
and its error.

Figure 2 describes mean values calculated with the sequential procedure based
on the data of time traces of heat fluxes in a quasistationary phase.

Figure 2: The standard error of the mean is calculated accounting to
discarding 15% of reading that do not belong to stationary phase, calculating
the autocorrelation time of the remaining time traces, taken a sub-sample
constituent of a single mean reading per autocorrelation window and
renormalising mean error using the effective size of the sample.

After the relative change of the mean and the relative change of the standard
error of the mean converge below a certain threshold, as shown in Figure 3, we
considered the duration of a simulation sufficient enough to have a reasonable
estimate of the average flux and its error and added this data to a dataset
describing the computational turbulence model behaviour. Subsequently, this
dataset is used to draw conclusions on the model’s parametric uncertainties
and train a surrogate model.

Figure 3: Sequential analysis of time traces with an interval of 150 time-
steps per window, defined by 6 hours of computation on a single MPCDF’s COBRA
node. After each step of time traces computation, we recalculate the
autocorrelation time, repopulate an effective data sample, and estimate new
values for the mean, standard deviation and standard error of the mean. With
the growth of the sample, the mean estimation converges to the true value, and
the standard error of the mean decreases to around $1\%$ of the respective
mean value after $2300$ time steps out of $4000$ in total. Here, one can
decide to stop the computation using a tolerance threshold on the relative
change of the mean estimate of $10^{-3}$ or by the error of the mean dropping
below $1.1\%$ Coefficient of Variation. The plots of the standard error of the
mean and autocorrelation time dependency on number simulation steps indicate a
convergence of the time traces analysis with simulation time.

###  3.4 Epistemic Uncertainty

To explore the epistemic uncertainties of the properties of the turbulence
model, we performed an uncertainty propagation, taking a particular input of
macroscopic parameter with a value taken from experimental measurement and,
assuming that on top of this value, there is uncertainty, described as a
normally distributed random variable with a Coefficient of Variation, a ratio
of the variance of a quantity to its mean value,
$\mathrm{CV}(x)=\sqrt{\mathbb{V}[x]}/\mathbb{E}[x]$ equal to $0.1$. The
statistics for average ion and electron heat flux
$\mathbf{y}=(Q_{i}(\rho),Q_{e}(\rho))$ were estimated using a PCE method with
Hermite polynomials of order $p=2$ applied for four input parameters of ion
and electron temperatures profile values and gradients:
$\mathbf{x}=T_{i}(\rho),T_{e}(\rho),\nabla T_{i}(\rho),\nabla T_{e}(\rho)$.

The uncertainty propagation was performed using EasyVVUQ[9] library, which
manages the definition of input uncertainties, sampling scheme, generation of
input data structures and collection of the output ones, as well as
calculating statistics. The batches of the turbulence code instances were run
on the MPCDF COBRA supercomputer using QCG-PilotJob[13] middle-ware library,
which manages HPC resources of a single computational allocation and
distributes it to run multiple instances of the code with a particular
parallelisation set-up.

One of the issues of uncertainty propagation through the model that exposes
chaotic behaviour on some scales is to untangle the epistemic and aleatoric
uncertainties. Before analysing the parametric uncertainty propagation
results, the time traces produced by the simulations are processed as
described in the previous subsection.

The result of this study with an error of $\mathrm{CV(\mathbf{x})}=0.1$ was
inconclusive, as there was a significant overlap of the error bounds for mean
values of quantities of interest for neighbouring points in parametric space
(Fig. 4). In this case, it was impossible to conclude with certainty whether
the ion heat flux trend was rising or falling with respect to the growth of
ion temperature and its gradient.

Figure 4: Dependency of the mean ion heat flux on the ion temperature. For the
small variation of input parameters, the mean estimation could not allow for
drawing good conclusions on the gradient of this dependency, and the means are
within uncertainty intervals of neighbouring points. For the larger variation
of input parameters, the relative errors of the mean, also stated next to the
points, are small enough to reconstruct the dependency with sufficient
certainty.

As the next step, we have taken an input parameter with normal distribution
with $\mathrm{CV(\mathbf{x})}=0.25$ and propagated this uncertainty through
the model. Enlarging the input variation allowed us to study the behaviour and
uncertainties of the model over the broader range of input parameter values
and recover the character of the dependency of ion heat flux on ion
temperature and its gradient (Fig. 4). Furthermore, longer simulation runs
allowed reaching lower levels of the aleatoric uncertainty estimates in the
quantities of interest,

The uncertainty analysis results per input component show a dominant parameter
for which the model is most sensitive (Fig. 5). For this particular prior PDF
over parametric space, the behaviour of the electron heat flux is most
sensitive to variation in electron temperature. Still, other parameters and
their interactions also play a significant role.

Figure 5: Sobol indices for electron and ion heat fluxes. Most of the
variation in electron heat flux is produced by the variation in the electron
temperature. Interactions of the input parameters explain some of the non-
negligible variation. For the ion heat flux, the output variation is dominated
by a single parameter in this regime. However, a fraction is still explained
by the nonlinear interaction of parameters.

Other metrics that were chosen to analyse the model’s epistemic uncertainty
are the Coefficient of Variation in the quantity of interest and the
Uncertainty Amplification Factor[5]. In this work, the coefficient of
variation of the model’s output ion heat flux is about
$\mathrm{CV}(Q_{i})\approx 2.96$, showing a significant range of possible
values that QoI can take in the study.

The Uncertainty Amplification Factor is the ratio of the Coefficient of
Variation of a QoI to the coefficient of variation of the input parameter
$\mathrm{UAF}(y|x)=\mathrm{CV}(y)/\mathrm{CV}(x)$. In the case of the studied
model and chosen input parameters prior distribution, the Uncertainty
Amplification Factor of ion heat flux is
$\mathrm{UAF}(Q_{i}|\mathbf{x})\approx 11.9$ showing a significant sensitivity
of the model to the parametric uncertainty.

###  3.5 Surrogate model

Having performed the original epistemic uncertainty propagation, we have
collected a substantial sample of the physical model’s input-output pairs.
This sample was used to fit a surrogate model able to capture the behaviour of
the physical model QoI means values and uncertainties over the selected input
parametric subspace. Here, the QoIs and the input parameters are the same as
described in the epistemic uncertainty propagation results section.

Figure 6: GPR surrogate model test results with the coefficient of
determination $R^{2}=0.988$. Inferred values of ion heat flux compared to the
observed values produced by the physical model. Vertical bars denote the $\pm
1.96$ of the standard deviation of the predictive posterior distribution, grey
bounds denote 15% relative error, and dotted lines denote the region of
interest. For most test samples, the true value is within $95\%$ prediction
interval, and the predicted mean is within $15\%$ relative error.

The GPR surrogate provides an intermediate model based on the assumptions on
the regularity of the actual dependency, primarily encoded in the covariance
kernel and the likelihood of the regression model. Unlike the PCE surrogate,
which depends on the predefined polynomial order and is difficult to adapt
given new information about the sought dependency, the uncertainties of the
GPR one produce strong indications that the information on the actual
dependency is sufficient for the regression and is easily adaptable through
sequential design and active learning algorithms.

Here we present the results of the Gaussian Process Regression surrogate
trained on $50\%$ of the collected simulation dataset (Fig. 6). We developed
functionality for the training and analysis of surrogates based on
EasySurrogate[6] library. After performing a hyper-parameter optimisation, a
particular set of categorical parameter values was chosen to judge surrogates
by the highest coefficient of determination of QoI values predicted by a
surrogate $y_{i}^{\mathrm{surr}}$ compared to the ground truth simulation
values $y_{i}^{\mathrm{g.t.}}$ for the model testing:
$R^{2}=1-{\sum_{i=1}^{n_{\mathrm{test}}}\left(y_{i}^{\mathrm{g.t.}}-y_{i}^{\mathrm{surr}}\right)^{2}}/{\sum_{i=1}^{n_{\mathrm{test}}}\left(y_{i}^{\mathrm{g.t.}}-{\overline{y^{\mathrm{g.t.}}}}\right)^{2}}$

The future plan is to apply a validation test for the applicability of such a
surrogate as a replacement for the solution of the turbulence model for given
inputs describing core profiles and inferring the average heat fluxes and
effective transport coefficients. The original simulation workflow utilising
the turbulence code GEM required around $27.0$ s on $1024$ cores of Cineca
MARCONI supercomputer[11] to predict fluxes for the next iteration of core
profiles advancement, whether the EasySurrogate GPR model requires a fraction
of a second on a single core for the same procedure. Furthermore, applying a
GPR model would allow the detection of points in parameter space where the
surrogate uncertainty is too high and should not be applied as a proxy for a
high-fidelity solver. When it is understood that a GPR surrogate performs sub-
par for a region of interest in parametric space, some work was done on
designing an Active Learning scheme based on Bayesian Optimisation and
Sequential Design of Experiments.

##  4 Discussion

In this paper, we presented a number of methods to treat uncertainties in
multiscale simulations. This includes ways to measure aleatoric and epistemic
uncertainties, to control and reduce the amount of computational resources
required to quantification of the uncertainties, as well as to identify
situations when uncertainties are hindering analysis of the physical model
behaviour. The suggested method for aleatoric uncertainty processing is
applicable when some model dependency could be decomposed into the mean and
fluctuation behaviour and allows for a trade-off between the computational
effort for the model solution and the accuracy of the mean behaviour
estimation. The surrogate modelling method proposes a way to make epistemic
uncertainty propagation cheaper by constructing a data-based machine learning
GPR proxy model for an expensive physical model solver that could be used as a
replacement in multiscale coupled simulations and for sampling in uncertainty
quantification methods.

We applied the described methods for the dependency of average heat fluxes
produced by plasma turbulence code GEM based on tokamak core temperature
profiles and their gradients. While studying such fluxes’ time traces, we
identified cases when the simulations converged early enough to allow reducing
the duration of runs and saving computational resources. The QoI mean levels
were estimated with sufficient certainty to draw conclusions about sought
dependency and use solution data as a basis for further analysis and surrogate
model training. The subsequently performed uncertainty propagation and
sensitivity analysis showed situations where aleatoric uncertainties could not
allow concluding on parametric dependencies. Also, we defined cases with a
strong indication that a particular set of variables influences output
variation dominantly. The surrogate modelling approach of fitting a Gaussian
Process Regression model to the dataset of flux computing turbulent
simulations showed that it is possible to train a metamodel that allows
capturing the mean behaviour of the fluxes well within the posterior
predictive bounds for a region of interest in parameter space of the model.

The development of a surrogate for the turbulence code opens up the much-
needed possibility of doing detailed UQ for the whole MFW workflow with the
accuracy of the underlying turbulence code but without the expense of the
turbulence code in the workflow.

Future work would include testing the surrogate as a proxy of the expensive
turbulence code, capable of inferring fluxes and subsequently defining the
effective transport coefficients. This includes implementing an expanded
version of a Multiscale Fusion Workflow, capable of switching between
different implementations of the turbulence model and validating the results
of the coupled turbulent transport simulations using a turbulence surrogate
against a version using the gyrofluid code. Subsequently, such a modified
workflow should include capacity for adaptive surrogate training through
Online Learning and Active Learning algorithms, as well as methods to select
turbulence model implementation in the course of a single simulation.

##  5 Acknowledgements

The authors of this paper would like to acknowledge the support of the Poznan
Supercomputer and Networking Centre (PSNC) and Max-Planck Computational Data
Facility (MPCDF). Research by Yehor Yudin is funded by the Helmholtz
Association under the ”Munich School for Data Science - MuDS”.

## References

  * [1] Borgdorff, J., et al.: Distributed multiscale computing with MUSCLE 2, the Multiscale Coupling Library and Environment. Journal of Computational Science 5(5), 719–731 (Sep 2014). https://doi.org/10.1016/j.jocs.2014.04.004
  * [2] Brajard, J., et al.: Combining data assimilation and machine learning to infer unresolved scale parametrization. Phil. Trans. of the Roy. Soc. A: Math., Phys. and Eng. Sc. 379(2194) (2021). https://doi.org/10.1098/rsta.2020.0086
  * [3] Coster, D., et al.: The European Transport Solver. IEEE Transactions on Plasma Science 38(9), 2085–2092 (Sep 2010). https://doi.org/10.1109/TPS.2010.2056707
  * [4] Coster, D., et al.: Building a turbulence-transport workflow incorporating uncertainty quantification for predicting core profiles in a tokamak plasma. Nuclear Fusion (2021)
  * [5] Edeling, W., et al.: The impact of uncertainty on predictions of the CovidSim epidemiological code. Nature Computational Science 1(2), 128–135 (Feb 2021). https://doi.org/10.1038/s43588-021-00028-9
  * [6] Edeling, W.: Easysurrogate github page, https://github.com/wedeling/EasySurrogate, last accessed 2 Mar 2023
  * [7] Falchetto, G., et al.: The European Integrated Tokamak Modelling (ITM) effort: achievements and first physics results. Nuclear Fusion 54(4), 043018 (Apr 2014). https://doi.org/10.1088/0029-5515/54/4/043018, 00005
  * [8] Imbeaux, F., et al.: A generic data structure for integrated modelling of tokamak physics and subsystems. Computer Physics Communications 181(6), 987–998 (Jun 2010). https://doi.org/10.1016/j.cpc.2010.02.001
  * [9] Jancauskas, V., Lakhlili, J., Richardson, R., Wright, D.: EasyVVUQ: Verification, validation and uncertainty quantification for HPC simulations (2021), https://github.com/UCL-CCS/EasyVVUQ
  * [10] Leiter, K., Barnes, B., Becker, R., Knap, J.: Accelerated scale-bridging through adaptive surrogate model evaluation. Journal of Computational Science 27, 91–106 (2018). https://doi.org/0.1016/j.jocs.2018.04.010
  * [11] Luk, O., et al.: ComPat framework for multiscale simulations applied to fusion plasmas. Computer Physics Communications 239, 126–133 (Jun 2019). https://doi.org/https://doi.org/10.1016/j.cpc.2018.12.021
  * [12] Lütjens, H., Bondeson, A., Sauter, O.: The chease code for toroidal mhd equilibria. Computer Physics Communications 97(3), 219–260 (1996)
  * [13] Piontek, T., et al.: Development of science gateways using qcg – lessons learned from the deployment on large scale distributed and hpc infrastructures. J. of Grid Comp. 14(4), 559–573 (2016). https://doi.org/10.1007/s10723-016-9384-9
  * [14] Preuss, R., von Toussaint, U.: Global optimization employing gaussian process-based bayesian surrogates. Entropy 20 (2018). https://doi.org/10.3390/e20030201
  * [15] Rasmussen, C.E., W.K.: Gaussian Processes for Machine Learning,. the MIT Press (2006). https://doi.org/10.7551/mitpress/3206.001.0001
  * [16] Scott, B.D.: Free-energy conservation in local gyrofluid models. Physics of Plasmas 12(10), 102307 (Oct 2005). https://doi.org/10.1063/1.2064968
  * [17] Sobol, I.: On sensitivity estimation for nonlinear mathematical models. Mathematical Modeling 2(1), 112–118 (1990)
  * [18] Sullivan, T.: Introduction to Uncertainty Quantification, TAM, vol. 63. Springer Cham (2015). https://doi.org/10.1007/978-3-319-23395-6
  * [19] Vassaux, M., et al.: Ensembles are required to handle aleatoric and parametric uncertainty in molecular dynamics simulation. Journ. of Chem. Th. and Comp. 17(8) (2021). https://doi.org/10.1021/acs.jctc.1c00526
  * [20] Veen, L., Hoekstra, A.: Easing multiscale model design and coupling with MUSCLE 3. In: Computational Science – ICCS 2020. pp. 425–438. Springer International Publishing, Cham (2020). https://doi.org/10.1007/978-3-030-50433-5-33