# Discrepancies and their means

M.M. SKRIGANOV  St. Petersburg Department of the Steklov Mathematical
Institute of the Russian Academy of Sciences, 27, Fontanka, St.Petersburg
191023, Russia<EMAIL_ADDRESS>
###### Abstract.

It is shown that the discrepancy function for point distributions on a torus
is expressed by an explicit formula in terms of its mean values on sub-tori.
As an application of this formula, a simple proof of a theorem of Lev [2] on
the equivalence of $L_{\infty}$– and shifted $L_{q}$–discrepancies is given.

###### Key words and phrases:

Point distributions, discrepancy theory

###### 2010 Mathematics Subject Classification:

11K38

The point distribution problem on the $d$-dimensional torus
${\mathbb{T}}^{d}={\mathbb{R}}^{d}/{\mathbb{Z}}^{d}$ is conveniently
considered as a periodic problem on the covering space ${\mathbb{R}}^{d}$. For
$X=(x_{1},\dots,x_{d})\in{\mathbb{R}}^{d}$ and
$Y=(y_{1},\dots,y_{d})\in{\mathbb{K}}^{d}=[0,1]^{d}$, we define the periodic
discrepancy function by

$\displaystyle L(X,Y)=\chi(X,Y)-v(Y),$ (1)

where
$\chi(X,Y)=\prod\nolimits_{j=1}^{d}\chi(x_{j},y_{j}),\,\,v(Y)=\prod\nolimits_{j=1}^{d}v(y_{j}),v(y)=y,$
and

$\displaystyle\chi(x,y)=\begin{cases}1,&\text{if $\\{x\\}<y$,}\\\
0,&\text{otherwise,}\end{cases}$ (2)

here $\\{x\\}$ is the fractional part of $x\in{\mathbb{R}}$, and $y\in[0,1]$.
It is clear that $\chi(X,Y)$ is the indicator function of the periodic
collection of rectangular boxes

$\displaystyle{\mathscr{B}}(Y)=\bigcup\nolimits_{(m_{1},\dots,m_{d})\in{\mathbb{Z}}^{d}}\prod\nolimits^{d}_{j=1}[m_{j},y_{j}+m_{j}).$

The mean value of the discrepancy function (1) has the form

$\displaystyle\lambda(X)=\int\nolimits_{{\mathbb{K}}^{d}}L(X,Y)\,\mathrm{d}Y=\prod\nolimits_{j=1}^{d}(1-\\{x_{j}\\})-2^{-d}.$
(3)

Let $[d]=(1,\dots,d)$ denote the set of coordinate indexes. For subsets
$J\subseteq[d]$ we introduce the partial discrepancy functions by

$\displaystyle L_{J}(X,Y)=\chi_{J}(X,Y)-v_{J}(Y),$ (4)

where $\chi_{J}(X,Y)=\prod\nolimits_{j\in
J}\chi(x_{j},y_{j}),\,\,v_{J}(Y)=\prod\nolimits_{j\in J}v(y_{j}),$ The
corresponding mean values have the form

$\displaystyle\lambda_{J}(X)=\int\nolimits_{{\mathbb{T}}^{|J|}}L_{J}(X,Y_{J})\,\mathrm{d}Y_{J}=\prod\nolimits_{j\in
J}(1-\\{x_{j}\\})-2^{-|J|},$ (5)

where ${\mathbb{T}}^{|J|}\subseteq{\mathbb{T}}^{d}$ denotes the sub-torus
corresponding to the subset $J\subseteq[d]$, and $|J|$ denotes the number of
elements of $J$.

The quantities (4), (5) depend on the projections $X_{J}=(x_{j})_{j\in
J}\in{\mathbb{R}}^{|J|},Y_{J}=(y_{j})_{j\in J}\in{\mathbb{K}}^{|J|}$, but not
on the additional variables $X_{J^{\prime}}=(x_{j})_{j\in
J^{\prime}}\in{\mathbb{R}}^{|J^{\prime}|},Y_{J^{\prime}}=(y_{j})_{j\in
J^{\prime}}\in{\mathbb{K}}^{|J^{\prime}|}$, where $J^{\prime}=[d]\setminus J$
denotes the compliment of $J$. For $J=[d]$, we write
$L_{[d]}(X,Y)=L(X,Y),\,\lambda_{[d]}(X)=\lambda(X)$, and for the empty set, we
put $L_{\emptyset}(X,Y)=0,\,\lambda_{\emptyset}(X)=0$.

Definition. For a periodic function
$f(X)=f(X_{J},X_{J^{\prime}}),X\in{\mathbb{R}}^{d}$, and a vector
$Y_{J}=(y_{j})_{j\in J}\in{\mathbb{K}}^{|J|}$, the alternant is defined by

$\displaystyle
f^{(alt)}(X\,|\,Y_{J})=f^{(alt)}(X_{J},X_{J^{\prime}}\,|\,Y_{J})=\sum\nolimits_{\Theta_{J}}(-1)^{|\Theta_{J}|}f(X_{J}-\Theta_{J}\cdot
Y_{J},X_{J^{\prime}}),$ (6)

where $\Theta_{J}=(\theta_{j})_{j\in J}\in\\{0,1\\}^{|J|}$ are the vertexes of
the cube ${\mathbb{K}}^{|J|}$, and summation in (6) is taken over all such
vertexes, $\Theta_{J}\cdot Y_{J}=(\theta_{j}y_{j})_{j\in J}$ and
$|\Theta_{J}|=\sum\nolimits_{j\in J}\theta_{j}$.

Remark. If $f(X)=\prod\nolimits_{j\in J}f_{j}(x_{j})$, then

$\displaystyle f^{(alt)}(X\,|\,Y_{J})=\prod\nolimits_{j\in
J}f^{(alt)}_{j}(x_{j}\,|\,y_{j})=\prod\nolimits_{j\in
J}(f(x_{j})-f(x_{j}-y_{j})),$ (7)

and if at least one of the functions $f_{j}$ is constant, then $f^{(alt)}=0$.

Theorem (Main Identity). The discrepancy function $L(X,Y)$ satisfies the
identity

$\displaystyle
L(X,Y)=\sum\nolimits_{J\subseteq[d]}\,v_{J^{\prime}}(Y)\,\lambda^{(alt)}_{J}(X\,|\,Y_{J}).$
(8)

###### Proof.

Let $\omega(x)=\frac{1}{2}-\\{x\\}$, then
$\int\nolimits^{1}_{0}\omega(x)\mathrm{d}x=0$. We put

$\displaystyle\omega_{J}(X)=\prod\nolimits_{j\in
J}\omega(x_{j})\qquad\text{and}\qquad\omega_{\emptyset}(X)=0.$ (9)

The indicator function (2) can be written in the form

$\displaystyle\begin{split}\chi(x,y)&=y-\\{x\\}+\\{x-y\\}\\\
&=y+\omega(x)-\omega(x-y)=y+\omega^{(alt)}(x\,|\,y),\end{split}$ (10)

This formula can be proved by considering the graph of the function
$\\{x\\}-\\{x-y\\},x\in{\mathbb{R}}$.

Substituting (10) into (1) and using (9) and (7), we obtain

$\displaystyle
L(X,Y)=\sum\nolimits_{J\subseteq[d]}\,v_{J^{\prime}}(Y)\,\omega^{(alt)}_{J}(X\,|\,Y_{J}).$

For the mean value (5), we find

$\displaystyle\lambda_{J}(X)=\prod\nolimits_{j\in
J}(2^{-1}+\omega(x_{j}))-2^{-|J|}=\sum\nolimits_{I\subseteq
J}\,2^{-|J\setminus I|}\,\,\omega_{I}(X).$ (11)

Let us calculate the alternant of the mean value (11). We have

$\displaystyle\lambda^{(alt)}_{J}(X\,|\,Y_{J})=\omega^{(alt)}_{J}(X)+\sum\nolimits_{I\subset
J}\,2^{-|J\setminus I|}\,\,\omega^{(alt)}_{I}(X\,|\,Y_{J}).$

By the above Remark $\omega^{(alt)}_{I}(X\,|\,Y_{J})=0$ for proper subsets
$I\subset J$, since $\omega_{J}(X)=\prod\nolimits_{j\in
I}\omega(x_{j})\prod\nolimits_{j\in J\setminus I}1$. Therefore,
$\lambda^{(alt)}_{J}(X\,|\,Y_{J})=\omega^{(alt)}_{J}(X\,|\,Y_{J})$, and
Theorem follows. ∎

We consider periodic point distributions ${\mathscr{D}}$ on
${\mathbb{R}}^{d},{\mathscr{D}}+M={\mathscr{D}},M\in{\mathbb{Z}}^{d},$ with a
finite set of residues ${\mathscr{D}}/{\mathbb{Z}}$. Notice that instead of
point distributions, arbitrary periodic complex Borel measures on
${\mathbb{R}}^{d}$ finite on ${\mathbb{T}}^{d}$ could be considered but we do
not consider such a generalization in order not to complicate the notation.

We define the local discrepancy

$\displaystyle
L[{\mathscr{D}},Y]=\sum\nolimits_{X\in{\mathscr{D}}/{\mathbb{Z}}^{d}}L(X,Y),$
(12)

and the $L_{q}$– discrepancies

$\displaystyle
L_{q}[{\mathscr{D}}]=\left(\int\nolimits_{{\mathbb{K}}^{d}}|L[{\mathscr{D}},Y]|^{q}\,\mathrm{d}Y\right)^{1/q},\,0<q<\infty,\quad\quad
L_{\infty}[{\mathscr{D}}]=\sup\nolimits_{Y\in{\mathbb{K}}^{d}}|L[{\mathscr{D}},Y]|.$

We also introduce the shifted discrepancies

$\displaystyle
L_{q}^{*}[{\mathscr{D}}]=\sup\nolimits_{Z\in{\mathbb{T}}^{d}}L_{q}[{\mathscr{D}}+Z],\quad\quad
L_{\infty}^{*}[{\mathscr{D}}]=\sup\nolimits_{Z\in{\mathbb{T}}^{d}}L_{\infty}[{\mathscr{D}}+Z],$

and their mean values

$\displaystyle\lambda_{J}[{\mathscr{D}}]=\sum\nolimits_{X\in{\mathscr{D}}/{\mathbb{Z}}^{d}}\lambda_{J}(X),\quad\quad\lambda^{*}_{J}[{\mathscr{D}}]=\sup\nolimits_{Z\in{\mathbb{K}}^{d}}\lambda_{J}[{\mathscr{D}}+Z].$

The above Theorem immediately implies the following.

Lemma 1. The local discrepancy $L[{\mathscr{D}},Y]$ satisfies the identity

$\displaystyle
L[{\mathscr{D}},Y]=\sum\nolimits_{J\subseteq[d]}\,v_{J^{\prime}}(Y)\,\lambda^{(alt)}_{J}\,[{\mathscr{D}}\,|\,Y_{J}],$
(13)

where
$\lambda^{(alt)}_{J}[{\mathscr{D}}\,|\,Y_{J}]=\sum\nolimits_{X\in{\mathscr{D}}/{\mathbb{Z}}^{d}}\,\lambda^{(alt)}_{J}(X\,|\,Y_{J})$.

The discrepancy $L_{\infty}[{\mathscr{D}}]$ satisfies the inequality

$\displaystyle
L_{\infty}[{\mathscr{D}}]\leq\sum\nolimits_{J\subseteq[d]}\,2^{|J|}\,\lambda^{*}_{J}[{\mathscr{D}}].$
(14)

###### Proof.

The identity (13) follows from (8). The definition (6) implies
$\lambda^{(alt)}_{J}[{\mathscr{D}}\,|\,Y_{J}]\leq
2^{|J|}\,\lambda^{*}_{J}[{\mathscr{D}}],$ and the inequality (14) follows,
since $0\leq v_{J^{\prime}}(Y)\leq 1$. ∎

The $L_{q}^{*}[{\mathscr{D}}]$-discrepancies can be easily estimated by means
from below.

Lemma 2. For $1\leq q\leq\infty$ and any subset $J\subseteq[d]$, the
discrepancy $L_{q}^{*}[{\mathscr{D}}]$ satisfies the inequality

$\displaystyle L_{q}^{*}[{\mathscr{D}}]\geq
2^{d-|J|}\,\lambda^{*}_{J}[{\mathscr{D}}].$ (15)

###### Proof.

We have

$\displaystyle L_{q}^{*}[{\mathscr{D}}]\geq L_{1}^{*}[{\mathscr{D}}]$
$\displaystyle=\sup\nolimits_{Z\in{\mathbb{T}}^{d}}\int\nolimits_{{\mathbb{K}}^{d}}|L[{\mathscr{D}}+Z,Y]|\,\mathrm{d}Y\geq\sup\nolimits_{Z\in{\mathbb{T}}^{d}}\left|\int\nolimits_{{\mathbb{K}}^{d}}L[{\mathscr{D}}+Z,Y]\,\mathrm{d}Y\right|$
$\displaystyle\geq\sup\nolimits_{Z\in{\mathbb{T}}^{d}}\left|\lambda[{\mathscr{D}}+Z]\right|=\lambda^{*}[{\mathscr{D}}].$

For simplicity, we put
$Z=(Z_{J},Z_{J^{\prime}})\in{\mathbb{T}}^{d},Z_{J}\in{\mathbb{T}}^{|J|},Z_{J^{\prime}}\in{\mathbb{T}}^{|J^{\prime}|}$,
and continue

$\displaystyle\lambda^{*}[{\mathscr{D}}]$
$\displaystyle=\sup\nolimits_{Z_{J}}\sup\nolimits_{Z_{J^{\prime}}}|\lambda[{\mathscr{D}}+Z]|\geq\sup\nolimits_{Z_{J}}\int\nolimits_{{\mathbb{K}}^{|J^{\prime}|}}|\lambda[{\mathscr{D}}+Z]|\,\mathrm{d}Z_{J^{\prime}}$
$\displaystyle\geq\sup\nolimits_{Z_{J}}\left|\int\nolimits_{{\mathbb{K}}^{|J^{\prime}|}}\lambda[{\mathscr{D}}+Z]\,\mathrm{d}Z_{J^{\prime}}\right|=2^{-|J^{\prime}|}\sup\nolimits_{Z_{J}}|\lambda[{\mathscr{D}}+Z_{J}]|=2^{d-|J|}\,\lambda^{*}_{J}[{\mathscr{D}}],$

that completes the proof. ∎

The next simple fact is well–known, see, for example, [2, p. 4]. For
completeness, we give a short proof.

Lemma 3. The $L_{\infty}$– and $L_{\infty}^{*}$– discrepancies are equivalent:

$\displaystyle
L_{\infty}[{\mathscr{D}}]\,\leq\,L_{\infty}^{*}[{\mathscr{D}}]\,\leq\,3^{d}L_{\infty}[{\mathscr{D}}].$
(16)

###### Proof.

The left inequality (16) is obvious. Let us prove the right. For
$y\in[0,1],\,z\in[0,1)$, we introduce the notation

$\displaystyle\delta_{y,z}=\begin{cases}1,&\text{if \,$y+z\geq 1$,}\\\
0,&\text{otherwise,}\end{cases}$

The shifted indicator function, see (2) and (10), can be written in the form

$\displaystyle\chi(x+z,y)=\chi(x,y+{z})[1-\delta_{y,z}]+\chi(x,y+{z}-1)\,\delta_{y,z}-\chi(x,{z})+\chi(x,1)\,\delta_{y,z},$

and similarly
$v(y)=v(y+z)[1-\delta_{y,z}]+v(y+z-1)\,\delta_{y,z}-v(z)+v(1)\,\delta_{y,z}.$
Moreover, each of these formulas contains at most three non-zero terms.
Substituting these formulas into the definitions (1) and (12), we find that
the shifted discrepancy can be written as the sum
$L[{\mathscr{D}}+Z,Y]=\sum\nolimits_{k}c_{k}\,L[{\mathscr{D}},V_{k}],$ with
some vectors $V_{k}=V_{k}(Y,Z)\in{\mathbb{K}}^{d}$ and coefficients
$c_{k}=c_{k}(Y,Z)$ equal either $\pm 1$ or $0$. Moreover, the sum contain at
most $3^{d}$ non-zero terms. This implies the right inequality (16). ∎

Lev [2] established the equivalence of $L_{\infty}$– and $L_{q}^{*}$–
discrepancies: $L_{\infty}[{\mathscr{D}}]\,\cong\,L_{q}^{*}[{\mathscr{D}}]$,
for all $q\geq 1$ with the implicit constants depending only on the dimension
$d$. Another proof of the equivalence was given later by Kolountzakis, we
refer to the survey article [1] by Chen for a detailed discussion of these
issues.

The equivalence of the $L_{\infty}$– and $L_{q}^{*}$– discrepancies can be
easily derived from the foregoing statements. We will formulate and prove the
corresponding result in the following somewhat more general form.

Corollary (Lev’s Equivalence). For $0<q<\infty$, the $L_{\infty}$– and
$L_{q}^{*}$– discrepancies are equivalent:

$\displaystyle
3^{-d}\,L^{*}_{q}[{\mathscr{D}}]\,\leq\,\,L_{\infty}[{\mathscr{D}}]\,\leq\,C_{d,q}\,L^{*}_{q}[{\mathscr{D}}],$
(17)

where the constant

$\displaystyle C_{d,q}=\begin{cases}\,(5/2)^{d},&\text{if $1\leq
q<\infty$,}\\\ \,(5/2)^{d/q}\,3^{d/q\,-d},&\text{if $0<q<1.$}\end{cases}$ (18)

We did not seek to obtain the best constant $C_{d,q}$ in (17), however, we
note that the formula (18) correctly reflects the order of the constant for
large and small $q$.

###### Proof of Corollary..

The proof consists of three steps.

(i) For $0<q<{\infty}$, the lower bound (17) follows from Lemma 3:
$L^{*}_{q}[{\mathscr{D}}]\,\leq
L_{\infty}^{*}[{\mathscr{D}}]\,\leq\,3^{d}L_{\infty}[{\mathscr{D}}],$ since
$L^{*}_{q}[{\mathscr{D}}]$ is a non-decreasing function of $q>0$.

(ii) For $1\leq q<\infty$, the upper bound (17) follows from Lemma 1 and Lemma
2. Substituting (15) into (14), we obtain

$\displaystyle L_{\infty}[{\mathscr{D}}]\leq
2^{-d}\,\left(\sum\nolimits_{J\subseteq[d]}\,\,2^{2|J|}\right)\,\,L_{q}^{*}[{\mathscr{D}}]=(5/2)^{d}\,\,L_{q}^{*}[{\mathscr{D}}].$

(iii) Finally, for $0<q<1$, the upper bound (17) follows from the logarithmic
convexity of $L_{q}^{*}[{\mathscr{D}}]$ as a function of $q>0$. By the
standard interpolation at points $q<1<p$, we obtain

$\displaystyle
L_{1}^{*}[{\mathscr{D}}]\leq(L_{q}^{*}[{\mathscr{D}}])^{q\frac{p-1}{p-q}}\,(L_{p}^{*}[{\mathscr{D}}])^{p\frac{1-q}{p-q}}.$

This inequality takes the form

$\displaystyle
L_{1}^{*}[{\mathscr{D}}]\leq(L_{q}^{*}[{\mathscr{D}}])^{q}\,(L_{\infty}^{*}[{\mathscr{D}}])^{1-q},$
(19)

as $p\to\infty$.

The bounds $(5/2)^{-d}\,L_{\infty}[{\mathscr{D}}]\leq
L_{1}^{*}[{\mathscr{D}}]$ and $L^{*}_{\infty}[{\mathscr{D}}]\leq
3^{d}\,L_{\infty}[{\mathscr{D}}],$ are already established. Substituting these
bounds into (19), we obtain the upper bound (17) with the constant (18).

The proof of Corollary is completed. ∎

## References

  * [1] W. W. L. Chen, Results and problems old and new in discrepancy theory, _Discrepancy Theory_ , pp. 21–42 (Radon Series on Computational and Applied Mathematics), de Gruyter, 2020.
  * [2] V. F. Lev, Translations of nets and relationship between supremum and $L^{\kappa}$-discrepancies, _Acta Math. Hungar._ , vol. 70 (1-2), pp. 1–12, 1996.


Free boundary regularity for the inhomogeneous one-phase Stefan problem
]Fausto Ferrari, Nicolò Forcillo, Davide Giovagnoli, David Jesus

Fausto Ferrari: Dipartimento di Matematica, Università di Bologna, Piazza di Porta
S.Donato 5, 40126, Bologna-Italy
Nicolo Forcillò: Departments of Mathematics, Michigan state university, 619 Red Cedar Road
East Lansing, MI 48824, USA
Davide Giovagnoli: Dipartimento di Matematica, Università di Bologna, Piazza di Porta
S.Donato 5, 40126, Bologna-Italy
David Jesus: Dipartimento di Matematica, Università di Bologna, Piazza di Porta
S.Donato 5, 40126, Bologna-Italy
This research is partially supported by PRIN 2022 7HX33Z - CUP J53D23003610006. F.F. is partially supported by  INDAM-GNAMPA project 2024: "Free boundary
problems in noncommutative structures and degenerate operators.
Keywords: Inhomogeneous Stefan problem, free boundary problems, flat free boundaries, parabolic operators, viscosity solutions.
MSC: 35R35, 35K55, 80A22
Abstract. In this paper, we prove that flat free boundaries of  solutions to inhomogeneous one-phase Stefan problem are $C^{1,\alpha}$.

§ INTRODUCTION

The goal of this paper is to introduce new results concerning the regularity of flat free boundaries  in a one-phase Stefan problem governed by an inhomogeneous equation. More precisely, we deal with nonnegative viscosity solutions $u: \Omega \times [0,T] \mapsto \mathbb{R}$ of the following  free boundary problem

\begin{equation}\label{eq:StefanNH}
\begin{cases}
\partial_{t} u - \Delta u  =  f& \text{in} \,\, \Omega^+(u):=\left(\Omega \times [0,T] \right) \cap \{ u > 0  \},\\
\partial_{t} u = \nor{\nabla u}^2& \text{on} \,\, F(u):=\left(\Omega \times [0,T] \right) \cap \partial\{ u > 0  \},
\end{cases}
\end{equation}

where $\Omega \subset \mathbb{R}^{n}$, $T > 0$,  $f: \Omega \times [0,T] \mapsto \mathbb{R}, \, f \in C(\Omega \times [0,T]) \cap L^{\infty }(\Omega \times [0,T])$, and $f\geq 0$.

The Stefan problem consistently attracts attention due to its modelling of fundamental physical phenomena, especially the processes of melting or freezing ice. This topic assumes a pivotal role within the context of free boundary problems.

Concerning this research field, we have to cite some pioneering papers that represent a reference point for people working on this subject. In particular,  the trilogy <cit.> collects the last fundamental contribution to the free boundary regularity in the homogeneous Stefan problem. Further achievements can be found in <cit.> regarding the  existence of analytic solutions, and in <cit.> about the relationship between the smoothness of solutions and the regularity of  initial  data, moreover, <cit.> focused on the existence of a solution, and <cit.> adopts the same techniques as in <cit.> to study operators with variable coefficients. In these papers, the authors obtained sharp results about Stefan problem, in one or two-phase settings, without the source term.

Recently, in <cit.>, a novel methodology for investigating the regularity of free boundaries has emerged. The authors address the regularity in a  parabolic homogeneous one-phase scenario employing techniques such as the hodograph transform and viscosity theory, and in particular, by building upon the concept of a flat free boundary within an evolving framework. This approach follows the spirit of the  elliptic one introduced in <cit.>, and then widely applied in many other different cases, see for instance <cit.> for a survey and <cit.> for some further recent results obtained in the elliptic degenerate nonlinear setting.
In the same direction, we  also point out the new contribution in <cit.>, where a parabolic homogeneous  two-phase problem with a free boundary condition of transmission-type was considered. There, the boundary condition only involves the spatial derivatives, differently from the one used in <cit.>.

In our paper, we deal with the inhomogeneous one-phase case (<ref>), adding a new contribution that, to our knowledge, hasn't appeared in the literature yet. We follow the approach suggested in <cit.>, instead of the one applied in <cit.>, due to the apparently wider flexibility of the tools involved.

The condition on the free boundary, given by $ \partial_{t} u = \nor{\nabla u}^2$   on  $F(u)$, prescribes the law of conservation of energy on the free boundary, which is moving from the positive phase to the zero level set.
The physical interpretation of the problem models the melting of ice. Indeed, we are in the simpler case encoded by the general law of conservation of energy, summarized by

\begin{equation}\label{energy_conservation_law}
-\frac{ \partial_{t} u^-}{\nor{\nabla u^-}} =\frac{ \partial_{t} u^+}{\nor{\nabla u^+}} = \nor{\nabla u^+}-\nor{\nabla u^-}, \quad\textit{on}\:\: F(u),
\end{equation}

see <cit.> for a detailed discussion as well as <cit.>, that characterizes the more general Stefan problem.
Here, $u\in C(\Omega\times [0,T])$, $u^+=\sup\{u,0\}$, $u^-=\sup\{-u,0\}$, and $u$ satisfies in a viscosity sense:

\begin{equation}\label{two_phase_equations}
\begin{cases}
\partial_tu - \Delta u=f,& \quad \textit{ in }\{(x,t)\in \Omega \times [0,T]:\:\:u>0\},\\ \partial_tu - \Delta u =f,& \quad \textit{ in } \{(x,t)\in \Omega\times [0,T]:\:\:u\leq 0\}^0.
\end{cases}
\end{equation}

(0,0) circle (1.73cm);

at (-0,1.3) $u  > 0$;
at (0,-0.9) $u  < 0$;
at (0,0.2) $u  = 0$;
[->,blue] (0,-0.2) – (0,0.05);
[->,red] (0,0.515) – (0,0.265);
[scale=1, domain=-1.481:1.481, smooth, variable=, black] plot (, 0.5**-0.2);
[scale=1, domain=-0.854:0.87, smooth, variable=, black] plot (, cos(r  -0.020252976382 r)-(0.75-0.2732240183828)));
[->,red] (-1.2,0.515) – (-1.5,0.265);
[->,red] (1.2,0.515) – (1.7,0);
at (1.35,0.15) $\nu_3$;
[->,blue] (-1.2,0.515) – (-0.5,1.1);
at (-1,0.9) $\nu_1$;
[->,blue] (1.2,0.515) – (0.95,0.765);
at (0.2,0.38) $\nu_2$;

Description of the free boundary movement

Figure <ref> describes a possible configuration of the  resulting contributions from the law of conservation of energy (<ref>) from each phase in the free boundary movement. For example, in $\nu_1$, since the contribution from the negative phase is greater than the one from the positive phase, that is, $\nor{\nabla u^-} >\nor{\nabla u^+}$, the free boundary will move in the direction of the positive phase. In $\nu_3$, instead, the opposite happens and thus the free boundary will move towards the negative phase. In $\nu_2$, the only contribution comes from the positive phase, hence the free boundary behaves like in the one-phase problem.

Since we are interested in the one-phase Stefan problem, we assume $u\geq 0$ and thus we recover equation (<ref>). We emphasize that the condition on the free boundary, $\partial_t u=|\nabla u|^2$, only describes the melting of ice, since the time derivative is positive, and so the the interface between the two states is pushed over time towards the ice side. Therefore, it is clear that we need to assume the sign condition $f\geq 0$, so that we are adding heat to the system.
Otherwise, it would be possible for the free boundary to move in a non-smooth way, or it could even appear suddenly in a different part of the domain, and it seems very difficult to obtain regularity results for the free boundary in these cases. Indeed, this informal discussion translates into a very precise technical difficulty that appears when $f$ is allowed to be negative, which is the fact that the Hopf Lemma, as it is stated in the homogeneous setting, is false in this case.

We recall that, while facing problem (<ref>) as well as (<ref>) - (<ref>), due to the different scalings in the PDE and the boundary condition, often emerges a dichotomy already remarked in <cit.> associated with the way of using the rescale argument.
In fact, there exist at least two basic ways of rescaling this problem: the hyperbolic and the parabolic one.

In the hyperbolic rescaling, for instance centered at $(0,0)\in F(u)$,  and given by

\begin{equation*}
(x,t) \mapsto u_{\lambda}^{H} := \frac{u(\lambda x, \lambda t )}{\lambda}, \quad (x,t) \in B_1 \times [-K^{-1},0],
\end{equation*}

we obtain that  $u_{\lambda}^{H}$ solves a Stefan problem having a rescaled invariant boundary condition, but with a possibly large diffusion coefficient $\lambda^{-1}$ in the equation, that is

\begin{equation}\label{eq:StefanHyper}
\begin{cases}
\lambda \partial_{t} u_{\lambda}^{H}-\Delta u_{\lambda}^{H}  = \lambda f(\lambda x, \lambda t)& \text{in} \,\, \left(B_1 \times [-K^{-1},0]  \right) \cap \{ u_{\lambda}^{H} > 0  \},\\
\partial_{t} u_{\lambda}^{H} = \nor{\nabla u_{\lambda}^{H}}^2& \text{on} \,\, \left(B_1 \times [-K^{-1},0]\right) \cap \partial\{u_{\lambda}^{H} > 0  \}.
\end{cases}
\end{equation}

On the other hand, in the parabolic rescaling,

\begin{equation*}
(x,t) \mapsto u_{\lambda}^{P} := \frac{u(\lambda x, \lambda^2 t )}{\lambda}, \quad (x,t) \in B_1 \times [-(K\lambda)^{-1},0]
\end{equation*}

we get an equation with the same form, but having a different boundary condition depending on the coefficient $\lambda$. In fact,

\begin{equation}\label{eq:StefanPara}
\begin{cases}
\partial_{t} u_{\lambda}^{P} - \Delta u_{\lambda}^{P}  = \lambda f(\lambda x, \lambda^2 t)& \text{in} \,\, \left(B_1 \times [-(K\lambda)^{-1},0]  \right) \cap \{ u_{\lambda}^{P} > 0  \},\\
\partial_{t} u_{\lambda}^{P} = \lambda \nor{\nabla u_{\lambda}^{P}}^2& \text{on} \,\, \left(B_1 \times [-(K\lambda)^{-1},0]\right) \cap \partial\{u_{\lambda}^{P} > 0  \}.
\end{cases}
\end{equation}

These two alternative representations of the problem highlight the main difference with respect to the elliptic situation, where only the hyperbolic rescaling appears.
Another crucial disparity lies in how the flatness of the free boundary is defined, as in the Stefan problem, the free boundary evolves over time.

Our main result is the following one, which establishes $C^{1,\alpha}$ regularity of flat free boundaries, that is, at each time, the free boundary is trapped in a small strip, see Definition <ref>.

Fix a large constant $K$  and let $u$ be a solution to the one-phase Stefan problem (<ref>) in $B_{\lambda} \times [-K^{-1}\lambda,0]$ for some $\lambda \leq 1$. Assume that

\begin{equation*}
\nor{u} \leq K \lambda, \qquad u(x_0,t) \geq K^{-1} \lambda \; \text{ for some }x_0 \in B_{\frac{3}{4}\lambda}.
\end{equation*}

There exists $\varepsilon_0$ depending only on $K$ and $n$ such that if, for each $t$, $\p_x \{ u(\cdot, t) >0 \}$ is $\varepsilon_0-$flat in $B_{\lambda}$, then the free boundary $\p\{ u >0 \}$ (and $u$ up to the free boundary) is $C^{1,\alpha}$ in $B_{\frac{\lambda}{2}} \times [-(2K)^{-1} \lambda, 0]$.

In order to prove Theorem <ref> we need  the following one that transfers conveniently the flatness assumption on the whole function $u$.  It is essential to highlight that the primary distinctions from <cit.> manifest here, mainly due to the inclusion of the source term $f$, which we assume, without loss of generality, to be very small, as elaborated upon later.

Fix a large constant $K$ and let $u$ be a solution to the one-phase Stefan problem (<ref>) in $B_{2\lambda} \times [-2\lambda,0]$ for some $\lambda \le
Assume that $0\in \p \{u>0\},$and

$$ a_n(t) \, \left (x_n - b(t)-\varepsilon_1 \lambda \right)^+ \, \, \le u \, \, \le a_n(t) \, \left(x_n - b(t)+\varepsilon_1 \lambda \right)^+,$$

$$K^{-1}\le a_n \le K, \quad \quad |a_n'(t)| \le \lambda^{-2}, \quad \quad b'(t)=-a_n(t),\quad\quad  \|f\|_{L^\infty}\leq \varepsilon_1^2,$$

for some small $\varepsilon_1$ depending only on $K$ and $n$. Then in $B_\lambda \times[-\lambda ,0]$ the free boundary $\p \{u>0\}$ is a $C^{1,\alpha}$ graph in
the $x_n$ direction.

In Subsection <ref> we include the definitions of flatness both for the free boundary and the solution and a discussion relating these two notions in the inhomogenenous case which implies the equivalence of Theorem <ref> and Theorem <ref>, here the sign condition $f \geq 0$ play a key role in the use of Hopf lemma.

As in <cit.>, our argument relies on applying the hodograph transform to solutions of the Stefan problem (<ref>), thus reducing the problem to studying multi-valued viscosity solutions to the following equation

\begin{equation}\label{eq:QuasiStefanNH}
\begin{cases}
\partial_{t} u - \tr( A(\nabla u) D^2 u) + (\partial_{x_n}u) \,  f(x',u,t)=0  &\text{in} \,\, \{ x_n > 0 \},\\
\partial_{t} u = g(\nabla  u) &\text{on} \,\,  \{ x_n = 0 \},
\end{cases}
\end{equation}

where $g$ encodes the application of the hodograph transform. Its precise definition can be found later on in (<ref>) and (<ref>). The problem (<ref>) is new due to the presence of the function $f$ which introduces a new gradient term.

The hodograph transform straightens the free boundary and transfers the difficulty of the problem to the fact that solutions may become multi-valued near the new fixed boundary. Thus the problem of studying regularity of the free boundary becomes equivalent to studying regularity up to the fixed boundary of solutions to (<ref>).  As a consequence, the proof of Theorem <ref> relies on an improvement of flatness result of solutions to (<ref>).

Note that the presence of the source term $f$ introduces, after the transformation, a new term depending on the derivative that changes considerably the equation. In particular, this introduces a challenge to recover the Harnack inequality away from the boundary. In the homogeneous case this follows simply from parabolic interior Harnack inequality, however the presence of the gradient term requires a more sophisticated Harnack for the following extremal equations

\begin{align*}
\partial_t u-\mathcal{M}^\pm(D^2u)\mp\Lambda|\nabla u|+h=0
\end{align*}

which can be found in <cit.>. This is the tool used obtain the  inhomogeneous H-property, see Subsection <ref>.

Equation (<ref>) does not fit into the framework of the homogeneous one and thus our main contribution is the study of solutions of this new equation. Nevertheless, we prove that the limiting profile we obtain is the same observed in the homogeneous case as the linearised equation

\begin{equation*}
\begin{cases}
\lambda \partial_t v = \tr(A_{\lambda}(t) D^2 v )&\text{in} \,\, \{ x_n > 0 \},\\
\partial_t v = \gamma_{\lambda} (t) \cdot \nabla v &\text{on} \,\,  \{ x_n = 0 \},
\end{cases}
\end{equation*}

\begin{equation*}
A_{\lambda}(t):= A(a(\lambda t)), \quad \gamma_{\lambda} (t) :=  \nabla g (a(\lambda t)).
\end{equation*}

Hence, we can exploit the Schauder-type estimates obtained in <cit.> with respect to an appropriate distance $d_\lambda$ which encodes features of both the parabolic and hyperbolic scalings of the problem. The improved regularity for the linearized problem is then transferred to (<ref>) by a comparison principle.

Now we briefly describe the structure of the paper. In Section <ref> we provide the necessary notation and definitions, and we include a discussion on the relation between the flatness of the free boundary and the solution, and how this gives the equivalence between Theorem <ref> and Theorem <ref>.

In Section <ref>, we apply the hodograph transform to change equation (<ref>) into (<ref>) and we discuss how the initial assumptions translate into this new context. Then the crucial inhomogeneous $H$-property is established. We finish this section by stating the main improvement of flatness result and by introducing the limiting profile. In this regard, we devote Section <ref>, to state regularity results for the limiting profile and obtain comparison principles which relate this profile to our equation.

In Section <ref>, we study properties of the error function, which is defined as $u-\ell$ where $\ell$ is the affine approximation of $u$, for each time $t$. This section contains most of the new and important results of the paper, in particular oscillation decay, Harnack inequality up to the boundary and almost Hölder estimates for the error. This allows in Section <ref> to give the proof of the improvement of flatness and of the main Theorem <ref>.

§ PRELIMINARIES

§.§ Notation

Given $(x,t)\in \R^{n+1}$, we call $x\in \mathbb{R}^n$ the space variable and $t \in \R$  the time variable. A point $x\in \mathbb{R}^n$ will sometimes be written as $x=(x',x_n)$, where $x'\in \R^{n-1}$ and $x_n \in \R$.
Given a function $u(x,t)$, we denote by $\nabla u$ the gradient of $u$ with respect to $x$, $D^2 u$ the Hessian with respect to $x$, $\nabla_{x'} u$ the gradient with respect to $x'$, $\partial_{x_n}u $  and $\partial_{t} u$ the partial derivatives of $u$ with respect to $x_n$ and $t$, respectively.

We call a parabolic cylinder any set of the form

$$\mathcal{C}=Q \times (t_1,t_2],$$

where $U$ is a smooth bounded domain in $\R^n$ and $t_1 < t_2$. We denote by $\partial_p \mathcal{C}$ the parabolic boundary of $\mathcal{C}$, i.e., $\partial_p \mathcal{C}=\left( Q\times \{t_1\}  \right)  \cup   \left( \partial Q \times [t_1,t_2] \right) $.

Given $(x,t), (y,s)\in \R^{n+1}$, we define the parabolic distance as

d_p((x,t), (y,s)) :=\big(|x-y|^2 + |t-s|\big)^{1/2}.

Let $C(\mathcal{C})$ be the set of continuous function in $\mathcal{C}$. For $\alpha \in (0,1]$ and $u\in C(\mathcal{C})$, define the parabolic Hölder semi-norm as

\begin{equation*}
[u]_{C^{\alpha}(\mathcal{C})} := \sup_{\substack{(x,t)\neq (y,s)}} \frac{\nor{u(x,t)-u(y,s)}}{d_p((x,t),(y,s))^{\alpha}}.
\end{equation*}

The parabolic Hölder norm is defined as

\begin{equation*}
\norm{u}_{C^{\alpha}(\mathcal{C})} := \norm{u}_{L^{\infty}(\mathcal{C})} +[u]_{C^{\alpha}(\mathcal{C})}.
\end{equation*}

Further we define the Hölder semi-norm in $t$ to be

\begin{equation*}
[u]_{C_t^{\alpha}(\mathcal{C})} := \sup_{ \substack{(x,t)\neq (x,s)}} \frac{|u(x,t)-u(x,s)|}{|t-s|^\alpha}.
\end{equation*}

This is necessary to introduce  the $C^{1,\alpha}(\mathcal{C})$-norm as

\begin{equation*}
\norm{u}_{C^{1,\alpha}(\mathcal{C})}:=\norm{u}_{L^{\infty}(\mathcal{C})} + \norm{\nabla u}_{L^{\infty}(\mathcal{C})} +[\nabla u]_{C^{\alpha}(\mathcal{C})} + [u]_{C_t^{\alpha/2}(\mathcal{C})}.
\end{equation*}

We denote by $C^{\alpha}(\mathcal{C})$ (resp. $C^{1,\alpha}(\mathcal{C})$) the space of continuous (resp. differentiable in space) functions on $\mathcal{C}$ which are bounded in the $\norm{\cdot}_{C^{\alpha}(\mathcal{C})} $ (resp. $\norm{\cdot}_{C^{1,\alpha}(\mathcal{C})}$) norm.

From this point forward, for $n \geq 2$, given $r > 0$ we set

\begin{equation*}
\begin{aligned}
&Q_{r} := (-r,r)^{n}, \;  Q_{r}^{+}:=  Q_{r} \cap \{ x_n \geq 0 \},  \;  Q_{r}(x_0):=  x_0 + Q_{r} , \\
&\mathcal{C}_r := (Q_r \cap \{x_n >0 \}) \times (-r,0], \\ &\mathcal{F}_r:=\{(x,t) | \,  x \in Q_r \cap \{x_n = 0 \}, \, t \in (-r,0] \}.
\end{aligned}
\end{equation*}

It will be useful to define the Dirichlet boundary of $\mathcal{C}_1 $ as

\begin{equation*}
\p_D \mathcal{C}_1:= \partial \mathcal{C}_1 \cap \left( \{ t=-1\} \cup \{ x_n=1\} \cup_{i=1}^{n-1} \{ \nor{x_i}=1\} \right),
\end{equation*}

which together with $\mathcal{F}_1$ represent the parabolic boundary of $\mathcal{C}_1$.

Recall the definition of extremal Pucci operators with ellipticity constants $K^{-1}$ and $K$

\begin{equation}\label{Pucci} \mathcal M_K^+ (N) = \max_{K^{-1} I \le A \le K I} \quad tr \, \,  A N, \quad \quad  \quad \mathcal M_K^- (N) = \min_{K^{-1} I \le A \le K I} \quad tr \, \, A N.\end{equation}

§.§ Viscosity solution

We define what we mean by a viscosity solution to (<ref>).

We say that $u\in C(\Omega\times [0,T])$ is a viscosity solution to (<ref>) if its graph cannot be touched by above (resp. below), at a point $(x_0,t_0)$ in a parabolic cylinder $B_r (x_0) \times (t_0 - r^2,t_0],$
by the graph of a classical strict supersolution $\varphi^{+}$ (resp. subsolution).
By a classical strict supersolution we mean that $\varphi(x,t) \in C^{2}$, $\nabla \varphi \neq 0$ and it solves

\begin{equation}\label{eq:StefanNH1}
\begin{cases}
\partial_{t} \varphi -\Delta \varphi > f& \text{in} \,\, \left(\Omega \times [0,T] \right) \cap \{ \varphi > 0  \},\\
\partial_{t} \varphi > \nor{\nabla \varphi}^2& \text{on} \,\, \left(\Omega \times [0,T] \right) \cap \partial\{ \varphi > 0  \}.
\end{cases}
\end{equation}

Similarly we can define a strict classical subsolution.

§.§ Flatness and smallness regime

To relate our inhomogeneous equation to the homogeneous one studied in <cit.> we start by performing a rescaling to ensure that we can assume, without loss of generality, that the source term $f$ is arbitrarily small. Indeed, if $u$ is the solution to (<ref>) in $B_1\times [-1,0]$, and calling

\[
\Lambda := \frac{\|f\|^\frac{1}{2}_{L^\infty}}{\varepsilon_1},
\]

then the function

\[
\bar u(x,t)=u(\Lambda^{-1}x, \Lambda^{-2} t), \quad (x,t)\in B_{\Lambda}\times [-\Lambda^2, 0]
\]

is a solution to

\begin{align*}
\begin{cases}
\partial_{t} \bar u - \Delta \bar u  =  \bar f& \text{in} \,\, \left(B_{\Lambda}\times [-\Lambda^2, 0] \right) \cap \{ \bar u > 0  \},\\
\partial_{t} \bar u = \nor{\nabla \bar u}^2& \text{on} \,\, \left(B_{\Lambda}\times [-\Lambda^2, 0] \right) \cap \partial\{ \bar u > 0  \},
\end{cases}
\end{align*}

with $\|\bar f\|_{L^\infty}\leq \varepsilon_1^2$. For ease of notation, we will refer to this function as $u$.

We now comment on the flatness assumption of the free boundary and how it relates to the flatness assumption on $u$. We start with the following definition.

Let $\partial_x\{u(\cdot,t)>0\}$ denote the spatial boundary in $\R^n$ of $\{u(\cdot,t)>0\}$. We say that $\partial_x\{u>0\}$ is $\varepsilon_0$-flat in $B_\lambda(x_0)$ if, for every time $t$, there exists a direction $\nu$ such that

\[
\partial_x\{u>0\}\cap B_\lambda(x_0)\subset \{|\ps{x-x_0}{\nu}|\leq \varepsilon_0\lambda\}
\]

\begin{align*}
u=0 \quad &\textit{ in } \quad \{|\ps{x-x_0}{\nu}|\leq -\varepsilon_0\lambda\},\\
u>0 \quad &\textit{ in } \quad \{|\ps{x-x_0}{\nu}|\geq \varepsilon_0\lambda\}.
\end{align*}

We can also repeat the argument in <cit.> to conclude that solutions with flat free boundaries are flat.

Assume that $u$ solves (<ref>),

\[
|u|\leq K, \quad u(x_0,t)\geq K^{-1} \quad \textit{for some } \quad x_0\in B_{3/4},
\]

\[
0\in \partial_x\{ u(\cdot,0)>0\}, \quad \textit{and}\quad \partial_x \{u(\cdot,t)>0\} \textit{ is $\varepsilon_0$-flat in } B_1.
\]

Then for all $\eta >0$ we have

\[
a_n(t)(x_n-b(t)-\eta^{1+\beta})^+\leq u \leq a_n(t)(x_n-b(t)+\eta^{1+\beta})^+
\quad \mbox{in }B_{\eta} \times [-\lambda^{-1}\eta,0],\]

with $\beta=1/20$ and for $c,C>0$ universal,

$$c \le a_n(t) \le C, \quad |a_n'(t)| \le \eta^{\beta-2}, \quad \quad b'(t)=- \lambda a_n(t), \quad b(0)=0,$$

provided that $\varepsilon_0 \le c(\eta,K)$.

The proof of this result follows the lines of <cit.>. The presence of a source term plays a role in the regularity up to the boundary <cit.>
and enters in the Hopf lemma. We point out that the sign condition on $f$ is crucial and allows the use of Hopf lemma.

Lemma <ref> shows that the graph of $u$ is $\eta^\beta$-flat in the domain $B_{\eta\lambda} \times [-\eta \lambda,0]$ provided that $\varepsilon_0 \leq c(\eta,K).$ Thus, selecting $\eta$ sufficiently small, the hypotheses of Theorem <ref> are fulfilled.

There are some special cases in which the behavior of the free boundary for the inhomogeneous problem turns out to be even simpler than the case without a given source. This situation occurs, for instance, if we consider $f=  a^2$. Here, instead of considering a perturbation in the form of a travelling wave solution

\begin{equation*}
(e^{a x_n +a^2 t} -1)^{+},
\end{equation*}

we can take into account directly

\begin{equation*}
(a x_n + a^2 t)^{+},
\end{equation*}

that,  for each $t$ fixed, exactly describes an hyperplane. Note that this is yet another instance where considering a negative source term significantly complicates the problem, since there is no easy solution if we consider $f=-1$.

§ THE HODOGRAPH TRANSFORM IN THE INHOMOGENEOUS CASE

In order to adapt the argument of the homogeneous case introduced in <cit.>, we need to study how the hodograph transform modifies (<ref>)  with the addition of the source term.
To avoid confusion in this section we will use $u_{S}$ to denote a solution of the Stefan problem (<ref>) and $u_H$ its hodograph transform.
Precisely, we look the graph of $u_S$ in $\mathbb{R}^{n+2}$

\begin{equation*}
\Gamma := \{ (x,x_{n+1},t) | \; x_{n+1}= u_S(x_1, x_2, \dots, x_n,t)  \}
\end{equation*}

as the graph of a possibly multi-valued function $u_H$ with respect to the $x_n$ coordinate

\begin{equation*}
\Gamma := \{ (x,x_{n+1},t) | \; x_n= u_H(x_1, x_2, \dots, x_{n-1}, x_{n+1}, t)  \}.
\end{equation*}

If we denote with $(y_1, \dots, y_n) $ the coordinates $(x_1, x_2, \dots, x_{n-1}, x_{n+1})$, then it is possible to have an expression of $\nabla u_S$ in terms of $\nabla u_H$ for a point on the graph $\Gamma$.
If $y' = (y_1, \dots, y_{n-1}) \in \mathbb{R}^{n-1}$ and $\nabla_{y'} $ denotes the gradient with respect to $y'$, we  define

\begin{equation*}
\begin{aligned}
F : \, &\mathbb{R}^n \rightarrow \mathbb{R}^n \\
& y \mapsto (y',\Bar{u}(y)).
\end{aligned}
\end{equation*}

Now, with this notation, we derive the expression
\begin{equation*}
x_{n+1}= u_S(y', u_H(y)) = (u_S \circ F) (y),
\end{equation*}

thus differentiating in the last coordinate gives

\begin{equation*}
1 = \partial_{x_{n+1}} (u_S \circ F) (y) = \sum_{i=1}^{n} \partial_{x_i} u_S(F(y)) \, \partial_{x_{n+1}} F_i(y) = \partial_{x_n} u_S \cdot \partial_{x_{n+1}} u_H
\end{equation*}

which implies

\begin{equation*} \displaystyle
\partial_{x_n} u_S = \frac{1}{\partial_{x_{n+1}} u_H}.
\end{equation*}

On the other hand, instead, for $j=1,\dots, n-1$ we have

\begin{equation*}
0 = \sum_{i=1}^{n} \partial_{x_i} u_S \, \partial_{x_j} F_i = \partial_{x_j} u_S + \partial_{x_n} u_S \cdot \partial_{x_j} u_H  \Longrightarrow  \partial_{x_j} u_S = - \frac{\partial_{x_j} u_H}{\partial_{x_{n+1}} u_H },
\end{equation*}

and thus

\begin{align*}
&\nabla u_S = - \frac{1}{\partial_{y_n}u_H}\left( \nabla_{y'} u_H, -1 \right), \quad \partial_{t} u_S = - \frac{\partial_{t}u_H}{\partial_{y_n} u_H}, \\
&D^{2}u_S= -\frac{1}{\partial_{y_n}u_H} A^T(\nabla u_H) \, D^2 u_H \, A(\nabla u_H),
\end{align*}

where, $A(\nabla u_H)$ is a square matrix which coincides with the identity except on the $n$-th row where the entries are given by the right hand side of $\nabla u_S $ above. Then, we can transform the Stefan problem (<ref>) in terms of $u_S$ to a quasilinear parabolic equation with respect to $u_H$  with oblique derivative boundary condition of the type

\begin{equation}\label{eq:HodographStefan}
\begin{cases}
\partial_{t} u_H - \tr( \Bar{A}(\nabla u_H) D^2 u_H) + \partial_{y_n}u_H \,  f(y',u_H,t)=0 & \text{in} \,\, \{ y_n > 0 \},\\
\partial_{t} u_H = g(\nabla  u_H)& \text{on} \,\,  \,\, \{ y_n = 0 \},
\end{cases}
\end{equation}

where, taking $p=(p',p_n) \in \R^n$

\begin{equation} \label{eq:g}
g(p)= - \frac{1}{p_n} \left( 1+ \nor{p'}^2 \right)
\end{equation}
and $\Bar{A}$ is the symmetrized of $A$:

\[
\Bar{A}(p)= \begin{pNiceArray}{c|c}
\text{ \begin{Large} $I_{n-1}$ \end{Large}} &   \displaystyle -\frac{p'}{p_n}\\
\hline
\displaystyle -\frac{p'}{p_n} &  \displaystyle \frac{1}{p_n^2} \left( 1+\nor{p'} \right)
\end{pNiceArray},
\]

with eigenvalues

\begin{equation*}
\spec(\Bar{A}) = \left\{1, \frac{1  + \nor{p}^2 - \sqrt{\left(1 + \nor{p}^2\right)^2 - 4p_n^2}}{2p_n^2},\frac{1  + \nor{p}^2 +
\sqrt{\left(1 + \nor{p}^2 \right)^2 - 4p_n^2}}{2p_n^2}\right\}.
\end{equation*}

The matrix $\Bar{A}$ is positive definite with the assumption that exists a $K>0$ such that $p \in \mathcal{R}_K$ with

\begin{equation*}
\mathcal{R}_K:= B_{K}(0) \cap \{ y_n  \geq K^{-1}\} \subset \mathbb{R}^{n}.
\end{equation*}

Notice also that if $p\in \mathcal{R}_K$ then $\partial_{y_n} g(p) > 0 $.

If $K > 0$ and $p \in $ is chosen large enough, we may assume $\tr(\Bar{A}(p)M) $ is uniformly elliptic in $M$ with ellipticity constants $K$ and $K^{-1}$.
Furthermore, we can suppose that

\begin{equation} \label{eq:boundg}
\nor{\nabla_{p} \tr(\Bar{A}(p)M)} \leq K \nor{M}, \quad \norm{g}_{C^1} \leq K, \quad \partial_{p_n} g(p) \geq K^{-1}.
\end{equation}

The hodograph transform is very useful because it flattens the free boundary, and thus the problem is transformed from a free boundary problem into one with fixed boundary but with possible multi-valued solution. This issue of multi-valued functions becomes the main subtlety with the transformed equation (<ref>).

Indeed one must adapt the definition of viscosity solution to include possibly multi-valued functions as solutions of the following fully nonlinear problem

\begin{equation}\label{eq:FullynonlinearProb}\begin{cases}\partial_t u = F(\nabla u, D^2u)  &\text{in $\mathcal C_{\lambda},$}\\
\partial_t u =g( \nabla u)& \text{on $\mathcal F_\lambda$}.
\end{cases}\end{equation}

We say that a  continuous function  $u : \overline{\mathcal{C}_\lambda} \rightarrow \R$ is a viscosity subsolution to (<ref>) if $u$ cannot be touched by above at points $\mathcal{C}_\lambda \cup \mathcal{F}_\lambda$ by a strict $C^2$ supersolution $\varphi$ of (<ref>). More precisely, we require that there do not exist points $(x_0,t_0) \in \mathcal{C}_\lambda \cup \mathcal{F}_\lambda $ and test functions $\varphi \in C^2(Q_r(x_0) \times (t_0-r^2,t_0])$ that satisfy

\begin{equation}\label{eq:TestFullynonlinear}\begin{cases}
\partial_t \varphi > F(\nabla \varphi, D^2\varphi)  &\text{in $Q_r(x_0) \times (t_0-r^2,t_0],$}\\
\partial_t \varphi >g( \nabla \varphi)& \text{on $\mathcal{F}_\lambda \cap  \left(Q_r(x_0) \times (t_0-r^2,t_0]\right)$,}
\end{cases}\end{equation}

such that

\begin{equation} \label{eq:TestFullynonlinear2}
u(x_0,t_0)= \varphi(x_0,t_0), \; u \leq \varphi \; \text{in  $Q_r(x_0) \times (t_0-r^2,t_0].$}
\end{equation}

Similarly we can define viscosity supersolutions and viscosity solutions to (<ref>).
These definitions can be generalized for multi-valued functions requiring that the comparison with a (single-valued) test functions $\varphi$ holds for all possible values.

Let $u : \overline{\mathcal{C}_\lambda} \rightarrow \R$ be a multi-valued function with compact graph in $\R^{n+2}$. We say that $u$ is a viscosity subsolution to (<ref>) if the definition above holds and (<ref>) is understood as $\varphi(x_0,t_0) \in u(x_0,t_0)$, while the inequality $u \leq \varphi$ in $Q_r(x_0) \times (t_0-r^2,t_0]$ means that $u(x,t) \leq \varphi(x,t)$ for all possible values of $u$ at $(x,t)$, and for all $(x,t) \in Q_r(x_0) \times (t_0-r^2,t_0]. $

Using the flatness assumption on the free boundary, we can localize the set where $u_H$ is possibly multi-valued. We will prove that $u_H$ is single-valued in the set $\{y_n \geq C_0\varepsilon_1\lambda\}$.  Let $t\in [t_0-\lambda^2,t_0+\lambda^2]$, then using the bounds for $|a'|$ and $|b'|$ in the assumptions of Theorem <ref>, we get

\begin{equation*}
\begin{aligned}
&\nor{a_{n}(t)-a_{n}(t_0)} \leq \nor{a_{n}'(t_1)} \nor{t-t_0} \leq c_0  \\
&\nor{b(t)-b(t_0)} \leq \nor{-a_{n}(t_2)} \nor{t-t_0}  \leq K\lambda^{2}
\quad \text{for }t_1,t_2 \in ]t_0,t[.
\end{aligned}
\end{equation*}

If $\lambda \leq \varepsilon_1, c_0 \leq \varepsilon_1$ then in $Q_{\lambda}(x_0) \times [t_0-\lambda^2, t_0 + \lambda^2]$ we have

\begin{align*}
&\nor{u_S-a_n(t_0) (x_n-b(t_0))^{+}} \\ \leq \,&\nor{u_S - a_n(t)(x_n - b(t))^+} + \nor{a_n(t_0)} \nor{(x_n -b(t))^{+} -(x_n -b(t_0))^{+}  } \\
+\,&\nor{a_n(t)-a_n(t_0)} \nor{(x_n-b(t))^{+}}\\
\leq\,& C \lambda \varepsilon_1 + K \nor{(x_n -b(t))^{+} -(x_n -b(t_0))^{+}  } + c_0 \lambda \leq  C \lambda \varepsilon_1 + K \nor{b(t) - b(t_0)} + c_0 \lambda \leq C \varepsilon_1 \lambda.
\end{align*}

If $(x_0,t_0)$ is such that $x_0 \in Q_{\lambda}, t_0 > -c\lambda$ and $C_0 \varepsilon_1 \lambda \leq u_S(x_0,t_0) \leq c\lambda$, for some $C_0$ to be chosen large. Note that $(u_S-a_n(t_0)(x_n - b(t_0))^+)$ still solves the heat equation with the same source term in cylinders away from the free boundary.

We start by  estimating $d:=\textit{dist}((x_0,t_0),\partial \{u>0\})$. Since

\begin{align*}
C_0\varepsilon_1 \lambda\leq u(x_0,t_0)\leq C \varepsilon_1 \lambda+a_n(t_0)((x_0)_n-b(t_0))^+
\end{align*}

we get

\begin{align*}
(x_0)_n\geq b(t_0)+(C_0-C)\lambda \varepsilon_1 a_n^{-1}(t_0).
\end{align*}

For the same time $t_0$, let $\bar x$ be a point in the free boundary, that is, $u(\bar x, t_0)=0$. If $\bar x \not \in Q_\lambda(x_0)$ then we are done because $d\geq \lambda$. Otherwise, it holds

\[
a_n(t_0)(\bar x_n-b(t_0))^+\leq C \varepsilon_1 \lambda
\]

and thus

\[
\bar x_n\leq b(t_0)+a_n^{-1}(t_0)C\varepsilon_1\lambda.
\]

\[
d\geq |(x_0)_n-\bar x_n|\geq (C_0-2C)\varepsilon_1\lambda a_n^{-1}(t_0).
\]

Combining this with the interior gradient estimates, see <cit.>, we get that there exists $C_e$ which depends only on $n$, such that

\begin{align*}
&\nor{\nabla u_S (x_0,t_0) - a_n(t_0)e_n} \\
\leq&\, C_e d^{-1}\left(\norm{u_S-a_n(t_0)(x_n - b(t_0))^+}_{L^{\infty}(Q_{\lambda})}+ \lambda \norm{f}_{L^{\infty}(Q_{\lambda})} \right)\\
\leq &\, C_e(C_0-2C)^{-1}(\varepsilon_1\lambda)^{-1}a_n(t_0)\left(C\varepsilon_1\lambda+\lambda \varepsilon_1^2\right)\\
=&\, C_e\frac{C+\varepsilon_1}{C_0-2C}a_n(t_0)
\end{align*}

Taking $C_0=(4C_e+2)C$ yields

\begin{align}\label{eq:nablaus}
\nor{\nabla u_S (x_0,t_0) - a_n(t_0)e_n}\leq \frac{1}{2}a_n(t_0),
\end{align}

which implies $|\nabla u_S(x_0,t_0)|\geq a_n(t_0)/2\geq (2K)^{-1}$ and thus by the inverse function theorem, the hodograph transform $u_H$ is well-defined in such points. Actually, we get $\nabla u_S(x_0,t_0)\in \mathcal{R}_{2K}$ which in terms of the hodograph transform still reads $\nabla u_H(x_0,t_0)\in \mathcal{R}_{2K}$.

Exploiting the inverse theorem, the main hypothesis of Theorem <ref> can be written in terms of the hodograph as

\begin{align*}
&\nor{u_H - (\bar a_n(t)y_n+\bar b(t))} \leq C\varepsilon_1 \bar \lambda\quad,   \textit{ in } \mathcal{C}_{\bar \lambda},  \\
&\bar b'(t)=g(\bar a_n(t) e_n), \quad  K^{-1} \leq \bar a_n \leq K, \quad \nor{\bar a_n'} \leq c_1 \bar \lambda^{-2},
\end{align*}

with $\bar \lambda\leq \lambda_1$ for some $\bar \lambda_1>0$ small.

Therefore we can reduce the initial problem of obtaining regularity of the free boundary, to proving an improvement of flatness directly to the hodograph transform $u_H$.

For agility of the discussion from now on we will denote by $u$ the transform function $u_H$ and we keep referring by $u_S$ the solution of (<ref>). Moreover, we drop the bar from $\bar A, \bar \lambda, \bar a, \bar b$ and we use $x$ rather than $y$.

§.§ The inhomogeneous $H$-property
In this Subsection, we explore an essential Harnack-type property of solutions to equation (<ref>). For this purpose, we introduce

\begin{equation}
l_{a,b}= a \cdot x + b
\end{equation}

where $a\in \mathbb{R}^{n}$, $\nor{a} \leq K$, $b \in \mathbb{R}$. Then  $\omega =u- l_{a,b}$ solves

\begin{equation*}
\partial_{t} \omega - \tr( \Bar{A}(\nabla \omega + a) D^2 \omega) + (\partial_{y_n}\omega + a_n) \,  f(y',\omega + a \cdot x - b,t)=0.
\end{equation*}

By uniform ellipticity, we have that

\begin{equation*}
\begin{aligned}
\partial_{t} \omega -\Mp(D^2 \omega) + \partial_{y_n}\omega f + a_n f \leq 0
\end{aligned}
\end{equation*}

and thus

\begin{align}\label{eq:M+omega}
\partial_{t} \omega -\Mp(D^2 \omega) - \nor{\nabla \omega} \norm{f}_{L^{\infty}} + a_n f \leq 0.
\end{align}

Analogously for $\Mm$,

\begin{equation}\label{eq:M-omega}
\partial_{t} \omega - \Mm(D^2 \omega) + \nor{\nabla \omega} \norm{f}_{L^{\infty}} + a_n f \geq 0
\end{equation}

Note that for every $Q_r(x_0)\subset \mathcal{C}_\lambda$ with $r>C\varepsilon_1\lambda$  we have that $u$ (and thus $\omega$) is single valued in $Q_{r/2}(x_0)\times [t_0-r^2,t_0+r^2]$. Thus by the Harnack inequality for the extremal equations (<ref>) and (<ref>) obtained in  <cit.>, we have

\begin{equation*}
\sup_{\substack{Q_{r/2}(x_0) \times [t_0-r^2, t_0]}} \omega \leq C_H \left( \inf_{\substack{Q_{r/2}(x_0) \times [t_0+\frac{r^2}{2}, t_0+r^2]}} \omega + \nor{a_n} \norm{f}_{L^{\infty}}\right)
\end{equation*}

where $C_H=C_H(n,K,\norm{f}_{L^{\infty}})$. Therefore, if
$$u - l_{a,b}=\omega \geq 0 \qquad \text{ in } Q_{r}(x_0) \times [t_0-r^2, t_0+r^2], \text{ for }r \geq C \varepsilon_1 \lambda$$

\begin{equation*}
(u - l_{a,b})(x_0,t_0)=\omega(x_0,t_0) \geq \mu, \qquad \text{for some } \mu \geq 0,
\end{equation*}

\begin{equation}
u - l_{a,b}=\omega \geq C_{H}^{-1}\inf_{\substack{Q_{r/2} \times [t_0+r^2/2, t_0+r^2]}} \omega  - \nor{a_n} \norm{f}_{L^{\infty}} \geq c \mu - K \norm{f}_{L^{\infty}} \geq c\mu - K\varepsilon_1^{2}.
\end{equation}

Similarly holds an analogous statement for $u \leq l_{a,b}$.
Now we are ready to state the $H(\sigma)$ property for general $\sigma >0$ small, which we have just proved for $\sigma=c\varepsilon_1$.
Given $\sigma$ positive small constant, we say that $u$ has property $H(\sigma)$ in $\mathcal{C}_\lambda$ if for any $r \in [\sigma \lambda, \lambda]$ and if

\begin{equation*}
u - l_{a,b} \geq 0 \qquad \text{ in } Q_{r}(x_0) \times [t_0-r^2, t_0+r^2] \subset \mathcal{C}_{\lambda},
\end{equation*}

\begin{equation*}
(u - l_{a,b})(x_0,t_0)\geq \mu, \qquad \text{for some } \mu \geq 0,
\end{equation*}

\begin{equation*}
u-l_{a,b} \geq c\mu \quad \text{ in } Q_{\frac{r}{2}}(x_0) \times [t_0+\frac{r^2}{2}, t_0+r^2].
\end{equation*}

Note that if $\sigma\leq\bar \sigma$ and $u$ satisfies $H(\sigma)$, then $u$ also satisfies $H(\bar \sigma)$. By the discussion above and Remark <ref>, we know that $u$ satisfies $H(C\varepsilon)$. Taking $\varepsilon$ small enough, we get that $u$ also satisfies $H(\varepsilon^{1/2})$.

§.§ Statement and discussion on main result

Let $\ell_{a,b}$ be the linear approximation of $u$ in $x$ for each $t$ fixed defined as

\begin{equation*}
\ell_{a,b}(x,t) := a(t) \cdot x + b(t)
\end{equation*}

\begin{equation*}
a(t):= (a_1, \dots, a_{n-1}, a_n(t))= (\hat{a},a_n(t)), \; \hat{a} \in \mathbb{R}^{n-1}
\end{equation*}

and in order to satisfy boundary condition of (<ref>) we further ask

\begin{equation} \label{eq:b1}
b'(t) = g(a(t)).
\end{equation}

The flatness assumption and Remark <ref> assure that

\begin{equation} \label{eq:flat}
\nor{u - (a_n(t_0) x_n + b(t_0)) }\leq C \varepsilon_1 \lambda \quad \text{ in } Q^{+}_{\lambda} \times [t_0-\lambda^2,t_0+\lambda^2].
\end{equation}

Now we are able to state the main result, concerning the improvement of flatness of $u$.

Fix $K > 0$ large, $u$ is a viscosity solution to the problem (<ref>)
and  assume that $\norm{f}_{L^{\infty}(\Bar{\mathcal{C}}_\lambda)} \leq \varepsilon_1^2$. Let $u$ be possibly a multi-valued function, satisfying the property $H(\varepsilon^{1/2})$ and

\begin{equation} \label{Hp:improv}
\begin{aligned}
\nor{u-\ell_{a,b}} \leq \varepsilon \lambda \text{ in } \bar{\mathcal{C}}_{\lambda} \text{ with } b'(t)=g(a(t)), \\
a(t) \in \mathcal{R}_{K}, \nor{a_n'(t)} \leq \delta \varepsilon \lambda^{-2},
\end{aligned}
\end{equation}

\begin{equation*}
\varepsilon \leq \varepsilon_{1}, \quad \lambda \leq \lambda_{1}, \quad \lambda \leq \delta \varepsilon.
\end{equation*}

Then there exists $\ell_{\Tilde{a},\Tilde{b}}$ such that

\begin{equation*}
\nor{u- \ell_{\Tilde{a},\Tilde{b}}} \leq \frac{\varepsilon}{2} \tau \lambda \; \text{ in } \Bar{\mathcal{C}}_{\tau \lambda }, \quad \Tilde{b}'(t)=g(\Tilde{a}(t)),
\end{equation*}

\begin{equation*}
\nor{a(t) - \Tilde{a}(t)} \leq C \varepsilon, \quad \nor{\Tilde{a}_n'(t)} \leq \frac{\delta \varepsilon}{2} (\tau \lambda)^{-2}.
\end{equation*}

Here the constants $\varepsilon_{1}, \lambda_{1}, \delta, \tau > 0$ small and $C$ large depend only on $n$ and $K$.

The improvement of flatness in this case is given by the Proposition <ref> that we are going to prove.

\begin{equation}\label{eq:w_error}
u(x,t) = \ell_{a,b}(x,t) + \varepsilon \lambda w\left( \frac{x}{\lambda}, \frac{t}{\lambda}\right), \quad (x,t) \in \mathcal{C}_{\lambda},
\end{equation}

where $w(x,t)$ is the rescaled error function, defined in $\mathcal{C}_1,$ possibly multi-valued near $\{ x_n = 0 \}$ and $\nor{w} \leq 1$ in $\mathcal{C}_1$, since $\nor{u - \ell_{a,b}} \leq \varepsilon \lambda.$

The linearized problem associated to (<ref>) becomes

\begin{equation}\label{eq:LinStefanNH}
\begin{cases}
\lambda^2 \varepsilon^{-1} a_n'(\lambda t) x_n + \lambda \varepsilon^{-1} b'(\lambda t) + \lambda \partial_t w(x,t) - \tr(A(a(\lambda t) + \varepsilon \nabla w)) D^2 w ) +\\ \left( \lambda \partial_{x_n}w+\lambda \varepsilon^{-1} a_n(\lambda t) \right) f(x',\ell_{a,b}+\varepsilon\lambda w,t) =0  & \text{in} \,\, \mathcal{C}_1,\\
\varepsilon^{-1} b'(\lambda t) + \partial_{t} w = \varepsilon^{-1}  g(a(\lambda t) + \varepsilon \nabla w)& \text{on} \,\,  \,\, \mathcal{F}_1.
\end{cases}
\end{equation}

From (<ref>) and (<ref>) we can write the boundary condition in a convenient way as

\begin{equation*}
\partial_{t} w =  \displaystyle\frac{(1 + \nor{\hat{a}}^2) \partial_{x_n}w}{a_{n}(\lambda t) (a_{n}(\lambda t) + \varepsilon \partial_{x_n}w)} - \frac{\varepsilon \nor{\nabla_{x'}w}^2+ 2\hat{a}\nabla_{x'}w}{a_{n}(\lambda t) + \varepsilon \partial_{x_n} w} \quad \text{on} \,\,  \,\, \mathcal{F}_1.
\end{equation*}

Now letting $\varepsilon \rightarrow 0, \delta \rightarrow 0 $ in (<ref>) and using that $\nor{a'} \leq \delta \varepsilon \lambda^{-2},a(\lambda t)  \in \mathcal{R}_K$ and $\lambda \varepsilon^{-1} \leq \delta \rightarrow 0$, we formally obtain that our equation becomes close to

\begin{equation}\label{eq:LimitStefanNH}
\begin{cases}
\lambda \partial_t v = \tr(A_{\lambda}(t) D^2 v )& \text{in} \,\, \mathcal{C}_1,\\
\partial v_t = \gamma_{\lambda} (t) \cdot \nabla v & \text{on} \,\,  \,\, \mathcal{F}_1,
\end{cases}
\end{equation}

\begin{equation*}
A_{\lambda}(t):= A(a(\lambda t)), \quad \gamma_{\lambda} (t) :=  \nabla g (a(\lambda t)).
\end{equation*}

Using the fact that $A,g \in C^2(\mathcal{R}_K)$  we have following conditions on $A_\lambda$ and $\gamma_\lambda$

\begin{align*}
K^{-1}I \leq A_\lambda(t) \leq K I,  \quad
\nor{A_\lambda'(t)} \leq \lambda^{-1}, \quad  \gamma_\lambda \in \mathcal{R_K},
\quad
\nor{\gamma_\lambda'(t)} \leq \lambda.
\end{align*}

§ PRELIMINARY RESULTS

In this section, we recall regularity results for solutions of (<ref>) obtained for the homogeneous case and then we obtain new comparison results relating this equation to (<ref>).

§.§ Regularity results for the limiting profile

The following regularity results can be found in <cit.>.

[Interior estimates]
Let $v$ be a viscosity solution to (<ref>) such that
$\|v\|_{L^\infty} \leq 1.$ Then

$$| \nabla v|, \, \, \,  |D^2 v| \, \, \leq C \quad \quad \textit{in} \quad \mathcal C_{1/2},$$

and there exists $l_{\bar a,\bar b}$  such that for every $\rho \leq 1/2$,

$$|v - l_{\bar a, \bar b}|\leq C \rho^{1+\alpha} \quad \quad \text{in $\mathcal C_{\rho}$},$$

$$\bar b'(t)= \gamma_\lambda(t) \cdot \bar a, \quad \quad |\bar a'_n| \leq C \rho^{\alpha-1}\lambda^{-1}, \quad \quad |\bar a| \le C,$$
with $\alpha$, $C$ universal.

[The Dirichlet problem]
Let $\phi$ be a continuous function on $\p_D \mathcal C_1$. Then there exists a unique classical solution $v \in C^{2,1} (\mathcal C_1) \cap C(\bar {\mathcal C_1})$ to the Dirichlet problem (<ref>) with $v=\phi$ on $ \p_D \mathcal C_1$. Moreover,

$$|\nabla v|, |D^2 v| \le C(\sigma) \|v\|_{L^\infty} \quad \textit{in} \quad C_1^\sigma:=\{ d_\lambda((x,t), \p_D \mathcal C_1) \ge \sigma\},$$

and if $\phi$ is $C^\alpha$ with respect to the distance $d_\lambda,$ then $v$ is also $C^\alpha$ up to the boundary and

\alpha_{d_\lambda}} \le C \|\phi\|_{C^\alpha_{d_\lambda}},$$

with $C(\sigma)$, $C$ universal constants, independent of $\lambda$.

§.§ Comparison principle

We establish that solutions of (<ref>) are close to solution of $\eqref{eq:LimitStefanNH}$. This is a consequence of the comparison principle stated in the following proposition

Let $v \in C^{2}(\Bar{\Omega})$ with $\Omega \subset \mathcal{C}_1$ satisfy

\begin{equation*}
\nor{\nabla v}, \, \nor{D^2 v} \leq M,
\end{equation*}

for some large constants $M$ and

\begin{equation}\label{eq:CompareLimitStefanNH}
\begin{cases}
\lambda \partial_t v \leq \tr(A_{\lambda}(t) D^2 v ) - C \delta& \text{in} \,\, \Omega,\\
\partial_t v \leq \gamma_{\lambda} (t) \cdot \nabla v - \delta & \text{on} \,\,  \,\, \mathcal{F}_1 \cap \Bar{\Omega},
\end{cases}
\end{equation}

with $A_{\lambda}(t), \gamma_{\lambda}(t)$ as above.
Then $v$ is a subsolution to (<ref>), as long as $C$ is sufficiently large, universal, and $\varepsilon \leq \varepsilon_1(\delta, M)$. In particular, if

\begin{equation*}
v \leq w \; \text{ on } \overline{\partial \Omega \setminus (\{ t=0\} \cup \{x_n = 0 \})}
\end{equation*}

\begin{equation*}
v \leq w \text{ in } \Omega.
\end{equation*}

Similarly holds an analogous result for supersolution replacing $\leq$ with $\geq$ and the $-$ in (<ref>) by $+$.

Note that since

\begin{equation*}
\lambda\delta \nor{a'(\lambda t)} + \delta\nor{b'(\lambda t)} \leq \delta(\lambda K +  K),
\end{equation*}

\[
\left| \lambda \partial_{x_n}v+\lambda \varepsilon_1^{-1} a_n(\lambda t) \right|\, |f(x',v,t)|\leq \left( \lambda M+\delta K\right)\|f\|_\infty
\]

\begin{equation*}
\left|\tr (A(a(\lambda t) + \varepsilon \nabla v) - A(a(\lambda t)) D^2 v )\right| \leq M K \varepsilon_1 \nor{D^2 v} \leq K \varepsilon_1 M^2,
\end{equation*}

since $a \in \mathcal{R}_K$.
Now, obtain that v is a solution of the first equation of (<ref>) if we choose $C > 0$ large. In fact

\begin{align*}
\lambda \partial_t v &\leq \tr(A_{\lambda}(t)D^2 v) - C\delta \leq \tr(A(a(\lambda t + \varepsilon_1 \nabla v))D^2 v) + K \varepsilon_1 M^2 - C \delta \\
& \leq \tr(A(a(\lambda t + \varepsilon_1 \nabla v))D^2 v) + K \varepsilon_1 M^2 - C \delta - (\lambda^2 \varepsilon_1^{-1} a_n'(\lambda t) x_n + \lambda \varepsilon_1^{-1} b'(\lambda t))\\ &+ \delta (\lambda K + K) - (\lambda \partial_{x_n} w + \lambda \varepsilon_1^{-1})f + (\lambda M + \delta K) \norm{f}_{L^{\infty}}.
\end{align*}

Then if $C$ is chosen such that

\begin{equation*}
C \delta \geq K \varepsilon_1 M^2 + \delta (\lambda K + K) +(\lambda M + \delta K) \norm{f}_{L^{\infty}}
\end{equation*}

we have

\begin{equation*}
\lambda \partial_t v \leq \tr(A(a(\lambda t + \varepsilon_1 \nabla v))D^2 v) - \lambda^2 \varepsilon_1^{-1} a_n'(\lambda t) x_n - \lambda \varepsilon_1^{-1} b'(\lambda t)- (\lambda \partial_{x_n} w + \lambda \varepsilon_1^{-1})f
\end{equation*}

which is exactly the first of (<ref>)
The same holds for the boundary condition using

\begin{equation*}
\nor{g(a(\lambda t) + \varepsilon_1 \nabla v) -g(a(\lambda t)) - \varepsilon_1 \nabla g(a(\lambda t)) \cdot \nabla v } \leq C \varepsilon_1^2 M^2.
\end{equation*}

Thus it follows  that if the rescaled error $w$ is close to a $C^2$ solution $v$ of (<ref>) on the Dirichlet boundary of $\Omega \subset \mathcal{C}_1$, then $v$ and $w$ remain close to each other in the whole domain $\Omega$

Let $w$ be a solution to (<ref>) and $v \in C^{2}$ be a solution to (<ref>)  in a domain $\Omega \subset \mathcal{C}_1$, with

\begin{equation*}
\nor{\nabla v}, \nor{D^2 v} \leq M.
\end{equation*}

If $\varepsilon \leq \varepsilon_1(\delta,M) $ and

\begin{equation*}
\nor{v-w} \leq \sigma \quad \text{ on } \overline{\partial \Omega \setminus (\{ t=0\} \cup \{x_n = 0 \})}
\end{equation*}

\begin{equation*}
\nor{v-w} \leq \sigma + C\delta \quad \text{in } \Omega.
\end{equation*}

This corollary follows from Proposition <ref> if instead of $v$ we take

\begin{equation*}
v \pm \left( C \delta(x_n^2 -t -2)-\sigma \right).
\end{equation*}

§ THE OSCILLATION DECAY PROPERTIES OF THE ERROR FUNCTION $W$

In this section it will be useful to consider the following distance in $\R^{n+1}$
which interpolates between the parabolic and standard distances depending on how far we are from the boundary $\{x_n=0\}$.

\begin{align*}
d_\lambda((x,t),(y,s)):=\min\left\{|x'-y'|+|x_n-y_n|+\lambda^{-\frac{1}{2}} |t-s|^{\frac{1}{2}}, |x'-y'|+|x_n| + |y_n| +|t-s|    \right\}.
\end{align*}

We call by $\mathcal{B}_{\lambda,r}$ and $C_{d_\lambda}^{\alpha}$ respectively the balls and Hölder space induced by the metric $d_\lambda$.

Let $w$ be given by (<ref>). We want to prove the following.

Assume $u$ satisfy the hypotheses of Proposition <ref> and let $w$ be given by (<ref>) with $-1\leq w\leq 1$. Then

\[
\osc_{\mathcal{C}_{1/2}}w\leq 2(1-c),
\]

with $c$ small universal, provided that $\delta \leq K^{-1}$ and $\varepsilon \leq c\delta$.

In order to prove Proposition <ref> we consider the function $\bar w $ defined as

\begin{equation}
\bar w:=w+1+C\delta\left(2+t-x_n^2\right)\geq 0,
\end{equation}

From (<ref>) follows immediately that $\bar w $ satisfies

\begin{equation} \label{eq:errorshifted}
\begin{cases}
\lambda \bar w_t = &\tr(A(a(\lambda t) + \varepsilon(\nabla \bar w + 2C\delta x_n e_n))(D^2 \bar w + 2C\delta (e_n \otimes e_n))) - \lambda^2 \varepsilon^{-1} a_n'(\lambda t) x_n + \\
&- \lambda \varepsilon^{-1} b'(\lambda t) + \lambda (\p_{x_n} \bar w + 2C \delta x_n + \varepsilon^{-1} a_n (\lambda t)) f +C \delta \lambda  \qquad \text{in }\mathcal{C}_1; \\
\bar w_{t} = &-  \varepsilon^{-1} g(a(\lambda t)) + \varepsilon^{-1} g(a(\lambda t) + \varepsilon(\nabla \bar w + 2C \delta x_n e_n)) +C \delta \qquad \text{on }\mathcal{F}_1.
\end{cases}
\end{equation}

Note that considering $v \in C^{2}(\mathcal{C}_1)$ with $\nor{\nabla v}, \nor{D^2 v} \leq M$  a solution of

\begin{equation*}
\begin{cases}
\lambda \partial_t v \leq \Mm(D^2 v) & \text{in} \,\, \mathcal{C}_1,\\
\partial_t v \leq K^{-1} v_n^+ - Kv_n^- - K \nor{\nabla_{x'} v}  & \text{on} \,\,  \,\, \mathcal{F}_1,
\end{cases}
\end{equation*}

is a subsolution to (<ref>) provided that $\delta \leq K^{-1}, \varepsilon \leq c \delta$.

Proof of Remark <ref>
Exploiting the properties of Pucci operators

\begin{equation*}
\begin{aligned}
\Mm(D^2 v + 2C\delta e_n \otimes e_n) \geq \Mm(D^2 v) + \Mm(2C\delta e_n \otimes e_n) \geq  \Mm(D^2 v) + \frac{2C\delta}{K}.
\end{aligned}
\end{equation*}

\begin{equation*}
\begin{aligned}
\lambda \p_t v \leq \Mm(D^2 v)  &\leq \Mm(D^2v  + 2C\delta e_n \otimes e_n) -\frac{2C\delta}{K}\\ &\leq \tr(A(a(\lambda t) + \varepsilon( \nabla v + 2C\delta x_n)) (D^2 v + 2C \delta (e_n \otimes e_n) ) ) - \frac{2C\delta}{K}.
\end{aligned}
\end{equation*}

Define the following quantity to lighten the notation

\begin{equation*}
H(v,\p_{x_n}v) :=  -\lambda^2 \varepsilon^{-1} a_n'(\lambda t) - \lambda \varepsilon^{-1} b'(\lambda t) + \lambda (\p_{x_n} v + 2C \delta x_n  + \varepsilon^{-1} a_n( \lambda  t) ) f +C\delta \lambda.
\end{equation*}

Then, as $\lambda \leq \lambda_1 \leq 1$ we have that $ H(v,\p_{x_n} v ) \geq - (3K + M + C\lambda) \delta $ which implies that

\begin{equation}
\lambda \p_t v  \leq \tr(A(a(\lambda t) + \varepsilon( \nabla v + 2C\delta x_n)) (D^2 v + 2C \delta (e_n \otimes e_n) ) )  + H(v,\p_{x_n} v)  - \frac{2C\delta}{K} + (3K + M +  C \lambda) \delta .
\end{equation}

Having $\delta \leq K^{-1}$ and choosing $C ( \geq K(3K+M))$ large enough leads to

\begin{equation}
\lambda \p_t v  \leq \tr(A(a(\lambda t) + \varepsilon( \nabla v + 2C\delta x_n)) (D^2 v + 2C \delta (e_n \otimes e_n) ) )  + H(v,\p_{x_n} v) .
\end{equation}

On the other hand for the boundary condition, we have

\begin{equation*}
g(a(\lambda t + \varepsilon \nabla v + 2 \varepsilon C \delta x_n )) \geq g (a(\lambda t))+ \nabla g (a(\lambda t)) \cdot (\varepsilon \nabla v + 2\varepsilon C \delta x_n) - \varepsilon^2 C \nor{\nabla v + 2C \delta x_n}^2
\end{equation*}

which leads to

\begin{align*}
&\nabla g(\lambda t) \cdot \nabla v \\
\leq\,& \varepsilon^{-1} \left( g(a(\lambda t + \varepsilon \nabla v + 2 \varepsilon C \delta x_n)) - g(a(\lambda t))\right) - 2C \delta x_n \nabla g(a(\lambda t)) \cdot e_n + C \varepsilon \left( \nor{\nabla v}^2 + 4 C^2 \delta^2\right)   \\
\leq \,& \varepsilon^{-1}g(a(\lambda t + \varepsilon \nabla v + 2 \varepsilon C \delta x_n)) - \varepsilon^{-1} b'(\lambda  t) + \delta \left(  2KC  + C\varepsilon  \delta^{-1} M^2 + 4 \varepsilon C\delta \right)
\end{align*}

assuming $\varepsilon \leq c\delta$

\begin{align*}
\p_{t} v \leq \,&K^{-1} v_n^+ - Kv_n^- - K \nor{\nabla_{x'} v} \leq \nabla g ( a (\lambda t)) \cdot \nabla v \\
\leq\,& \varepsilon^{-1}g(a(\lambda t + \varepsilon \nabla v + 2 \varepsilon C \delta x_n)) - \varepsilon^{-1} b'(\lambda  t) + C\delta
\end{align*}

where $C=C(M,K,n)$.

Now, we state an alternative version of  <cit.> involving $\bar w$ which describes its evolution over time when compared to an explicit barrier. We start by introducing such function considering $\Omega$ be a smooth domain in $\mathbb{R}^n$, $n \geq 2,$ such that

$$\bar Q^+_{3/4} \subset \bar \Omega \subset \bar Q^+_{7/8},$$

and call

$$T:= \{x_n=0\} \cap Q_{3/4} \subset \partial\Omega.$$

Define $\eta(x')$ a standard bump function supported on $Q'_{5/8}$ and equal 1 on $Q'_{1/2}$ (here the prime denotes cubes in $\mathbb{R}^{n-1}$).
Let $\phi$ satisfy

$$\mathcal M^-_{K}(D^2\phi)= 0 \quad \text{in $\Omega$},$$

$$\phi= 0 \quad \text{on $\partial \Omega \setminus T$}, \quad \phi=\eta \quad \text{on $T$},$$
and notice that $0 \leq \phi \leq 1$, $\phi \geq c$ on $Q_{1/2}^+$, and by Hopf lemma $\phi_n>0$ on $\{x_n=0\} \cap \{\phi=0\}$. Furthermore since $\Mm$ is concave, then by  <cit.> we have that $\phi \in C^{2,\alpha}(\overline{\Omega})$ with universal estimates.

Let $\bar w \geq 0$ satisfy (<ref>)
in the viscosity sense. If for some $t_0 \in (-1,0]$,

$$\bar w(x, t_0) \geq s_0 \, \, \phi(x) \quad \text{in $Q^+_1$}, \quad s_0 \geq 0,$$

$$\bar w(x,t) \geq s(t) \, \phi(x) \quad \text{in $Q^+_1 \times [t_0, 0],$}$$

$$s'(t)= - C_0 s(t), \quad s(t_0)=s_0, \quad \quad \textit{$C_0$ large universal}.$$

Moreover, if $s_0 \le c_0$ with $c_0$ small universal, and

\begin{equation}\label{eq:center} \bar w \left(\frac 1 2 e_n, t_0+\lambda/4 \right) \geq \frac {1}{2},\end{equation}

$$\bar w(x,t_0+\lambda) \geq (s_0 + c_0\lambda) \phi(x).$$

Note that $\bar w$ is possibly multi-valued near $\mathcal{F}_1$, therefore the inequalities hold for every possible value. The proof of this result is very similar to <cit.> since it still satisfies the interior Harnack inequality given by the $H$ property of $w$ after a rescaling in time to size $\lambda$. This fact together with Remark <ref> allows to perform the comparison with the explicit same barriers constructed in the proof given by the authors.

We are ready now to provide a proof for Proposition <ref>:

Proof of Proposition <ref>  We can assume that $w(e_n /2 , t_k + \lambda/4) \geq 0$ for more than half of the values of $t_k=-1+k\lambda \in [-1,-1/2]$ because if not we consider $-w$ that satisfies an equation for which is possible to replicate Lemma <ref>. Therefore, it holds

\begin{equation} \label{eq:morethanhalf}
\bar w(e_n /2 , t_k + \lambda/4) \geq 1,\quad \textit{ for more than half of the values of } \quad t_k\in [-1,-1/2].
\end{equation}

We want to use Lemma <ref> at each time $t_k$. We start by explaining the argument. Assume we have proved $\bar w\geq s_k\phi(x)$. If (<ref>) holds (good case), then we use the second part of Lemma <ref> to improve

\[
\]

On the other hand, if (<ref>) does not hold (bad case), then we have to decrease our control and get

\[
s_{k+1}\geq s_k(1-C_0\lambda).
\]

Since (<ref>) holds half of the time, the final control we get on $s_k$ will depend on the order of the times $t_k$ where (<ref>) holds. Clearly, the worst case scenario is when it holds for the first half of times $t_k$. We present only the estimate in this case.

We start with $s_0=0$. Then we apply the second part of Lemma <ref> $\bar k$ times, where $\bar k=1/(4\lambda)$ and get

\[
s_{\bar k}=\bar kc_o\lambda =\frac{c_0}{4}.
\]

Then, in the remaining $\bar k$ times we have the bad case, and we use the first part of the lemma. We get the final estimate

\begin{align*}
s_{2\bar k}\geq s_{\bar k}(1-C_0\lambda)^{\bar k}\geq \frac{c_0}{4}4^{-C_0/4}=:\bar c,
\end{align*}

provided we choose $\lambda\leq 1/(2C_0)$. Thus, the lemma implies

\begin{align*}
\bar w(x,t)>\bar c\, \phi(x),\quad  (x,t) \in Q^+_1\times[-1/2,0].
\end{align*}

Since $\phi>c$ in $Q_{1/2}^+$, we obtain the desired estimate.

Using the result of Proposition <ref> we can prove that $w$ has an almost Hölder modulus of continuity which will play a key role in establishing the improvement of flatness.
[Almost Hölder continuity for $w$]
Assume that $u$ satisfies the hypotheses of the Proposition <ref> and $w$ is always defined with (<ref>) as the error. Then for any $(x_0,t_0) \in \mathcal{C}_{\frac{1}{2}}$

\begin{equation} \label{eq:Pro:w_holder}
\osc_{\substack{\mathcal{B}_{\lambda,r}(x_0,t_0)}} w \leq C r^{\alpha}
\end{equation}

for $r > C(\delta) \varepsilon^{1/2}$ with $\delta \leq c' $ universal.

The proof follows from Proposition <ref> as in the homogeneous case.  For sake of completeness we provide the details of the proof which entails the need of working with the metric $d_\lambda$.

Note that we can rephrase the thesis of Proposition <ref> in terms of $u$ in $\mathcal{C}_{\lambda/2}$. In fact

\begin{equation*}
\osc_{\mathcal{C}_{\lambda /2}} (u- l_{a,b}) \leq 2\varepsilon \lambda (1-c).
\end{equation*}

In particular, from the proof of the same result we have shown that $(1-c)\varepsilon \lambda \leq u-l_{a,b} \leq \varepsilon \lambda$.
If we now consider

\begin{equation*}
\tilde w \left( \frac{x}{\lambda},\frac{t}{\lambda}\right):= \frac{u-l_{a,b}}{(1-c) \varepsilon \lambda} - \frac{c}{2-c}
\end{equation*}

it appears that $\tilde w \in \left( -\frac{2-c}{2-2c}, \frac{2-c}{2-2c}  \right)$. So, in order to prepare the iteration we are keen to take $w_1 := \tilde w \cdot \frac{2-2c}{2-c}$ which stays in $(-1,1)$. In this way, apart from relabeling the constants, we obtain that

\begin{equation*}
w_1 = \frac{u- \left(l_{a,b}+c\lambda \varepsilon/2\right)}{\lambda \left( \varepsilon (1-c)\right)} \quad \text{and} \quad  \nor{u - \left(l_{a,b} + \frac{c\lambda\varepsilon}{2}\right)} \leq \lambda \varepsilon (1-c).
\end{equation*}

Thus the hypotheses of Proposition <ref> are fulfilled in $\mathcal{C}_{\lambda/2}$ if one takes, instead of $\lambda$, $\varepsilon$, $\delta$, $a(t)$, $b(t)$, the quintet $\lambda/2$, $2(1-c)\varepsilon$, $\delta$, $a(t)$, $b(t)-c\varepsilon\lambda/2$ and this starts the iteration.
The original $H(\varepsilon^{1/2})$ property in $\mathcal{C}_{\lambda/2^{k}}$ scales to $H(2^k \varepsilon^{1/2})$. Hence, we are able to iterate Proposition $\ref{prop:wharnack}$ $k$ times as long as

\begin{equation*}
2^{k} \varepsilon^{1/2} \leq c'', \quad \delta \leq K^{-1}, \quad 2^{k} (1-c)^k \varepsilon \leq c \delta.
\end{equation*}

This means that we can continue the iteration $k$ times if $2^k \varepsilon^{1/2} \leq C \delta. $ In terms of $w$ we obtain that its oscillation in $\mathcal{C}_{2^{-k}}$ is bounded by $2(1-c)^k$. Now note that in order to make use of $H(c\varepsilon^{1/2})$ property of $u$ one considers parabolic cylinders centered in $(y,s)$ such that $\nor{y_n} \geq c\varepsilon^{1/2}$. Furthermore, the construction of $\mathcal{B}_{\lambda,r}(y,s)$ coincides with $\mathcal{C}_{r}(y,s)$ if $\nor{y_n} \leq r$, so the iteration leads to the desired result if the radius of the ball is greater than the distance of its center from $\{x_n = 0 \}$. For interior balls  $\mathcal{B}_{\lambda,r}(y,s)$ (say $\nor{y_n} \geq \varepsilon^{1/2}$ and $2r \leq \nor{y_n}$) the same diminishing behavior of the oscillation holds using directly the $H$ property for $w$.

[$w$ is close to a Holder function]
For every $\bar \varepsilon > 0$ there exists a $\bar \delta >0$ such that if $w$ satisfies the thesis of Proposition <ref> for $r \geq \bar \delta$, $-1\leq w \leq 1$, its graph is compact in $\mathbb{R}^{n+2}$ and is single-valued for $x_n \geq \bar \delta$ then there exists a (single-valued) function $\phi \in C_{d_{\lambda}}^{\alpha}( \mathcal{C}_{1/2})$ such that

\begin{equation} \label{eq:w_near_holder}
\nor{w- \phi} \leq \bar \varepsilon, \qquad \norm{\phi}_{C_{d_{\lambda}}^{\alpha}( \mathcal{C}_{1/2})} \leq C.
\end{equation}

In order to obtain this result we report hereafter the following lemma which is a slight variation of <cit.> for multivalued functions with compact graph.  In the following $d_{\mathcal{H}}(\cdot,\cdot)$ will denote the Hausdorff distance of two nonempty sets defined as

\begin{equation*}
d_{\h}(K_1,K_2) = \max \left\{  \sup_{ y\in K_1} d( y, K_2), \sup_{ y\in K_2} d( y, K_1) \right\}
\end{equation*}

where $K_1,K_2 \in \R^{m}$ are nonempty sets.
Moreover, if $(M,d)$ is a metric space, then $\mathcal{K}(M)$ denotes the set of all nonempty compact subsets of $M$. Furthermore, if $M$ is compact then so is $\mathcal{K}(M)$.

Let $\{w_k\}_{k \geq 1}$ a sequence of possibly multivalued continuous functions on  $\overline{\mathcal{C}_{1/2}} \subset \mathbb{R}^{n+1}$ such that

(i) $\nor{w_k} \leq 1$ on  $ \overline{\mathcal{C}_{1/2}}$ for all $k \geq 1$;

(ii) there exists a sequence $\delta_k \rightarrow 0$ as $k \rightarrow \infty$ such that $w_k$ are single-valued functions in the region where $x_n \geq \delta_k$ and exists a modulus of continuity $\omega : [0,\infty) \rightarrow [0,\infty) $ such that

\begin{equation} \osc_{\substack{\mathcal{B}_{\lambda,d_p((x,t),(y,s))}(x,t)}} w_k \leq  \omega(d_p((x,t),(y,s)))
\end{equation}

for all  $(x,t),(y,s) \in \overline{\mathcal{C}_{1/2}}$  such that $ d_p((x,t),(y,s)) \geq \delta_k$.

For all $k \geq 1$ let $G(w_k):=\left\{ (x,t,x_{n+1}) : (x,t) \in \overline{\mathcal{C}_{1/2}}, \, x_{n+1}= w_{k}(x,t) \right\}$ be graph in $\mathbb{R}^{n+2}$ assumed to be compact. Then exists a single-valued function $\phi : \overline{\mathcal{C}_{1/2}}\rightarrow [-1,1]$ such that

\begin{equation} \label{eq:Lm:w_approx}
d_{\mathcal{H}}\left(G(w_k), G(\phi) \right) \rightarrow 0 \; \text{for }k \rightarrow \infty.
\end{equation}

Moreover for all  $(x,t),(y,s) \in \overline{\mathcal{C}_{1/2}}$ we have

\begin{equation*}
\nor{\phi(x,t)-\phi(y,s)} \leq \omega\left( d_p((x,t),(y,s))\right).
\end{equation*}

Proof of Lemma <ref>
The proof of this result follows the lines of the proof of Lemma 9.1 in <cit.> being $G(w_k)$ a compact subset of  $\overline{\mathcal{C}_{1/2}} \times [-1,1]$.

By the continuity of $w_k$ and $(i)$, we have that $G_k\equiv G(w_k)$ is a compact subset of $\overline{\mathcal{C}_{1/2}} \times [-1,1]$. Since $\mathcal{K}(\overline{\mathcal{C}_{1/2}} \times [-1,1])$ is a compact metric space with the Hausdorff distance, there exists a compact set $K\subset \overline{\mathcal{C}_{1/2}} \times [-1,1]$ such that

\begin{equation}\label{eq:convhauss}
d_{\h}(G_k, K)\to 0 \quad\hbox{as}~ k \to \infty,
\end{equation}

up to a subsequence.
Next, we show that $K$ must be the graph of some function. If not, there exist two points $(x,t,p), (x,t,q) \in K$ with $p\neq q$. By (<ref>) for all $\varepsilon>0$, there is $k_0\geq 1$ such that $(x,t,p) \in B_\varepsilon(G_k)$, for all $k\geq k_0$, i.e.,
$d((x,t,p), G_k) \leq \varepsilon$ for all  $k\geq k_0$ (here $d$ denote the distance in $\mathbb{R}^{n+2}$).
In particular, since $G_k$ is compact, the distance is attained at some point $(x_k,t_k, p_k)\in G_k$. Hence,

\begin{equation*}
(x_k,t_k, p_k) \to (x,t,p) \quad\ \hbox{as}~ k \to \infty.
\end{equation*}

Similarly, there are points $(y_k,s_k, q_k)\in G_k$ such that

(y_k, s_k , q_k) \to (x,t, q) \quad\ \hbox{as}~ k \to \infty.

Consider first the case in which $x_n \neq 0$, choosing $k$ large enough we can always assume $(x_k)_n \geq \delta_k$.
We distinguish two cases: If  $d_p((x_k,t_k),(y_k,s_k))\geq \delta_k$  for infinitely many $k$'s, then by $(ii)$, we get

\begin{align*}
|p_k-q_k| & \leq  \osc_{\substack{\mathcal{B}_{\lambda,d_p((x_k,t_k),(y_k,s_k))}(x_k,t_k)}} w_k \leq  \omega (d_p((x_k,t_k),(y_k,s_k))).
\end{align*}

Letting $k\to\infty$, we get that $p=q$, which is a contradiction.
Otherwise, there is some $(z_k,r_k) \in \overline{\mathcal{C}_{1/2}} $ such that $\delta_k \leq |d_p((x_k,t_k),(z_k,r_k))|\leq 3\delta_k$ and $\delta_k \leq |d_p((z_k,r_k),(y_k,s_k))|\leq 3\delta_k$.

\begin{align*}
|p_k-q_k| \leq \osc_{\substack{\mathcal{B}_{\lambda,3\delta_k}(x_k,t_k)}} w_k  + \osc_{\substack{\mathcal{B}_{\lambda,3\delta_k}(y_k,s_k)}} w_k \leq 2 \omega (3\delta_k)
%|p_k-q_k| &\leq  |w_k(x_k,t_k)-w_k(z_k,r_k)| +  |w_k(z_k,r_k)-w_k(y_k,s_k)| \\
%&\leq  \omega(d_p((x_k,t_k),(z_k,r_k))) + \omega(d_p((z_k,r_k),(y_k,s_k)))\\
%&\leq 2 \omega (3\delta_k).
\end{align*}

Letting $k\to\infty$, we get the contradiction.
Observe that for any $x\in \overline{\mathcal{C}_{1/2}}$, there is some $p\in [-1,1]$ such that $(x,,t,p)\in K$. Indeed, let $x\in \overline{\mathcal{C}_{1/2}}$ and consider $(x,t,w_k(x,t))\in G_k$. By (<ref>), for all $\varepsilon>0$, there is $k_0\geq 1$ such that $(x,t,w_k(x,t)) \in B_\varepsilon(K)$ for all $k\geq k_0$. Taking $\varepsilon= 1/j$, by a diagonalization argument, we obtain that there is a subsequence of $\{(x,t,w_k(x))\}_{k\geq 1}$ converging to $(x,t,p)\in K$, where $|p|\leq 1$ by $(i)$.
We define the function $\phi: \overline{\mathcal{C}_{1/2}}\to [-1,1]$ as $\phi(x,t)=p$ if $x_n \neq 0$. To include also the case $x_n=0$ we point out that for every $\varepsilon >0$ $w_k(x',\varepsilon,t) \rightarrow \phi(x',\varepsilon,t)$ as $k \rightarrow \infty$, since $\phi$ is uniformly continuous and bounded in $\mathcal{C}_{1/2}$ we define $\phi(x',0,t):= \lim_{\varepsilon \rightarrow 0^+} \phi(x',\varepsilon,t)$.
Let $(x,t), (y,s) \in \overline{\mathcal{C}_{1/2}}$ with $(x,t)\neq (y,s)$. Then

|\phi(x,t)-\phi(y,s)| = \lim_{k\to\infty} \osc_{\substack{\mathcal{B}_{\lambda,\delta_k}(x_k,t_k)}} w_k \leq \omega(d_p((x,t),(y,s))),

where we used that $d_p((x,t),(y,s))\geq \delta_k>0$ for $k$ large enough.

By construction $\phi$ is a single-valued function with a modulus of continuity in $\overline{\mathcal{C}_{1/2}}$.

Proof of Lemma <ref>
We argue by contradiction by assuming that there exist sequences of $\delta_k \rightarrow 0$ and $w_k$ satisfying the assumption with $\delta_k$ instead of $\bar \delta$; nevertheless there exists a $\varepsilon_0 > 0$ such that for every $\phi \in C_{d_{\lambda}}^{\alpha}( \mathcal{C}_{1/2})$   it holds

\begin{equation}
\norm{w_k - \phi}_{L^{\infty}(\mathcal{C}_{1/2})} \geq \varepsilon_0.
\end{equation}

However this immediately leads to a contradiction by applying Lemma <ref>.

§ PROOF OF THE IMPROVEMENT OF FLATNESS AND MAIN RESULT

Proof of Proposition <ref>
Since the limiting profile (<ref>) coincides with the one in the homogeneous case this proof is identical to the one presented in <cit.>, having obtained the analogous results for the error function $w$ in the inhomogeneous case, see Section <ref>. For the sake of completeness we include the proof.
We proceed in two steps:

Step 1.
First we prove that there exists a function $v$ solving (<ref>) which approximates well $w$ in $\mathcal{C}_{1/2}$ that is

\begin{equation*}
\nor{v-w} \leq C \delta \quad \textit{in } \mathcal{C}_{1/2},
\end{equation*}

provided $\varepsilon \leq \varepsilon_1(\delta)$.

We want to apply Lemma <ref> to $w$ with $\bar \varepsilon=\delta$, which gives us $\bar \delta$. We thus take $\varepsilon$ so small that $\varepsilon^{1/2}\leq\bar \delta$, depending only on $\delta$ and universal constants.  By Remark <ref>, this implies that we are in the conditions of the lemma, which guarantees the existence of a (single-valued) function $\phi \in C_{d_{\lambda}}^{\alpha}( \mathcal{C}_{1/2})$ such that

\begin{equation*}
\nor{w- \phi} \leq \delta, \qquad \norm{\phi}_{C_{d_{\lambda}}^{\alpha}( \mathcal{C}_{1/2})} \leq C.
\end{equation*}

In view of Proposition <ref>, let $v$ be the solution of (<ref>) in $\mathcal{C}_{1/2}$ with $v = \phi $ on $\p_D \mathcal{C}_{1/2}$ which satisfies $\norm{v}_{C_{d_{\lambda}}^{\alpha}( \mathcal{C}_{1/2})} \leq C.$
Then if $d_\lambda((x,t), \p_D \mathcal{C}_{1/2}) \leq \delta^{1/\alpha}$ , there exists $(y,s) \in \p_D \mathcal{C}_{1/2}$ such that

\begin{align*}
\nor{v(x,t)- \phi (y,s)} &\leq C\delta,  \\
\nor{w(x,t)- \phi(y,s)} &\leq \delta + C d_\lambda((x,t), (y,s))^{\alpha} \leq C \delta,
\end{align*}

Therefore in $\mathcal{C}_{1/2} \cap \{d_\lambda((x,t)), \p_D \mathcal{C}_{1/2}) \leq \delta^{1/\alpha} \} $ we have

\[
|v-w|\leq C\delta.
\]

In particular, for $\Omega := \mathcal{C}_{1/2} \cap \{d_\lambda((x,t)), \p_D \mathcal{C}_{1/2}) > \delta^{1/\alpha} \} $, it holds

\[
|v-w|\leq C \delta, \quad \textit{ on } \,\partial_D \Omega.
\]

Using the estimates of Proposition <ref> we have $\nor{\nabla v }, \nor{D^2 v} \leq C(\delta)$ in $\Omega.$ Thus, by Corollary <ref> we have

\begin{equation*}
\nor{v- w} \leq C \delta \quad \textit{in } \Omega,
\end{equation*}

which concludes Step $1$.

Step 2. Now, applying Proposition <ref> to the solution $v$ above we obtain that

\begin{equation*}
\nor{w-\ell_{\bar a, \bar b}} \leq \nor{w-v}+ \nor{v-\ell_{\bar a, \bar b} } \leq C\rho^{1+\alpha} + C \delta \quad \textit{in } \mathcal{C}_{\rho}.
\end{equation*}

We then choose $\rho = \tau$ small, universal and $\delta= \tau^{1+\alpha/2}$ so that $\delta \leq c'$ in Proposition <ref> and

\begin{equation*}
\nor{w- \ell_{\bar a, \bar b}} \leq \frac{1}{4}\tau \quad \textit{in }\mathcal{C}_{\tau}, \quad \nor{\bar a_n'} \leq \frac{1}{4}\tau^{-2} \delta \lambda^{-1}.
\end{equation*}

Rescaling back to the original function $u$, we get

\begin{equation*}
\nor{u - \ell_{a, b} - \varepsilon \lambda \ell_{\bar a, \bar b}\left(\frac{x}{\lambda}, \frac{t}{\lambda}\right)} = \varepsilon\lambda \nor{w\left(\frac{x}{\lambda}, \frac{t}{\lambda}\right) - \ell_{\bar a, \bar b}\left(\frac{x}{\lambda}, \frac{t}{\lambda}\right)} \leq \frac{\varepsilon}{4}\tau \lambda \quad \textit{ in } \mathcal{C}_{\tau \lambda}.
\end{equation*}

Now set

\begin{equation*}
\tilde a(t):= a(t) + \varepsilon \bar a \left( \frac{t}{\lambda} \right), \quad \hat b(t) := b(t) + \varepsilon \lambda \bar b\left( \frac{t}{\lambda} \right),
\end{equation*}

for which  we can rewrite the previous inequality as

\begin{equation*}
\nor{u - \ell_{\tilde a, \hat b}} \leq \frac{\varepsilon}{4}\tau \lambda \quad \textit{ in } \mathcal{C}_{\tau \lambda}, \quad \nor{\tilde a_n'} \leq \frac{\varepsilon \delta }{\lambda^2 } \left( 1 + \frac{1}{4\tau^2}\right) \leq \frac{\varepsilon \delta}{2 (\tau \lambda)^2}.
\end{equation*}

In order to have a control on the behavior of $b$ we define $\tilde b$ as the solution to the ODE

\begin{equation*}
\tilde b' = g(\tilde a), \quad \tilde b(0) = \hat b(0),
\end{equation*}

for which we have

\[
\hat b'(t)=b'(t)+\varepsilon\bar b'\left(\frac{t}{\lambda}\right)=g(a(t))+\varepsilon\nabla g(a(t))\cdot \bar a\left(\frac{t}{\lambda}\right)=g(\tilde a(t))+O(\varepsilon^2)=\tilde b'(t)+O(\varepsilon^2).
\]

Therefore, if $t\in [-\tau \lambda,0]$ then

\[
|(\tilde b- \hat b)(t)|\leq C\varepsilon^2|t|\leq C\varepsilon^2\tau \lambda.
\]

This allows us to estimate

\begin{align*}
&\,\nor{u-\ell_{\tilde a, \tilde b}}=\nor{u-\tilde a(t)(x-\tilde b(t))}\\
\leq&\,\nor{u-\tilde a(t)(x- \hat b(t))}+\nor{\tilde a(t)(\tilde b(t) - \hat b(t))}\\
\leq&\, \frac{\varepsilon}{4}\tau \lambda+(K+C\varepsilon)(C\varepsilon^2\tau \lambda)\\
\leq&\, \frac{\varepsilon}{2}\tau \lambda
\end{align*}

for $\varepsilon\leq \varepsilon_1$ small enough, as intended. This concludes the proof.

Proof of Theorem <ref>
As already pointed out in Remark <ref>, Theorem <ref> is equivalent to obtaining $C^{1,\alpha} $ estimates up to the boundary $\{x_n = 0 \}$ for the hodograph transform. In the same remark we highlighted how the hypotheses of Theorem <ref> translate to Proposition <ref>, in particular we can take $\lambda \leq \min\{\delta \varepsilon_1, \lambda_1\}$ and we can assume that $a_0(t)=(0,\hdots, 0, (a_0)_n(t))\in \mathcal{R}_{K/2} $. We can thus iterate indefinitely Proposition <ref> in cylinders $\mathcal{C}_{\lambda \tau^k}$ with $\varepsilon=\varepsilon_k:=\varepsilon_1 2^{-k}.$ If we write $k=(\log_\tau \lambda_k - \log_\tau \lambda)$ then

\begin{equation*}
\varepsilon_1 2^{-k}=\varepsilon_1 \left( \frac{\lambda_k}{\lambda} \right)^{-\frac{1}{\log_2 \tau}} = \varepsilon_1 \lambda^\frac{1}{\log_2 \tau} \lambda_{k}^{-\frac{1}{\log_2 \tau}}=: C(\lambda) \lambda_k^{\alpha}.
\end{equation*}

Furthermore, the given coefficients $a_k(t)$ are still in $\mathcal{R}_{K/2}$, for every $k$,  since it holds

\[
\nor{a_k(t)-a_0(t)}\leq \sum_{j=1}^k\nor{a_j(t)-a_{j-1}(t)}\leq C\sum_{j=1}^k\varepsilon_k\leq C\varepsilon_1
\]

by possibly choosing $\varepsilon_1$ even smaller.

\[
|u-\ell_{a_k,b_k}|\leq \varepsilon_k\lambda_k= C(\lambda)\lambda_k^{1+\alpha} \quad \textit{in } \mathcal{C}_{\lambda_k}
\]

for all $k\geq 0$ which implies pointwise $C^{1,\alpha}(0)$  regularity in space. Moreover

\begin{align*}
\nor{a_{k+n}- a_k} &\leq \sum_{j=1}^{n} \nor{a_{k+j} - a_{k+j-1}} \leq C \lambda^{\alpha} \sum_{j=0}^{n-1} \tau^{(k+j+1) \alpha}  \\
&\leq  C \lambda^{\alpha}\tau^{(k+1) \alpha}  \sum_{j=0}^{\infty} \tau^{ j \alpha} = C \lambda_k^{\alpha} \frac{\tau^\alpha}{1-\tau^{\alpha}} \leq C \lambda_k^{\alpha},
\end{align*}

for $\tau \leq 2^{-\alpha}$.
So $a_k$ is a Cauchy sequence and for $k\rightarrow \infty$ we have $a_k \rightarrow a_\infty$ with $a_\infty= \nabla u (0,t)$.

Expanding with Taylor in space we get, for sufficiently small $r$,

\begin{align*}
\sup_{B_r^{+}}\nor{u(0,t) + \nabla u(0,t)\cdot  x + o(\nor{x}^{1+\alpha}) - a \cdot x -b} \leq C r^{1+\alpha}
\end{align*}

which implies

\begin{align*}
\sup_{B_r^{+}}\nor{ (u(0,t) - b)  + (\nabla u(0,t)- a)\cdot x +o(\nor{x}^{1+\alpha})  } \leq C r^{1+\alpha}
\end{align*}

Letting $r \rightarrow 0$ we get $u(0,t) = b$. Now, assume $\nabla u(0,t)\neq a$ and take

\[
\bar x=r\frac{\nabla u(0,t)-a}{\nor{\nabla u(0,t)-a}}\in \overline{B_r}.
\]

If $\bar x_n<0$, consider instead $-\bar x$, so that $\bar x\in \overline{B_r^+}$. Then, letting $r\to 0$ we get a contradiction.

\begin{equation} \label{eq:00}
\nor{\nabla u (0,t) - a_k(t)} \leq \sum_{j=0}^{\infty} \nor{a_j(t) - a_{j+1}(t)} \leq C \lambda_k^{\alpha}.
\end{equation}

Now from (<ref>) and (<ref>), letting $r =\nor{t-s} \in (0,2\lambda_k^2)$ for $t,s \in (t_0-\lambda_k^2, t_0+ \lambda_k^2)$

\begin{align*}
\nor{\nabla u (0,t) - \nabla u (0,s)} &\leq \sup_{t_1 \in [t-r,t+r]} \left( \nor{\nabla u (0,t) - \nabla u (0,t_1)} \right) \\
& \leq   \sup_{t_1 \in [t-r,t+r]}  \left( \nor{\nabla u (0,t) - a_k (t_0)} + \nor{\nabla u (0,t_1) - a_k (t_0)} \right) \\
&\leq 4\,r^{\alpha} = C \nor{t-s}^{\alpha}.
\end{align*}

Thus for any $t,s \in [-\lambda_k, 0]$

\begin{equation}\label{eq:01}
\begin{aligned}
\nor{a_k(t)-a_k(s)} \leq \,&\nor{a_k(t) - \nabla u (0,t)} + \nor{\nabla u (0,t) -\nabla u (0,s)} + \nor{\nabla u (0,s)- a_k(s)}\\
\leq\,&  C\lambda_k^\alpha + C\lambda_k^{\alpha/2}+C\lambda_k^\alpha  \leq 3C\lambda_k^{\alpha/2}.
\end{aligned}
\end{equation}

Being $b_k'=g(a_k) \in C^{1}$, by
(<ref>) and (<ref>), we have

\begin{equation*}
\nor{b''(t)} \leq \norm{g}_{C^1} \nor{a_k'(t)} \leq K \delta \varepsilon_k \lambda_k^{-2} = C \lambda_k^{\alpha-2}.
\end{equation*}

For $t \in (t_0 - \lambda_k^2, t_0 + \lambda_k^2)$, we get

\begin{equation} \label{eq:02}
\nor{b_k(t)- b_k(0) - b_k'(0) t} \leq  C \lambda_k^{\alpha+2}.
\end{equation}

Now in $C_{\lambda_k},$ using (<ref>) and (<ref>)

\begin{align*}
\nor{u-(a_k(0) \cdot x +b_k'(0)t + b_k(0))} &\leq \nor{u - \ell_{a_k,b_k} }  + \nor{  \ell_{a_k,b_k}  -  \ell_{a_k,b_k}(0) - b_k'(0)t} \\
& \leq  C \lambda_k^{1+\alpha} + \nor{(a_k(t)-a_k(0))\cdot x } + \nor{b_k(t)- b_k(0) - b_k'(0) t}  \\
& \leq  C \lambda_k^{1+\alpha} +  C \lambda_k^{1+\alpha/2} + C \lambda_k^{1+\alpha/2}\leq  C(\lambda) \lambda_{k}^{1+\alpha/2},
\end{align*}

which gives the desired result.

Declaration.  The authors declare that there are no relevant financial or non-financial competing interests to report.