# MINIMAL BOSONIZATION OF DOUBLE-GRADED QUANTUM MECHANICS

C. QUESNE
Physique Nucléaire Théorique et Physique Mathématique, Université Libre de
Bruxelles,
Campus de la Plaine CP229, Boulevard du Triomphe, B-1050 Brussels, Belgium
E-mail address<EMAIL_ADDRESS>
###### Abstract

The superalgebra of $\mathbb{Z}_{2}^{2}$-graded supersymmetric quantum
mechanics is shown to be realizable in terms of a single bosonic degree of
freedom. Such an approach is directly inspired by a description of the
corresponding $\mathbb{Z}_{2}$-graded superalgebra in the framework of a
Calogero-Vasiliev algebra or, more generally, of a generalized deformed
oscillator algebra. In the case of the $\mathbb{Z}_{2}^{2}$-graded
superalgebra, the central element $Z$ has the property of distinguishing
between degenerate eigenstates of the Hamiltonian.

PACS: 03.65.Fd; 11.30.Pb

Keywords: Supersymmetric quantum mechanics; $\mathbb{Z}_{2}^{2}$-graded
superalgebras; deformed oscillator algebras

##  1 Introduction

In recent years, there has been a renewed interest in the application of a
color superalgebra of
$\mathbb{Z}_{2}^{2}\equiv\mathbb{Z}_{2}\otimes\mathbb{Z}_{2}$ grading (and
more generally of $\mathbb{Z}_{2}^{n}$ grading) to quantum mechanics. Such
superalgebras were introduced in [1, 2] and extended the notion of ordinary
($\mathbb{Z}_{2}$-graded) Lie superalgebras appearing in [3, 4] (for an
introduction to color superalgebras see, e.g., [5]). Up to now, two
independent, but related, lines of research have been followed: the study of
physical models possessing a $\mathbb{Z}_{2}^{2}$-graded symmetry (see, e.g.,
[6, 7, 8, 9, 10]) and the investigation of parastatistics (see, e.g., [11, 12,
13]).

In the present paper, we are interested in an extension of
$\mathbb{Z}_{2}$-graded supersymmetric quantum mechanics (SUSYQM) [14, 15] to
$\mathbb{Z}_{2}^{2}$-graded SUSYQM, as recently proposed by Bruce and Duplij
[7] (see also [13, 16, 17]). As shown elsewhere [8], in contrast with standard
$\mathbb{Z}_{2}$-graded SUSYQM, which can be described in terms of ordinary
bosons and fermions, $\mathbb{Z}_{2}^{2}$-graded mechanics necessitates
ordinary bosons, two classes of fermions (commuting among themselves), and
exotic bosons that anticommute with the fermions of both classes.

Some years ago, it was shown that standard SUSYQM can be realized in terms of
a single bosonic degree of freedom [18]. Such a minimal bosonization
(different from the Nicolai-Witten construction in terms of two bosons [19,
20]) is obtained by using a generalized deformed oscillator algebra (GDOA)
[21, 22] and imposing a $\mathbb{Z}_{2}$-grading structure on the deformed
bosonic oscillator Fock space, which can be done, for instance, by restricting
oneself to a Calogero-Vasiliev algebra [23]. In such a framework, SUSYQM turns
out to be fully reducible, its irreducible components providing two sets of
bosonized operators realizing either broken or unbroken SUSYQM. Such an
approach in terms of GDOAs was latter on applied [24, 25, 26] to several
variants of SUSYQM, namely parasupersymmetric [27, 28, 29],
orthosupersymmetric [30], pseudosupersymmetric [31], and fractional
supersymmetric quantum mechanics [32].

The purpose of the present work is to extend the approach in terms of GDOAs to
$\mathbb{Z}_{2}^{2}$-graded SUSYQM.

##  2 $\mathbb{Z}_{2}^{2}$-graded SUSYQM

In $\mathbb{Z}_{2}^{2}$-graded SUSYQM, we are going to deal with a
$\mathbb{Z}_{2}^{2}$-graded vector space (over $\mathbb{C}$) that is the
direct sum of four homogeneous subspaces, labelled by
$\overrightarrow{a}=(a_{1},a_{2})$, where $a_{1}$, $a_{2}\in\\{0,1\\}$,

$g=g_{(0,0)}\oplus g_{(1,0)}\oplus g_{(0,1)}\oplus g_{(1,1)}.$ (1)

Homogeneous elements of $g_{\overrightarrow{a}}$ are denoted by
$X_{\overrightarrow{a}}$, $Y_{\overrightarrow{a}}$, …. The vector space $g$
becomes a $\mathbb{Z}_{2}^{2}$-graded Lie superalgebra if it admits a bilinear
form, denoted by $[[\cdot,\cdot]]$, satisfying the identities

$\displaystyle[[X_{\overrightarrow{a}},Y_{\overrightarrow{b}}]]\in
g_{\overrightarrow{a}+\overrightarrow{b}},$ (2)
$\displaystyle[[X_{\overrightarrow{a}},Y_{\overrightarrow{b}}]]=-(-1)^{\overrightarrow{a}.\overrightarrow{b}}[[Y_{\overrightarrow{b}},X_{\overrightarrow{a}}]],$
(3)
$\displaystyle(-1)^{\overrightarrow{a}.\overrightarrow{c}}[[X_{\overrightarrow{a}},[[Y_{\overrightarrow{b}},Z_{\overrightarrow{c}}]]\,]]+(-1)^{\overrightarrow{b}.\overrightarrow{a}}[[Y_{\overrightarrow{b}},[[Z_{\overrightarrow{c}},X_{\overrightarrow{a}}]]\,]]$
$\displaystyle\quad{}+(-1)^{\overrightarrow{c}.\overrightarrow{b}}[[Z_{\overrightarrow{c}},[[X_{\overrightarrow{a}},Y_{\overrightarrow{b}}]]\,]]=0,$
(4)

where

$\overrightarrow{a}+\overrightarrow{b}=(a_{1}+b_{1},a_{2}+b_{2}),\qquad\overrightarrow{a}.\overrightarrow{b}=\sum_{k=1}^{2}a_{k}b_{k}$
(5)

are calculated modulo 2. Equation (4) is called the
$\mathbb{Z}_{2}^{2}$-graded Jacobi relation. In practice, the general Lie
bracket is given by

$[[X_{\overrightarrow{a}},Y_{\overrightarrow{b}}]]=X_{\overrightarrow{a}}Y_{\overrightarrow{b}}-(-1)^{\overrightarrow{a}.\overrightarrow{b}}Y_{\overrightarrow{b}}X_{\overrightarrow{a}}$
(6)

and therefore coincides with either a commutator or an anticommutator.

In the $\mathbb{Z}_{2}^{2}$-graded version of SUSYQM, the
$\mathbb{Z}_{2}^{2}$-graded superalgebra has four generators $H_{00}$,
$Q_{10}$, $Q_{01}$, and $Z_{11}$,111From now on, for simplicity’s sake, we
skip all parentheses. which are Hermitian operators acting in a Hilbert space
${\cal H}={\cal H}_{00}\oplus{\cal H}_{10}\oplus{\cal H}_{01}\oplus{\cal
H}_{11}$ and satisfying the relations

$\\{Q_{10},Q_{10}\\}=\\{Q_{01},Q_{01}\\}=2H_{00},\qquad[Q_{10},Q_{01}]=2{\rm
i}Z_{11}.$ (7)

As a consequence of the latter, we also have

$[H_{00},Q_{10}]=[H_{00},Q_{01}]=[H_{00},Z_{11}]=0,\qquad\\{Z_{11},Q_{10}\\}=\\{Z_{11},Q_{01}\\}=0.$
(8)

The operator $H_{00}$ is the Hamiltonian, $Q_{10}$ and $Q_{01}$ are two
supercharges, and $Z_{11}$ is a central element. From now on, we will
abbreviate $H_{00}$ and $Z_{11}$ in $H$ and $Z$, respectively.

Instead of $Q_{10}$ and $Q_{01}$, we may consider their linear combinations

$Q^{\dagger}=\frac{1}{2}(Q_{10}+{\rm i}Q_{01}),\qquad
Q=\frac{1}{2}(Q_{10}-{\rm i}Q_{01}),$ (9)

from which we get

$Q_{10}=Q^{\dagger}+Q,\qquad Q_{01}=-{\rm i}(Q^{\dagger}-Q).$ (10)

On inserting these relations in Eq. (7), the latter gives rise to the
equivalent relations

$\displaystyle\\{Q^{\dagger},Q\\}=H,$ (11)
$\displaystyle(Q^{\dagger})^{2}+Q^{2}=0,$ (12)
$\displaystyle[Q^{\dagger},Q]=Z.$ (13)

Furthermore, Eq. (8) becomes

$[H,Q^{\dagger}]=[H,Q]=[H,Z]=0,\qquad\\{Z,Q^{\dagger}\\}=\\{Z,Q\\}=0.$ (14)

This reformulation of $\mathbb{Z}_{2}^{2}$-graded SUSYQM is closer to a usual
presentation of standard SUSYQM and will enable us to directly apply some
results known for the latter. In the next section, we will proceed to consider
the simplest approach of SUSYQM in terms of the Calogero-Vasiliev algebra,
leaving another one in terms of a more general GDOA for Sec. 4.

##  3 Bosonization in Terms of a Calogero-Vasiliev Algebra

Let us first recall that the SUSYQM superalgebra, generated by $H$,
$Q^{\dagger}$, and $Q$, such that

$\\{Q^{\dagger},Q\\}=H,\qquad(Q^{\dagger})^{2}=Q^{2}=0,\qquad[H,Q^{\dagger}]=[H,Q]=0,$
(15)

can be realized in terms of the Calogero-Vasiliev algebra, defined by [23]

$[N,a^{\dagger}]=a^{\dagger},\qquad[N,a]=-a,\qquad[a,a^{\dagger}]=1+\kappa T,$
(16)

where $N=N^{\dagger}$, $a=(a^{\dagger})^{\dagger}$, $T=e^{{\rm i}\pi
N}=(-1)^{N}$, and $\kappa$ is some real deformation parameter. The operator
$T$ is known as the Klein operator. On assuming the existence of a Fock-type
vacuum state $|0\rangle$, such that

$N|0\rangle=a|0\rangle=0,\qquad T|0\rangle=|0\rangle,$ (17)

the Calogero-Vasiliev algebra admits a Fock space representation, whose
carrier space is spanned by the orthonormalized vectors

$|n\rangle=\frac{1}{\sqrt{[n]_{\kappa}!}}(a^{\dagger})^{n}|0\rangle,\qquad\langle
n^{\prime}|n\rangle=\delta_{n^{\prime},n},$ (18)

where $[n]_{\kappa}!=\prod_{i=1}^{n}[i]_{\kappa}$,
$[i]_{\kappa}=i+\frac{\kappa}{2}\bigl{(}1-(-1)^{i}\bigr{)}$. These vectors are
such that $N|n\rangle=n|n\rangle$,
$a|n\rangle=\sqrt{[n]_{\kappa}}|n-1\rangle$, and
$a^{\dagger}|n\rangle=\sqrt{[n+1]_{\kappa}}|n+1\rangle$. One may write

$a^{\dagger}a=N+\kappa P_{1},\qquad aa^{\dagger}=N+1+\kappa P_{0},\qquad
N=\frac{1}{2}\\{a^{\dagger},a\\}-\frac{1}{2}(\kappa+1),$ (19)

where $P_{0}=\frac{1}{2}(1+T)$ and $P_{1}=\frac{1}{2}(1-T)$ are the projectors
on even and odd number states, respectively.

The superalgebra (15) may be realized in terms of the Calogero-Vasiliev
algebra in two different ways [18], namely

$Q^{\dagger}_{0}=aP_{0},\qquad Q_{0}=a^{\dagger}P_{1},\qquad
H_{0}=\frac{1}{2}\\{a^{\dagger},a\\}-\frac{1}{2}T[a,a^{\dagger}]=N+\frac{1}{2}-\frac{1}{2}T,$
(20)

and

$Q^{\dagger}_{1}=aP_{1},\qquad Q_{1}=a^{\dagger}P_{0},\qquad
H_{1}=\frac{1}{2}\\{a^{\dagger},a\\}+\frac{1}{2}T[a,a^{\dagger}]=N+\kappa+\frac{1}{2}+\frac{1}{2}T.$
(21)

The corresponding (linear) spectrum of $H_{0}$ and $H_{1}$ is given by

${\cal E}^{(0)}_{n}=2\left[\frac{n+1}{2}\right]_{\kappa},\qquad
n=0,1,2,\ldots,$ (22)

and

${\cal E}^{(1)}_{n}=2\left[\frac{n}{2}\right]_{\kappa}+1+\kappa,\qquad
n=0,1,2,\ldots,$ (23)

respectively. This means that ${\cal E}^{(0)}_{0}=0$, ${\cal
E}^{(0)}_{2k+1}={\cal E}^{(0)}_{2k+2}=2k+2$, $k=0$, 1, 2, …, whereas ${\cal
E}^{(1)}_{2k}={\cal E}^{(1)}_{2k+1}=2k+1+\kappa$, $k=0$, 1, 2, …. The former
case corresponds to unbroken SUSYQM, while the latter is associated with
broken SUSYQM. This has led to the conclusion [18] that in the framework of
Calogero-Vasiliev algebra, SUSYQM is both minimally bosonized and fully
reducible.

On comparing Eq. (15) with Eqs. (11) and (12), it is obvious that the sets
$Q^{\dagger}_{0}$, $Q_{0}$, $H_{0}$ and $Q^{\dagger}_{1}$, $Q_{1}$, $H_{1}$ of
Eqs. (20) and (21) also provide a realization of the
$\mathbb{Z}_{2}^{2}$-graded superalgebra. It only remains to use Eq. (13) to
determine the corresponding $Z$ operator. The results read

$Z_{0}=-TH_{0}\qquad\text{and}\qquad Z_{1}=TH_{1},$ (24)

respectively. As a consequence, in both scenarii, the central element
distinguishes between the degenerate eigenstates of $H$ in the Fock space
representation, since its eigenvalues are given by

${\cal Z}^{(0)}_{0}=0,\qquad{\cal Z}^{(0)}_{2k+1}=-{\cal
Z}^{(0)}_{2k+2}=2k+2,\qquad k=0,1,2,\ldots,$ (25)

and

${\cal Z}^{(1)}_{2k}=-{\cal Z}^{(1)}_{2k+1}=2k+1+\kappa,\qquad
k=0,1,2,\ldots,$ (26)

respectively.

##  4 Bosonization in Terms of a GDOA

A more general realization of the SUSYQM superalgebra (15) can be obtained by
replacing the Calogero-Vasiliev algebra by a GDOA, defined by [21, 22]

$[N,a^{\dagger}]=a^{\dagger},\qquad[N,a]=-a,\qquad[a,a^{\dagger}]=G(N),$ (27)

where $N=N^{\dagger}$, $a=(a^{\dagger})^{\dagger}$, and
$G(N)=[G(N)]^{\dagger}$ is some Hermitian function of $N$. We restrict
ourselves to GDOAs possessing a bosonic Fock space representation. In the
latter, we may write $a^{\dagger}a=F(N)$, $aa^{\dagger}=F(N+1)$, where the
structure function $F(N)=[F(N)]^{\dagger}$ is such that

$G(N)=F(N+1)-F(N)$ (28)

and is assumed to satisfy the conditions

$F(0)=0,\qquad F(n)>0\qquad\text{if $n=1,2,3,\ldots.$}$ (29)

The carrier space ${\cal F}$ of such a representation can be constructed from
a vacuum state $|0\rangle$ (such that $a|0\rangle=N|0\rangle=0$) by successive
applications of $a^{\dagger}$ and its orthonormalized basis states

$|n\rangle=\left(\prod_{i=1}^{n}F(i)\right)^{-1/2}(a^{\dagger})^{n}|0\rangle,\qquad
n=0,1,2,\ldots,$ (30)

satisfy the relations $N|n\rangle=n|n\rangle$,
$a^{\dagger}|n\rangle=\sqrt{F(n+1)}\,|n+1\rangle$, and
$a|n\rangle=\sqrt{F(n)}\,|n-1\rangle$.

A $\mathbb{Z}_{2}$-grading structure can be imposed on $\cal F$ by introducing
a grading operator $T=e^{{\rm i}\pi N}=(-1)^{N}$ with the two eigenvalues
$(-1)^{\mu}$, $\mu=0$, 1. The two corresponding eigenspaces ${\cal F}_{0}$ and
${\cal F}_{1}$, made of even or odd number states, respectively, are such that
${\cal F}={\cal F}_{0}\oplus{\cal F}_{1}$, the corresponding projectors being
$P_{0}=\frac{1}{2}(1+T)$ and $P_{1}=\frac{1}{2}(1-T)$.

On starting from the general results obtained for order-$p$ parasupersymmetric
quantum mechanics in [24] and setting $p=1$, which gives back SUSYQM, we
arrive at two different realizations of superalgebra (15) in terms of the GDOA
(27), namely

$Q^{\dagger}_{0}=f(N)a^{\dagger}P_{1},\quad Q_{0}=f(N+1)aP_{0},\quad
H_{0}=f^{2}(N)F(N)P_{0}+f^{2}(N+1)F(N+1)P_{1},$ (31)

and

$Q^{\dagger}_{1}=f(N)a^{\dagger}P_{0},\quad Q_{1}=f(N+1)aP_{1},\quad
H_{0}=f^{2}(N)F(N)P_{1}+f^{2}(N+1)F(N+1)P_{0}.$ (32)

Here $f(N)$ is an arbitrary real function of $N$ defined on the set
$\\{1,2,3,\ldots\\}$. The eigenvalues of $H_{0}$ and $H_{1}$ are given by

${\cal E}^{(0)}_{0}=0,\qquad{\cal E}^{(0)}_{2k+1}={\cal
E}^{(0)}_{2k+2}=f^{2}(2k+2)F(2k+2),\qquad k=0,1,2,\ldots,$ (33)

and

${\cal E}^{(1)}_{2k}={\cal E}^{(1)}_{2k+1}=f^{2}(2k+1)F(2k+1),\qquad
k=0,1,2,\ldots,$ (34)

respectively. and therefore correspond to either unbroken or broken SUSYQM
again. This time, however, the spectra are in general nonlinear.

As already observed in Sec. 3, the operators (31) and (32) also provide
realizations of the $\mathbb{Z}_{2}^{2}$-graded superalgebra equations (11)
and (12). The calculation of $Z$, given in Eq. (13), now leads to

$Z_{0}=-f^{2}(N+1)F(N+1)P_{1}+f^{2}(N)F(N)P_{0}$ (35)

and

$Z_{1}=-f^{2}(N+1)F(N+1)P_{0}+f^{2}(N)F(N)P_{1}.$ (36)

By comparing Eqs. (35) and (36) with (31 and (32), it is then straightforward
to show that

$Z_{\mu}=(-1)^{\mu}TH_{\mu},\qquad\mu=0,1.$ (37)

Hence, the eigenvalues of $Z_{0}$ and $Z_{1}$ are given by

${\cal Z}^{(0)}_{0}=0,\qquad{\cal Z}^{(0)}_{2k+1}=-{\cal
Z}^{(0)}_{2k+2}=-f^{2}(2k+2)F(2k+2),\qquad k=0,1,2,\ldots,$ (38)

and

${\cal Z}^{(1)}_{2k}=-{\cal Z}^{(1)}_{2k+1}=-f^{2}(2k+1)F(2k+1),\qquad
k=0,1,2,\ldots,$ (39)

respectively. We conclude that the central element distinguishes between
degenerate eigenstates of $H$ again.

Whenever the GDOA reduces to the Calogero-Vasiliev algebra, i.e.,
$G(N)=1+\kappa T$, Eq. (19) leads to $F(N)=N+\kappa P_{1}$ and
$F(N+1)=N+1+\kappa P_{0}$. On choosing in addition $f(N)=1$, Eqs. (31), (32),
and (37) give back Eqs. (20), (21), and (24) with $Q\leftrightarrow
Q^{\dagger}$ and $Z\leftrightarrow-Z$. Such discrepancies are of course
irrelevant.

##  5 Conclusion

In this work, we have shown that the minimal bosonization of the
$\mathbb{Z}_{2}$-graded SUSYQM superalgebra in terms of the Calogero-Vasiliev
algebra or, more generally, of a GDOA can be extended to the
$\mathbb{Z}_{2}^{2}$-graded one in a straightforward way. For the latter, the
central element $Z$ has the property of distinguishing between degenerate
eigenstates of the Hamiltonian.

An interesting open question for future investigation would be the search for
a bosonization of the $\mathbb{Z}_{2}^{3}$-graded SUSYQM superalgebra in the
framework of GDOAs.

## Acknowledgements

This work was supported by the Fonds de la Recherche Scientifique–FNRS under
Grant No. 4.45.10.08.

## References

  * [1] V. Rittenberg and D. Wyler, Nucl. Phys. B139, 189 (1978).
  * [2] V. Rittenberg and D. Wyler, J. Math. Phys. 19, 2193 (1978).
  * [3] V. G. Kac, Adv. Math. 26, 8 (1977).
  * [4] M. Scheunert, J. Math. Phys. 20, 712 (1979).
  * [5] N. Aizawa, Adv. Appl. Clifford Algebras 28, 28 (2018).
  * [6] N. Aizawa, Z. Kuznetsova, H. Tanaka and F. Toppan, Prog. Theor. Exp. Phys. 123A01 (2016).
  * [7] A. J. Bruce and S. Duplij, J. Math. Phys. 61, 063503 (2020).
  * [8] N. Aizawa, Z. Kuznetsova and F. Toppan, Eur. Phys. J. C80, 668 (2020).
  * [9] N. Aizawa, Z. Kuznetsova and F. Toppan, Nucl. Phys. B967, 115426 (2021).
  * [10] A. J. Bruce, J. Phys. A53, 455201 (2020).
  * [11] V. N. Tolstoy, Phys. Part. Nuclei Lett. 11, 933 (2014).
  * [12] N. I. Stoilova and J. Van der Jeugt, J. Phys. A51, 135201 (2018).
  * [13] F. Toppan, J. Phys. A54, 115203 (2021).
  * [14] F. Cooper, A. Khare and U. Sukhatme, Phys. Rep. 251, 267 (1995).
  * [15] G. Junker, Supersymmetric Methods in Quantum and Statistical Physics (Springer, Berlin, 1996).
  * [16] A. J. Bruce, Symmetry 11, 116 (2019).
  * [17] N. Aizawa, K. Amakawa and S. Doi, J. Math. Phys. 61, 052105 (2020).
  * [18] M. S. Plyushchay, Ann. Phys. (N. Y.) 245, 339 (1998).
  * [19] E. Witten, Nucl. Phys. B188, 513 (1981).
  * [20] H. Nicolai, J. Phys. A9, 1497 (1976).
  * [21] C. Daskaloyannis, J. Phys. A24, L789 (1991).
  * [22] C. Quesne and N. Vansteenkiste, Helv. Phys. Acta 69, 141 (1996).
  * [23] M. A. Vasiliev, Int. J. Mod. Phys. A6, 1115 (1991).
  * [24] C. Quesne and N. Vansteenkiste, Mod. Phys. Lett. A17, 839 (2002).
  * [25] C. Quesne and N. Vansteenkiste, Int. J. Mod. Phys. A18, 271 (2003).
  * [26] C. Quesne, Mod. Phys. Lett. A18, 515 (2003).
  * [27] V. A. Rubakov and V. P. Spiridonov, Mod. Phys. Lett. A3, 1337 (1988).
  * [28] A. Khare, J. Math. Phys. 34, 1277 (1993).
  * [29] J. Beckers and N. Debergh, Nucl. Phys. B340, 767 (1990).
  * [30] A. Khare, A. K. Mishra and G. Rajasekaran, Int. J. Mod. Phys. A8, 1245 (1993).
  * [31] J. Beckers and N. Debergh, Int. J. Mod. Phys. A10, 2783 (1995).
  * [32] S. Durand, Phys. Lett. B312, 115 (1993).


# Equivalent formulations of the oxygen depletion problem, other implicit free
boundary value problems, and implications for numerical approximation

Xinyu Cheng111School of mathematical sciences, Fudan University ✉
<EMAIL_ADDRESS>Zhaohui Fu222University of British Columbia, ✉
<EMAIL_ADDRESS>Brian Wetton333University of British Columbia, ✉
<EMAIL_ADDRESS>
###### Abstract

The Oxygen Depletion problem is an implicit free boundary value problem. The
dynamics allow topological changes in the free boundary. We show several
mathematical formulations of this model from the literature and give a new
formulation based on a gradient flow with constraint. All formulations are
shown to be equivalent. We explore the possibilities for the numerical
approximation of the problem that arise from the different formulations. We
show a convergence result for an approximation based on the gradient flow with
constraint formulation that applies to the general dynamics including
topological changes. More general (vector, higher order) implicit free
boundary value problems are discussed. Several open problems are described.

##  1 Introduction

The Oxygen Depletion (OD) problem is a free boundary value problem of implicit
type. Implicit here means that the free boundary is specified implicitly by an
extra boundary condition rather than explicitly as an interface normal
velocity as for a Stefan problem [26, 30, 22].

The OD problem was introduced as a model of oxygen consumption and diffusion
in living tissue but several other problems have similar structure. Some of
the early work is described in [6] with a great deal of subsequent interest
from the analysis and numerical research communities in [26, 24, 8, 2, 20].
Reference [20] has a review of much of the previous work. In the current work,
we pursue an understanding of the analysis of the OD problem as the simplest
example of an implicit moving free boundary value problem. We are motivated by
an interest in the analysis and computation of a general class of implicit
free boundary value problems.

By way of introduction, we present the OD problem in 1D for an unknown
$u(x,t)$ for $x\in[0,s(t)]$ with a single free boundary $x=s(t)$ and a no flux
condition $u_{x}=0$ at $x=0$. At the free boundary, $u=0$ and additionally
$u_{x}=0$. These two conditions implicitly define the free boundary $x=s(t)$.
The solution obeys

$u_{t}=u_{xx}-1$ (1.1)

for $x\in[0,s(t)]$ and it is natural to extend $u\equiv 0$ for $x>s(t)$ in a
$C^{1}$ continuous way. We consider positive initial conditions for $u$ in
$[0,s(0))$. This is one of the forms of the OD problem foreshadowed by the
title. We consider several formulations in the literature and in new results
we show they are all equivalent. We introduce a new formulation as the $L_{2}$
gradient flow with constraint on the energy from the elliptic obstacle
problem. The obstacle problem has had considerable interest in the literature
[10, 5, 31, 16, 21]. Some discussion of the numerical methods that follow from
the different formulations is given.

######  Remark 1.1.

We highlight that the extra boundary condition for implicit free boundary
value problems does not explicitly contain the interface velocity, hence this
velocity is determined implicitly. For steady state free boundary problems,
where the interface velocity is zero, the difference between implicit and
explicit formulations disappears.

######  Remark 1.2.

The Oxygen Depletion problem is also known as the Oxygen Diffusion problem in
the literature. We prefer the former name as it is the depletion, the second
term on the RHS of (1.1), rather than the diffusion, the first term on the RHS
of (1.1), that leads to the formation of free boundaries.

We invite the reader to view computational examples of the dynamics in Figures
2, 3, and 4. Solutions of (1.1) can go negative, but physically relevant
values of concentration have $u\geq 0$. In the 1D case, preserving
nonnegativity results in the break up or merger of intervals where $u>0$ as
shown in Figure 3. Topological change can be more complex in higher dimensions
as seen in Figure 4. Some of the problem formulations we consider can handle
topological changes while others cannot.

The paper is organized as follows. In Section 2 we present the different
formulations and show their equivalence. In Section 3 we present two numerical
schemes. One scheme gives high accuracy solutions to the 1D problem without
topological change. The other scheme, based on our new gradient formulation of
the problem, can be applied in higher dimensions and can handle topological
changes. A convergence proof for this new scheme is given. In Section 4 we
present some other implicit free boundary value problems of interest and
indicate how our results can be extended to them, with some open questions.
The analysis of a biharmonic problem with gradient flow structure follows
directly from our new formulation of the OD problem. We end with a short
Summary that includes a list of open problems.

## Notation

We define the space $H^{1}_{+}(\Omega)\coloneqq\\{u\in H^{1}(\Omega):u\geq 0\
a.e.,\frac{\partial u}{\partial n}|_{\partial\Omega}=0\\}$. For simplicity we
consider $d-$dimensional open connected bounded domain $\Omega=(0,1)^{d}$ with
homogeneous Neumann boundary conditions, where $d=1,2,3.$ We further denote
$\mathcal{J}$ to be collection of functions $v\in L^{2}(0,T;H^{1}(\Omega))$
such that $v(t)\in H^{1}_{+}(\Omega)$ for a.e. $t\in(0,T)$. In some instances,
we denote the time derivative by $\dot{u}(t)$ and the space derivative in 1D
case by $u^{\prime}(x)$. s Given two quantities $A$ and $B$, we use $A\lesssim
B$ to denote that there exists a constant $C>0$ such that $A\leq C\cdot B$.

##  2 Equivalent Formulations

###  2.1 1D formulations without topological change

####  2.1.1 Standard formulation in 1D

The one-dimensional oxygen depletion problem with associated free boundary and
initial conditions is as follows:

$\left\\{\begin{aligned} &u_{t}=u_{xx}-1,&&\ 0\leq x\leq s(t)\\\ &u(x,t)=0,&&\
x>s(t)\\\ &u_{x}(0,t)=0,&&\ t>0\\\ &u(s(t),t)=u_{x}(s(t),t)=0,&&\ t>0\\\
&u(x,0)=u_{0}(x),&&\ 0\leq x\leq 1\\\ &s(0)=1.&&\end{aligned}\right.$ (2.1)

We assume here that $u_{0}$ satisfies all necessary smoothness and
compatibility assumptions needed in the analyses cited below. By literature
convention, we consider here a problem with a fixed, no-flux boundary
condition at $x=0$ and only one free boundary $s(t)>0$. Uniqueness and lack of
topological change when $u_{0}^{\prime}\leq 0$ follows from a modified maximum
principle argument [8].

Existence can be seen by considering $v=u_{t}$ which satisfies a standard
Stefan problem [6] with explicit interface velocity:

$\left\\{\begin{aligned} &v_{t}=v_{xx},&&\ 0\leq x\leq s(t)\\\
&v_{x}(0,t)=0,&&\ t>0\\\ &v(s(t),t)=0,\ v_{x}(s(t),t)=-\dot{s}(t),&&\ t>0\\\
&s(0)=1.&&\end{aligned}\right.$

One can check the function $u=\int_{0}^{t}v\ d\tau$ solves the oxygen
depletion problem. To prove existence and uniqueness of Stefan problem, one
can verify that the map

$\mathcal{T}(s)(t)\coloneqq 1-\int_{0}^{t}v_{x}(s(\tau),\tau))\ d\tau,\ T\geq
t\geq 0,$

defines a contraction map [18].

######  Remark 2.1.

The reformulation in $v=u_{t}$ to an explicit free boundary problem with
interface velocity equal to $-v_{x}$ can be reinterpreted as a normal velocity
for the problem for $u$ with velocity equal to $-v_{x}=-u_{tx}=-u_{xxx}$. The
authors are not aware of any analysis or computational methods based on this
velocity expression with higher order spatial derivatives.

####  2.1.2 Mapped domain formulation in 1D

Considering the same smooth solutions without topological change in 1D
discussed above, we consider $s(t)>0$ in $t\in[0,T]$, take $y=x/s(t)$, and
reformulate the oxygen depletion problem as

$u_{yy}+\dot{s}syu_{y}-s^{2}u_{t}-s^{2}=0$ (2.2)

with boundary conditions $u_{y}(0,t)=u(1,t)=u_{y}(1,t)=0$. Over a short time
period, we assume that $\dot{s}(t)$ and $s(t)$ are uniformly bounded, thus the
linear operator is parabolic. Assuming $s(t)$ is known, uniqueness of $u$ is
not an issue; however, to prove uniqueness of the solution pair
$(\tilde{u},s)$, we introduce the map $\mathcal{G}$: $X\to Y$, where $X$ is
the closed subspace of $H^{1}(H^{2}([0,s(t)]);[0,T])\times C^{1}([0,T])$ that
solves OD system and $Y$ is the closed subspace of
$H^{1}(H^{2}([0,1]);[0,T])\times C^{1}([0,T])$ that solves the reformulated
system:

$\displaystyle\mathcal{G}((u(x,t),s(t))=(u(y,t),s(t))\coloneqq(U(y\cdot
s(t),t),s(t)).$

where $U(x,t)$ is the solution from the previous section. One can check that
the map $\mathcal{G}$ is a bijection and so all solutions of (2.2) are
equivalent to the solutions in the standard formulation of Section 2.1.1.

A numerical method based on this formulation is presented in Section 3.1. The
computations in Figure 2 are done with a method based on this formulation.

######  Remark 2.2.

A direct analysis of this formulation would be useful as a stepping stone to a
convergence proof for the numerical approximation in Section 3.1 and an
analysis of the general class of problems in Section 4. We have not been able
to make progress on such an analysis. There are subtleties in the problem:
note that changing $-s^{2}$ to $+s^{2}$ makes the problem ill defined as
$s(t)=+\infty$ for $t>0$ in that case.

###  2.2 Higher dimensional formulations that allow topological change

A weak form of the solution can be introduced using a variational inequality
approach (2.3)[19, 23]. This is described in Section 2.2.1 below. We use this
formulation as the basis for equivalence to the others. This formulation is
amenable to approximation using the Augmented Lagragian Method [14, 13]. We
then introduce a new formulation as $L^{2}$ gradient flow with constraint on
the energy from the elliptic obstacle problem in Section 2.2.2. The
computations in Figures 3 and 4 are done with a method based on this
formulation. We show a regularized approach with parameter $\epsilon$, similar
to the approach in [2], in Section 2.2.3.

####  2.2.1 A parabolic variational inequality formulation

To proceed with the discussion of the problem in higher dimensions with
topological changes, we consider the standard approach to weak solutions in
this setting: a variational inequality formulation [17, 19]. This approach has
been well studied and we describe results in the literature.

We consider the following problem: find a function $u\in\mathcal{J}$ with
$u(0)=u_{0}\in H_{+}^{1}(\Omega)$ that solves

$\displaystyle\int_{0}^{t}\int_{\Omega}u_{t}\cdot(v-u)+\int_{0}^{t}\int_{\Omega}\nabla
u\cdot\nabla(v-u)\geq\int_{0}^{t}\int_{\Omega}u-v;\ \mbox{for all}\
v\in\mathcal{J},\ \mbox{a.e.}\ t\in(0,T).$ (2.3)

######  Proposition 2.3.

The variational inequality (2.3) has at most one solution and in fact suppose
$u_{1}$ and $u_{2}$ solves (2.3) with distinct initial conditions $u_{1_{0}}$
and $u_{2_{0}}$ then

$\|u_{1}-u_{2}\|_{L^{\infty}(0,T;L^{2}(\Omega))}\leq\|u_{1_{0}}-u_{2_{0}}\|_{L^{2}(\Omega)}.$
(2.4)

###### Proof.

Note $u_{j}\in\mathcal{J}$ satisfies (2.3) for j=1,2, in particular

$\displaystyle\int_{0}^{t}\int_{\Omega}\partial_{t}u_{1}\cdot(u_{2}-u_{1})+\int_{0}^{t}\int_{\Omega}\nabla
u_{1}\cdot\nabla(u_{2}-u_{1})\geq\int_{0}^{t}\int_{\Omega}u_{1}-u_{2},$
$\displaystyle\int_{0}^{t}\int_{\Omega}\partial_{t}u_{2}\cdot(u_{1}-u_{2})+\int_{0}^{t}\int_{\Omega}\nabla
u_{2}\cdot\nabla(u_{1}-u_{2})\geq\int_{0}^{t}\int_{\Omega}u_{2}-u_{1}.$

Summing the two inequalities above and denote $w=u_{1}-u_{2}$, one has

$\displaystyle\int_{0}^{t}\int_{\Omega}\partial_{t}w\cdot
w+\int_{0}^{t}\int_{\Omega}\nabla w\cdot\nabla w\leq 0$
$\displaystyle\implies$ $\displaystyle\int_{0}^{t}\int_{\Omega}(w^{2})_{t}\leq
0\implies\|w\|_{L^{\infty}(0,T;L^{2}(\Omega))}\leq\|w_{0}\|_{L^{2}(\Omega)}.$

∎

######  Theorem 2.4.

There exists a unique solution to the variational inequality (2.3).

Note that this can be done by a standard monotone operator argument and we
refer to [19].

We show equivalance to the 1D formulations. Any smooth solution $u$ to (2.1)
must solve (2.3) in the 1D case and by uniqueness the solution to (2.3)
therefore solves OD. To see this, we first observe that $u\geq 0$ and
therefore $u\in\mathcal{J}$. Resulting from that, for any $v\in\mathcal{J}$
and for a.e. $t\in(0,T)$ by applying integration by parts we obtain that

$\displaystyle\int_{0}^{t}\int_{0}^{1}u_{t}\cdot(v-u)+\int_{0}^{t}\int_{0}^{1}u_{x}\cdot(v_{x}-u_{x})=$
$\displaystyle\int_{0}^{t}\int_{0}^{s(\tau)}(u_{t}-u_{xx})(v-u)\ dxd\tau$
$\displaystyle=$ $\displaystyle\int_{0}^{t}\int_{0}^{s(\tau)}u-v\
dxd\tau\geq\int_{0}^{t}\int_{\Omega}u-v.$

####  2.2.2 A gradient flow formulation

In this section, we formulate the OD problem as the $L^{2}$ gradient of the
energy from the elliptic obstacle problem. A formal calculation with

$\mathcal{E}(t)\coloneqq\int_{\Omega}\frac{1}{2}|\nabla u|^{2}+u$

leads to

$\frac{d\mathcal{E}}{dt}=-\int_{\Omega}(\Delta u-1)^{2}.$

It is convenient to present the equivalence of the gradient flow formulation
as the limit of implicit time steps as this gets us half way to the
convergence result for the fully discrete method described in Section 3.2. The
spatially continuous, time discrete solutions $u_{n}$ approximate
$u(\cdot,nk)$, where $k$ is a time step. We consider the following
minimization problem for $u=u_{n+1}$ to the following energy functional:

$E[u]=\int_{\Omega}\frac{1}{2}|\nabla u|^{2}+\frac{1}{2k}(u-u_{n})^{2}+u,$
(2.5)

where $u\in H^{1}_{+}(\Omega)$. Existence and uniqueness of the minimizer is
guaranteed by the standard calculus of variation technique and convexity of
the energy functional [27].

######  Remark 2.5.

By defining the discrete energy
$\mathcal{E}^{n}\coloneqq\int_{\Omega}\frac{1}{2}|\nabla u_{n}|^{2}+u_{n},$ we
can see that $\mathcal{E}^{n+1}\leq\mathcal{E}^{n}.$ This can be derived by
considering $E[u]=\int_{\Omega}\frac{1}{2}|\nabla
u|^{2}+u+\frac{(u-u_{n})^{2}}{2k}$ with $E[u_{n+1}]\leq E[u_{n}]$. This gives
the discrete gradient flow structure.

We will derive the corresponding Euler-Lagrange equation for the minimizing
problem following the idea from [10]. We give an adapted proof in our case for
completeness.

######  Theorem 2.6.

Suppose $u_{n+1}$ is the unique minimizer to the energy minimizing problem
(2.5), then $u_{n+1}$ is the (weak) solution to the following modified
backward Euler scheme:

$\frac{u_{n+1}-u_{n}\cdot\chi_{\\{u_{n+1}>0\\}}}{k}=\Delta
u_{n+1}-\chi_{\\{u_{n+1}>0\\}}.$

To begin with, we consider an equivalent energy minimizing problem:

$\tilde{E}[u]\coloneqq\int_{\Omega}\frac{1}{2}|\nabla
u|^{2}+\frac{1}{2k}u^{2}+(1-\frac{u_{n}}{k})u^{+}+\frac{1}{2k}u_{n}^{2},$
(2.6)

subject to

$u\in\tilde{\mathcal{K}}\coloneqq\\{v\in H^{1}(\Omega):\frac{\partial
v}{\partial n}|_{\partial\Omega}=0\\}$

where $u^{+}=\max(u,0)$.

######  Lemma 2.7.

There exists a unique $\tilde{u}\in\tilde{\mathcal{K}}$ such that

$\tilde{E}[\tilde{u}]=\min_{v\in\tilde{\mathcal{K}}}\tilde{E}[v];$

moreover such $\tilde{u}$ is the unique minimizer to (2.5).

###### Proof.

Firstly, by similar argument, the existence and uniqueness of this energy
minimizing problem can be proved.

Now to show the equivalence of these two minimizing problems, we recall that
the minimizer $u\geq 0$, so

$\min_{v\in
H^{1}_{+}(\Omega)}E[v]=E[u]=\tilde{E}[u]\geq\min_{v\in\tilde{\mathcal{K}}}\tilde{E}[v];$

On the other hand, in order to show

$\min_{v\in\mathcal{K}}E[v]\leq\min_{v\in\tilde{\mathcal{K}}}\tilde{E}[v],$

we note that for any $v\in\tilde{\mathcal{K}}$, the corresponding $v^{+}\in
H^{1}_{+}(\Omega)$. As a result,

$E[u]\leq E[v^{+}]\leq\tilde{E}[v]$

for any $v\in\tilde{\mathcal{K}}$, therefore we get

$E[u]\leq\min_{v}\tilde{E}[v].$

Now since $E[u]=\tilde{E}[u]=\min\tilde{E}[v]$, we have $\tilde{u}=u$ by the
uniqueness.

∎

It remains to derive the Euler-Lagrange equation for this new energy
minimizing scheme.

######  Proposition 2.8.

Suppose $u$ is the unique solution to the minimization problem (2.6), then $u$
is the (weak) solution to the following modified backward Euler scheme:

$\frac{u-u_{n}\cdot\chi_{\\{u>0\\}}}{k}=\Delta u-\chi_{\\{u>0\\}}.$

The proof of this result is found in Appendix A.

We follow the idea in [3] to formulate the minimization problem (2.5) as a
variational inequality:

$u\in H^{1}_{+}(\Omega):\int_{\Omega}\nabla
u\cdot\nabla(v-u)+\frac{u}{k}(v-u)\
dx\geq\int_{\Omega}\left(\frac{u_{n}}{k}-1\right)\cdot(v-u)\mbox{for all }v\in
H^{1}_{+}(\Omega).$ (2.7)

To see the equivalence of the energy minimization and elliptic variational
inequality we now state the proposition.

######  Proposition 2.9.

Any solution to the minimization problem (2.5) is also a solution to the
variational inequality (2.7) and vice versa.

###### Proof.

Suppose $u$ is an energy minimizer to (2.5). Let $v\in H^{1}_{+}(\Omega)$,
note that $H^{1}_{+}(\Omega)$ is convex then $(1-\lambda)u+\lambda v\in
H^{1}_{+}(\Omega)$ for any $\lambda\in[0,1]$. Using $(1-\lambda)u+\lambda v$
as a competitor in $E[u]\leq E[(1-\lambda)u+\lambda v]$, we can derive from
the order $O(\lambda)$:

$\displaystyle\int_{\Omega}\nabla u\cdot\nabla(v-u)+\frac{u}{k}(v-u)\
dx\geq\int_{\Omega}\left(\frac{u_{n}}{k}-1\right)\cdot(v-u),\ \forall\ v\in
H^{1}_{+}(\Omega).$

The reverse can be proved similarly.

∎

Note that this formulation uses convexity of $H^{1}_{+}(\Omega)$; the optimal
regularity of $u$ is $C^{1,1}_{loc}$:

######  Theorem 2.10 (regularity).

Suppose u is a solution to (2.5) (or (2.7)), then there exists a positive
constant $C$ such that

$\|\Delta u\|_{\infty}\leq C\Big{(}1+\frac{1}{k}\|u_{n}\|_{\infty}+\|\Delta
u_{n}\|_{\infty}\Big{)}.$

Moreover, for each compact $K\subset\Omega$ there exists a positive constant
$c(K)>0$ such that

$\sup_{i,j}\sup_{x\in K}|D_{ij}u(x)|\leq c.$

The proof follows from [3], where penalty argument is applied together with a
non-degeneracy argument, which we refer to Lemma 1.2 from [3].

######  Remark 2.11.

This upper bound can be improved by applying energy gradient flow. By
competing $u$ with $u_{n}$ in $E[u]\leq E[u_{n}]$, we have

$\frac{1}{k}\|u-u_{n}\|_{2}^{2}\leq\|\nabla u_{n}\|_{2}^{2}+\|u_{n}\|_{1}.$

By applying the penalty argument in [3], we can derive that

$\|\Delta u\|_{2}\leq C(\|u_{n}\|_{H^{2}}+1).$

Now we will show the energy minimization scheme has a limit, as the time step
$k\rightarrow 0$, that solves the parabolic variational inequality (2.3). We
study the energy minimization scheme as in previous sections and by
Proposition 2.9, it suffices to show the following lemma.

######  Lemma 2.12 (Rothe’s method).

Recall that $k$ is the small time step and suppose for each $j=1,\cdots M$,
where $M=T/k$, $u_{j}$ is the unique minimizer of $E_{j}(u)$ in
$H^{1}_{+}(\Omega)$. Here $E_{j}(u)$ are defined similarly as in Theorem 2.6:

$\displaystyle E_{j}(u)\coloneqq\int_{\Omega}\frac{1}{2}|\nabla
u|^{2}+\frac{1}{2k}(u-u_{j-1})^{2}+u.$

Then $u=\lim_{k\rightarrow 0}u_{M}(x,t)$ exists and solves the parabolic
variational inequality (2.3). Here $u_{M}$ is the associated linear
interpolation defined by

$u_{M}(x,t)\coloneqq(1-\theta)\cdot u_{j}(x)+\theta\cdot u_{j+1}(x)\ ,\
\mbox{for}\ t=(j+\theta)k,\ \theta\in[0,1).$

###### Proof of Lemma 2.12.

Our proof follows [15, 23]. First, note that $u_{j}$ is the unique minimizer
of $E_{j}$ and as discussed earlier in Proposition 2.9, it satisfies the
elliptic variational inequality (2.5):

$\int_{\Omega}\nabla u_{j}\cdot\nabla(v-u_{j})+\frac{u_{j}}{k}(v-u_{j})\
dx\geq\int_{\Omega}\left(\frac{u_{j-1}}{k}-1\right)\cdot(v-u_{j})\mbox{for all
}v\in H^{1}_{+}(\Omega).$ (2.8)

Taking $v=u_{j-1}$, one can derive

$\langle\nabla u_{j},\nabla(u_{j-1}-u_{j})\rangle+\frac{1}{k}\langle
u_{j},u_{j-1}-u_{j}\rangle\geq\frac{1}{k}\langle
u_{j-1},u_{j-1}-u_{j}\rangle-\langle 1,u_{j-1}-u_{j}\rangle.$

Similarly we take $v=u_{j}$ in the $j-1$-th inequality:

$\langle\nabla u_{j-1},\nabla(u_{j}-u_{j-1})\rangle+\frac{1}{k}\langle
u_{j-1},u_{j}-u_{j-1}\rangle\geq\frac{1}{k}\langle
u_{j-2},u_{j}-u_{j-1}\rangle-\langle 1,u_{j}-u_{j-1}\rangle.$

Adding the two inequalities above it follows that

$\frac{1}{k}\|u_{j}-u_{j-1}\|_{2}^{2}+\|\nabla(u_{j}-u_{j-1})\|_{2}^{2}\leq\frac{1}{k}\langle
u_{j}-u_{j-1},u_{j-1}-u_{j-2}\rangle.$

Note that when $j=1$, we choose $v=u_{0}$ and hence

$\frac{1}{k}\|u_{1}-u_{0}\|_{2}^{2}+\|\nabla(u_{1}-u_{0})\|_{2}^{2}\leq\left|\langle\nabla
u_{0},\nabla(u_{1}-u_{0})\rangle\right|+\left|\langle
1,u_{0}-u_{1}\rangle\right|\leq\left(\|\Delta
u_{0}\|_{2}+1\right)\cdot\|u_{1}-u_{0}\|_{2}.$

Therefore we obtain that

$\|\frac{u_{j}-u_{j-1}}{k}\|_{2}\leq C$

for any $j=1,\cdots,M$ and a positive absolute constant $C$. Note that
$\dot{u_{M}}(t)=\frac{u_{j+1}-u_{j}}{k}$, therefore by Arzelà-Ascoli Theorem,
$u_{M}(t)$ converges to some function $u$ in $C([0,T],L^{2}(\Omega))$. Then we
can define

$\widetilde{u_{M}}(t)=u_{j}\ ,\ \mbox{for}\ t\in[jk,(j+1)k),$

similar to Lemma 3.6 and Remark 3.2, $\widetilde{u_{M}}$ converges to the same
$u$. Indeed, $\nabla u_{M}$ converges to $\nabla u$ weakly in
$L^{2}((0,T),L^{2}(\Omega))$. As a result we can rewrite (2.8) as follows: for
any $v\in H^{1}_{+}(\Omega)$, we have

$\langle\dot{u_{M}}(t),v(t)-\widetilde{u_{M}}(t)\rangle+\langle\nabla\widetilde{u_{M}},\nabla(v-\widetilde{u_{M}})\rangle\geq-\langle
1,v-\widetilde{u_{M}}\rangle,$ (2.9)

for a.e. $t\in(0,T)$. It then implies that for arbitrary $\tau_{1}<\tau_{2}$
in $[0,T]$,

$\int_{\tau_{1}}^{\tau_{2}}\langle\dot{u_{M}}(t),v(t)-\widetilde{u_{M}}(t)\rangle+\langle\nabla\widetilde{u_{M}},\nabla(v-\widetilde{u_{M}})\rangle\
dt\geq-\int_{\tau_{1}}^{\tau_{2}}\langle 1,v-\widetilde{u_{M}}\rangle\ dt,$
(2.10)

letting $k\to 0$, we derive the desired result

$\int_{\tau_{1}}^{\tau_{2}}\langle\dot{u}(t),v(t)-u(t)\rangle+\langle\nabla
u,\nabla(v-u)\rangle\ dt\geq-\int_{\tau_{1}}^{\tau_{2}}\langle 1,v-u\rangle\
dt,$ (2.11)

for almost every $\tau_{1}<\tau_{2}$ in $[0,T]$.

∎

####  2.2.3 A regularized formulation

We introduce a formulation using a regularization method with parameter
$\epsilon$ proposed first in [2]. Here, we will see the convergence in
regularized solutions $u_{\epsilon}(x,t)$ as $\epsilon\to 0$ to the other OD
formulations. Its analysis is simplified since the approximating problems
avoid handling the free interfaces directly. We include this approach for
completeness. While there is theoretical insight to be gained from this
formulation, it is unattractive for numerical approximation for application
purposes as free interface locations are not easily identified from
$\epsilon>0$ results.

$\partial_{t}u_{\epsilon}=\Delta u_{\epsilon}-f_{\epsilon}(u_{\epsilon}),$
(2.12)

where

$f_{\epsilon}(u_{\epsilon})=\left\\{\begin{aligned}
1&&u_{\epsilon}>\epsilon\\\
\frac{u_{\epsilon}}{\epsilon}&&u_{\epsilon}\leq\epsilon,\end{aligned}\right.$
(2.13)

with same initial condition $u_{0}(x)$. Note that $f_{\epsilon}(x)$ is a
Lipschitz function and as a result $u_{\epsilon}$ exists as a smooth solution
for each $\epsilon>0$ with $u_{\epsilon}(x,t)>0$ for all $x\in\Omega$ and
$t>0$.

We consider $u_{\epsilon_{1}}$ and $u_{\epsilon_{2}}$ with
$\epsilon_{1}<\epsilon_{2}$. Denote their difference by
$w=u_{\epsilon_{1}}-u_{\epsilon_{2}}$, then

$\displaystyle\partial_{t}w-\Delta
w=-f_{\epsilon_{1}}(u_{\epsilon_{1}})+f_{\epsilon_{2}}(u_{\epsilon_{2}}).$

Note that,

$-f_{\epsilon_{1}}(u_{\epsilon_{1}})+f_{\epsilon_{2}}(u_{\epsilon_{2}})=\left\\{\begin{aligned}
0\ ,\ &&\mbox{if}\ u_{\epsilon_{1}}>\epsilon_{1}\ ,\
u_{\epsilon_{2}}>\epsilon_{2}\\\ -\frac{u_{\epsilon_{1}}}{\epsilon_{1}}+1\ ,\
&&\mbox{if}\ u_{\epsilon_{1}}\leq\epsilon_{1}\ ,\
u_{\epsilon_{2}}>\epsilon_{2}\\\ -1+\frac{u_{\epsilon_{2}}}{\epsilon_{2}}\ ,\
&&\mbox{if}\ u_{\epsilon_{1}}>\epsilon_{1}\ ,\
u_{\epsilon_{2}}\leq\epsilon_{2}\\\
-\frac{u_{\epsilon_{1}}}{\epsilon_{1}}+\frac{u_{\epsilon_{2}}}{\epsilon_{2}}\
,\ &&\mbox{if}\ u_{\epsilon_{1}}\leq\epsilon_{1}\ ,\
u_{\epsilon_{2}}\leq\epsilon_{2}\end{aligned}\right..$ (2.14)

We observe that

$\displaystyle-1+\frac{u_{\epsilon_{2}}}{\epsilon_{2}}\leq 0\ ,\ $
$\displaystyle\mbox{if}\ u_{\epsilon_{1}}>\epsilon_{1}\ ,\
u_{\epsilon_{2}}\leq\epsilon_{2}$
$\displaystyle-\frac{u_{\epsilon_{1}}}{\epsilon_{1}}+\frac{u_{\epsilon_{2}}}{\epsilon_{2}}=-\frac{w}{\epsilon_{1}}+u_{\epsilon_{2}}\cdot\left(\frac{1}{\epsilon_{2}}-\frac{1}{\epsilon_{1}}\right)\
,\ $ $\displaystyle\mbox{if}\ u_{\epsilon_{1}}\leq\epsilon_{1}\ ,\
u_{\epsilon_{2}}\leq\epsilon_{2}$ $\displaystyle
u_{\epsilon_{1}}<u_{\epsilon_{2}}\ ,\ $ $\displaystyle\mbox{if}\
u_{\epsilon_{1}}\leq\epsilon_{1}\ ,\ u_{\epsilon_{2}}>\epsilon_{2};$

so if we assume the maximal value of $w$ is achieved at $(x_{0},t_{0})$ with
$x_{0}\in\Omega$ and $t_{0}>0$, then $w(x_{0},t_{0})>0$,
$\partial_{t}w(x_{0},t_{0})=0$ and $\Delta w(x_{0},t_{0})\leq 0$. Then the
standard maximum principle gives a partial result of the following statement:

######  Theorem 2.13.

Suppose a sequence of classical functions $\\{u_{\epsilon}\\}$ solve (2.12),
then $u_{\epsilon}$ is monotonically decreasing as $\epsilon$ decreases to 0.
Moreover, the limiting function

$\lim_{\epsilon\to 0}u_{\epsilon}=u$

holds pointwisely. This limiting function $u$ solves the variational
inequality (2.3).

The details of the proof can be found in Appendix B. An alternate convergence
statement and proof is given in Appendix C.

###  2.3 Conjecture on the general dynamics in 1D

We make the following plausible conjecture for the general dynamics (including
topological changes) of the Cauchy problem in 1D with initial conditions
$u_{0}(x)\in H^{1}_{+}(\Omega)$ with compact support. Here, we consider the
problem for all space rather than half space with a no flux condition at
$x=0$.

######  Conjecture 2.14.

Assume $u_{0}$ has a finite ${\cal S}(0)$ where ${\cal S}(t)$ counts the
number of free boundary points:

${\cal S}(t)=\left\\{x:u(x,t)=0\mbox{\ and $u(y,t)>0$ for some $y$ in every
neighbourhood of $x$}\right\\}.$

Then

(i)

${\cal S}(t)$ is finite for every $t>0$.

(ii)

There exists a finite increasing sequence of times $t_{j}$, $j=0,\ldots,M$
with $t_{0}=0$ and $\operatorname{card}{\cal S}(t):=n_{j}$ constant on every
interval $(t_{j},t_{j+1})$ and $u\equiv 0$ for $t\geq t_{M}$.

(iii)

${\cal S}(t)=\\{s_{1}(t),s_{2}(t),\ldots s_{n_{j}}(t)\\}$ for $s_{l}(t)$
smooth on $(t_{j},t_{j+1})$.

(iii)

$u(x,t)$ is $C^{1}$ for $t>0$ and $C^{\infty}$ except at free boundary points.

Recent related results have been shown for the Stefan problem [9]. Similar
analysis of the OD problem is complicated by the reaction term that allows the
formation of new zones of constraint ($u\equiv 0$).

##  3 Numerical Approximation

We consider two numerical methods. The first, suitable for 1D dynamics without
topological change, is based on the mapped domain formulation described in
Section 2.1.2. The second, suitable for dynamics in higher dimensions
including topological change, is based on our new gradient formulation
described in Section 2.2.2. We prove convergence of this scheme.

###  3.1 Mapped domain ($Y$ formulation) method

We consider the discretization of the mapped domain formulation (2.2) in space
using cell centred finite differences. We first discretize in space, leaving
time continuous (known as a Method of Lines – MoL – discretization) with
approximations $u^{j}(t)\approx u((j-1/2)h,t)$, $j=1\ldots N$ where $h$ is the
uniform grid spacing with $N$ subintervals of $y\in[0,1]$. The interface
location $s(t)$ is approximated by $S(t)$.

Figure 1:  Cell centered finite difference spatial approximation of the 1D
mapped domain formulation

Boundary conditions are implemented using ghost points [29] $u^{0}(t)\approx
u(-h/2,t)$ and $u^{N+1}(t)\approx u(1+h/2,t)$ depicted in Figure 1. Boundary
conditions at $y=1$ are implemented using second order averages and
differences:

$\displaystyle(u^{N+1}+u^{N})/2$ $\displaystyle=$ $\displaystyle 0$ (3.1)
$\displaystyle(u^{N+1}-u^{N})/h$ $\displaystyle=0$ (3.2)

which implies that $u^{N}=u^{N+1}=0$. The no-flux boundary condition at $y=0$
is approximated similarly. The MoL discretization for the interior equations
is

$D_{2}u^{j}+S\dot{S}yD_{1}u^{j}-S^{2}\dot{u}^{j}-S^{2}=0.$ (3.3)

where $D_{2}$ and $D_{1}$ are the standard centered second order finite
difference operators. The system (3.1,3.2,3.3) is a Differential Algebraic
Equation (DAE) [1] and has index one. For computational results, we use
Implicit (Backward) Euler time stepping with Newton iterations for the
resulting nonlinear system at each time step. In a computational study, we
observe errors of size $O(h^{2})+O(k)$ where $k$ is the time step, as expected
for a second order spatial and first order temporal discretization.

######  Remark 3.1.

The convergence of the method has not been proved. The missing direct analysis
discussed in Remark 2.2 could give insight.

####  3.1.1 Computational results

Examples of the dynamics computed with the DAE formulation in the mapped
region are shown in Figure 2. The left figure shows the solution with initial
conditions $u_{0}(x)=(1-x)^{2}/2$ for $x\in[0,1]$ considered often in the
literature. It is the steady state of the problem forced with flux condition
$u_{x}=-1$ at $x=0$ [6]. In this solution, $s(t)$ moves monotonically to the
left. The solution in this formulation ends when $s(T)=0$ ($u\equiv 0$). A
specialized method in this general framework was developed in [20] to
accurately compute both the solution and the end time $T$ of the dynamics. Our
mapped formulation breaks down as $t\rightarrow T$. The right computation of
Figure 2 with initial conditions $u_{0}=-x^{4}+3x^{3}-5x^{2}/2+1/2$ for
$x\in[0,1]$ shows that $s(t)$ does not have to be monotone decreasing. Here,
$s(t)$ initially moves to the right driven by diffusion and then to the left
as $u$ values decrease due to the consumption term.

Figure 2:  1D solutions of the OD problem without topological change. Left:
initial conditions $u_{0}(x)=(1-x)^{2}/2$ (standard problem in the
literature), $s(t)$ decreases monotonically. Right: initial conditions
$u_{0}(x)=-x^{4}+3x^{3}-5x^{2}/2+1/2$, $s(t)$ initially moves to the right
driven by diffusion and then to the left as $u$ values decrease due to the
consumption term.

###  3.2 Gradient flow method

In this section, we continue the discretization of the gradient flow
formulation from Section 2.2.2 and discretize in space with
$u^{i,j}_{n}\approx u(ih,jh,nk)$. We consider the discretization in two
spatial dimension for ease of presentation but the argument extends to other
dimensions. The energy minimization problem (2.5) is approximated by the
discrete minimization of

$E_{n+1}^{N}=\frac{h^{2}}{2}\sum_{i,j=0}^{N-1}(\left(\frac{u^{i+1,j}-u^{i,j}}{h}\right)^{2}+\left(\frac{u^{i,j+1}-u^{i,j}}{h}\right)^{2})+h^{2}\sum_{i,j=0}^{N}\left(\frac{1}{2k}\left(u^{i,j}-u_{n}^{i,j}\right)^{2}+u^{i,j}\right),$
(3.4)

where $N$ the number of grid points with $N=1/h$, assuming without loss of
generality that all positive values of $u$ are captured in $(0,1)\times(0,1)$.
We solve this minimization problem subject to all non-negative discrete data

$\vec{u}_{n+1}^{N}\coloneqq(u^{1,1}_{n+1},u^{2,1}_{n+1},\cdots,u^{N,1}_{n+1},u^{1,2}_{n+1},\cdots,u^{N,N}_{n+1}).$

This is a convex, quadratic minimization problem with linear, inequality
constraints and so has a unique global minimum. We show below that the
solution to the discrete optimization problem converges to the OD solutions as
$h,k\rightarrow 0$. In Section 3.2.1 we discuss the technique we use to solve
the discrete optimization problem. Denoting $M=T/k$, we use
$\\{\vec{u}_{n}\\}_{n=1}$ to define an approximate solution:

$u_{N,M}(x,t)\coloneqq\left\\{\begin{aligned} &&(1-t)\cdot u_{0}(x)+t\cdot
u_{1}^{N}(x)&&\mbox{for}\ t\in[0,k)&&\\\ &&\cdots&&\\\ &&(1-t)\cdot
u_{m}^{N}(x)+t\cdot u_{m+1}^{N}(x)&&\mbox{for}\ t\in[mk,(m+1)k)&&\\\
&&\cdots&&\\\ &&(1-t)\cdot u_{M-1}^{N}(x)+t\cdot u_{M}^{N}(x)&&\mbox{for}\
t\in[(M-1)k,T],&&\end{aligned}\right.$ (3.5)

where $u_{0}(x)$ is the initial condition and for $1\leq m\leq M$ and
$u_{m}^{N}$ is the linear (bilinear) approximation function. In particular
when $x=(x_{1},x_{2})\in[ih,(i+1)h)\times[jh,(j+1)h)$ we define $u_{m}^{N}$ as
follows:

$\displaystyle u_{m}^{N}(x_{1},x_{2})=$ $\displaystyle
h^{-2}\Big{(}u_{m}^{i,j}((i+1)h-x_{1})((j+1)h-x_{2})+u_{m}^{i,j+1}((i+1)h-x_{1})(x_{2}-jh)$
(3.6)
$\displaystyle+u_{m}^{i+1,j}(x_{1}-ih)((j+1)h-x_{2})+u_{m}^{i+1,j+1}(x_{1}-ih)(x_{2}-jh)\Big{)}.$

Therefore we have when $(x_{1},x_{2})\in(ih,(i+1)h)\times(jh,(j+1)h)$

$\displaystyle\partial_{1}u_{m}^{N}(x_{1},x_{2})=h^{-2}\left((u^{i+1,j}_{m}-u^{i,j}_{m})((j+1)h-x_{2})+(u^{i+1,j+1}_{m}-u^{i,j+1}_{m})(x_{2}-jh)\right),$
$\displaystyle\partial_{2}u_{m}^{N}(x_{1},x_{2})=h^{-2}\left((u^{i,j+1}_{m}-u^{i,j}_{m})((i+1)h-x_{1})+(u^{i+1,j+1}_{m}-u^{i+1,j}_{m})(x_{1}-ih)\right).$

Therefore, by direct computation we obtain that

$\displaystyle\int_{jh}^{(j+1)h}\int_{ih}^{(i+1)h}(\partial_{1}u_{m}^{N})^{2}\
dx_{1}dx_{2}=$
$\displaystyle\frac{1}{3}(u_{m}^{i+1,j}-u^{i,j}_{m})^{2}+\frac{1}{3}(u_{m}^{i+1,j+1}-u_{m}^{i,j+1})^{2}$
(3.7)
$\displaystyle+\frac{1}{3}(u_{m}^{i+1,j+1}-u_{m}^{i,j+1})(u_{m}^{i+1,j}-u^{i,j}_{m}).$

Moreover the pointwise limit

$u_{m}=\lim_{N\to\infty}u_{m}^{N}$ (3.8)

exists. For convenience we also define

$\displaystyle\widetilde{u_{m}^{N}}(x)\coloneqq\sum_{i,j=1}^{N}u^{i,j}_{m}\cdot\chi_{((i-1)h,ih)}(x_{1})\chi_{((j-1)h,jh)}(x_{2}),$
(3.9)

where $\chi_{I}$ is the characteristic function on the interval I. We also
define

$u_{M}(x,t)\coloneqq(1-t)\cdot u_{m}(x)+t\cdot u_{m+1}(x)\ ,\ \mbox{for}\
t\in[mk,(m+1)k).$ (3.10)

######  Remark 3.2.

We observe that both approximations $u_{m}^{N}(x)$ and
$\widetilde{u_{m}^{N}}(x)$ will converge to the same limit in $L^{2}(\Omega)$
as $N\to\infty$ and similarly for $u_{N,M}(x,t)$ and
$\widetilde{u_{N,M}}(x,t)$ as $N\to\infty$. Arzelà-Ascoli theorem and the
finite energy assumption then imply the uniform convergence as in Lemma 2.12.

######  Theorem 3.3.

Suppose $\vec{u}_{n+1}^{N}$ solves the discrete minimization problem (3.4),
then

$u(x,t)=\lim_{M\to\infty}\lim_{N\to\infty}u_{N,M}$

with $h=1/N$ and $k=T/M$ exists in $\mathcal{J}$ and $u$ is the solution to
the variational inequality (2.3) that is

$\displaystyle\int_{0}^{t}\int_{\Omega}u_{t}\cdot(v-u)+\int_{0}^{t}\int_{\Omega}\nabla
u\cdot\nabla(v-u)\geq\int_{0}^{t}\int_{\Omega}u-v;\ \mbox{for all}\
v\in\mathcal{J},\ \mbox{a.e.}\ t\in(0,T).$

The proof relies on two lemmas. To start with, we give definitions of gamma
convergence of energy functionals shown in Lemma 3.6 as given in [7]:

######  Definition 3.4 (Gamma convergence).

We say that the sequence of functionals
$\\{\mathcal{E}_{l}\\}:X\to\mathbb{R}\cup\\{-\infty,+\infty\\}$ where $X$ is a
metric space, $\Gamma$-converges to $\mathcal{E}$ if the following conditions
are satisfied:

i

whenever $x_{l}\to x$, $\mathcal{E}(x)\leq\liminf_{l}\mathcal{E}_{l}(x_{l})$;

ii

for any $x\in X$, there exists $x_{l}\to x$ in $X$ such that
$\limsup_{l}\mathcal{E}_{l}(x_{l})\leq\mathcal{E}(x)$.

The following is a relevant property of $\Gamma$-convergence:

######  Proposition 3.5.

Given a metric space $X$ and suppose a sequence of functionals defined
$\mathcal{E}_{l}$ defined in $X$ $\Gamma$-converges to $\mathcal{E}$. Assume
that for each $l$, $x_{l}$ is a minimizer of $\mathcal{E}_{l}$, and if
$\overline{x}$ is a cluster point of $\\{x_{l}\\}$, then $\overline{x}$ is a
minimizer of $\mathcal{E}$.

We refer the proof to [7] (Corollary 7.20.). We then consider the following
energy functional:

$E_{n+1}(u_{n+1})=\int_{\Omega}\frac{1}{2}|\nabla
u_{n+1}|^{2}+\frac{1}{2k}(u_{n+1}-u_{n})^{2}+u_{n+1}$

where $u_{n+1}(x)$ is defined in (3.8).

######  Lemma 3.6 (Gamma convergence of discrete functionals).

For each $n$, $E_{n+1}^{N}$ $\Gamma$-converges to $E_{n+1}$ as $N\to\infty$ or
equivalently $h\to 0$ in $L^{2}(\Omega)$.

###### Proof of Lemma 3.6.

We follow the proof in [7].

To show (i): let $u_{n+1}^{N}\in L^{2}(\Omega)$ such that $\liminf
E_{n+1}^{N}(u_{n+1}^{N})<+\infty$ and therefore there exists a subsequence
$u_{n+1}^{N_{l}}$ such that $\lim E_{n+1}^{N_{l}}(u_{n+1}^{N_{l}})=\liminf
E_{n+1}^{N}(u_{n+1}^{N})$. For each $l$, there exists a mesh of grid points
and a vector
$\vec{u}_{n+1}^{N_{l}}$($\in\mathbb{R}^{N_{l}+2}\times\mathbb{R}^{N_{l}+2}$ in
the 2D Neumann boundary condition case) such that the corresponding
$u_{n+1}^{N_{l}}(x)$ is defined in (3.5)-(3.6). Then by the previous Remark
3.2, both $u_{n+1}^{N_{l}}$ and $\widetilde{u_{n+1}^{N_{l}}}$ converge to the
same limit $u$ in $L^{2}$. By (3.7) we also have

$\displaystyle
h^{2}\sum_{i,j=0}^{N_{l}-1}(\left(\frac{u_{n+1}^{i+1,j}-u_{n+1}^{i,j}}{h}\right)^{2}+\left(\frac{u_{n+1}^{i,j+1}-u_{n+1}^{i,j}}{h}\right)^{2})\geq\int_{\Omega}|\nabla
u_{n+1}^{N_{l}}|^{2}.$

Thus

$\int_{\Omega}|\nabla u|^{2}\leq\lim_{l}\int_{\Omega}|\nabla
u_{n+1}^{N_{l}}|^{2}\leq\liminf_{N}h^{2}\sum_{i,j=0}^{N-1}(\left(\frac{u^{i+1,j}_{n+1}-u^{i,j}_{n+1}}{h}\right)^{2}+\left(\frac{u^{i,j+1}_{n+1}-u^{i,j}_{n+1}}{h}\right)^{2}).$

On the other hand,

$h^{2}\cdot\sum_{i,j=0}^{N_{l}}\left(\frac{1}{2k}\left(u^{i,j}_{n+1}-u_{n}^{i,j}\right)^{2}+u_{n+1}^{i,j}\right)=\int_{\Omega}\frac{1}{2k}(\widetilde{u_{n+1}^{N_{l}}}-\widetilde{u_{n}^{N_{l}}})^{2}+\widetilde{u_{n+1}^{N_{l}}}.$

Applying the uniform convergence we obtain that

$\int_{\Omega}\frac{1}{2k}(u-u_{n})^{2}+u\leq\lim_{k}\int_{\Omega}\frac{1}{2k}(\widetilde{u_{n+1}^{N_{l}}}-\widetilde{u_{n}^{N_{l}}})^{2}+\widetilde{u_{n+1}^{N_{l}}}\leq\liminf_{N}h^{2}\cdot\sum_{i,j=0}^{N}\left(\frac{1}{2k}\left(u_{n+1}^{i,j}-u_{n}^{i,j}\right)^{2}+u_{n+1}^{i,j}\right).$

These two estimates lead to $E_{n+1}(u)\leq\liminf E_{n+1}^{N}(u_{n+1}^{N})$.

It remains to prove (ii): suppose $u\in L^{2}$ with $E_{n+1}(u)<+\infty$, so
$u\in H^{1}$. We then define $u^{i,j}_{n+1}\coloneqq
u(\frac{i}{N},\frac{j}{N})$ which defines the vector $\vec{u}_{n+1}^{N}$ with
the piecewise linear (bilinear) approximation $u_{n+1}^{N}(x)$ and piecewise
constant approximation $\widetilde{u_{n+1}^{N}}(x)$. By the finite energy
assumption, Arzelà-Ascoli theorem then guarantee the uniform convergence as in
Remark 3.2. It then follows that

$\limsup E_{n+1}^{N}(u^{N}_{n+1})\leq E_{n+1}(u).$

The 1D and 3D cases can be treated similarly. ∎

As a result of Lemma 3.6 and Proposition 3.5, we obtain the following
corollary immediately:

######  Corollary 3.7.

Suppose $u^{N}_{n+1}$ are minimizers of $E_{n+1}^{N}$ then $u_{n+1}^{N}$
converges to a function $u_{n+1}$ in $L^{2}(\Omega)$ up to a subsequence as
$h\to 0$ and such $u_{n+1}$ is the minimizer of $E_{n+1}$.

Now that $u_{m}\coloneqq\lim_{N}u_{m}^{N}$ is the minimizer of the continuous
functional $E_{m}$ for $m=1,\cdots,M$; it remains to show that
$u(x,t)=\lim_{M\to\infty}u_{M}(x,t)$ solves the variational inequality (2.3).
Recalling the Rothe’s Method (Lemma 2.12) and combining results of Lemma 3.6
and Lemma 2.12, we therefore complete the proof of Theorem 3.3.

####  3.2.1 Discrete Optimization Scheme

We consider the details of the discrete optimization problem (3.4) and present
the scheme in the 2D case. (Note that this scheme holds in 1D and 3D
similarly.) The corresponding Lagrangian problem is

$\displaystyle-\Delta_{h}u+\frac{u}{k}+\lambda=\frac{u_{n}}{k}-1,$
$\displaystyle\lambda_{(j_{1},j_{2})}<0,\ u^{(j_{1},j_{2})}=0$
$\displaystyle\forall(j_{1},j_{2})\in J$
$\displaystyle\lambda_{(i_{1},i_{2})}=0,\ u^{(i_{1},i_{2})}\geq 0$
$\displaystyle\forall(i_{1},i_{2})\in I,$

where $I$ and $J$ are a disjoint partition of the grid points and $\Delta_{h}$
is the finite difference Laplacian. The $u$ in the problem is the grid vector
at the next time step $u_{n+1}$. The partitions divide those points $J$ where
the values are at the constraint and those points $I$ (“$I$” for inactive
constraint) with positive solution values where the corresponding derivative
of $E^{N}$ must be zero. The method is an active set method, where the sets
$J$ and $I$ are updated iteratively at each time step. Note that
$\lambda_{(j_{1},j_{2})}<0$ for $(j_{1},j_{2})\in J$ corresponds to $\partial
E^{N}/\partial u^{(j_{1},j_{2})}>0$, a necessary and sufficient condition for
optimality (the KKT conditions [28]). There are many techniques available to
solve such quadratic optimization problems with linear inequality constraints.
We take advantage of the simple structure of the problem and the fact that
there is little change in the index sets from one time step to the next in the
following algorithm. It is an iterative algorithm with vectors $u^{(m)}$,
$\lambda^{(m)}$ at each iteration. The matrix $A=I/k-\Delta_{h}$, where $I$ is
the identity.

##### Algorithm

Step 1

Initialize $u^{(0)}\geq 0$ (component-wise),
$\lambda^{(0)}=\min\\{0,\frac{u_{n}}{k}-1-Au^{(0)}\\}$. Set $m=0$. Repeat
steps 2-5 until the convergence criteria in step 3 is reached.

Step 2

Construct the index sets

$\displaystyle J^{(m)}=\\{(j_{1},j_{2}):\lambda^{(m),(j_{1},j_{2})}<0\\},$
$\displaystyle I^{(m)}=\\{(j_{1},j_{2}):\lambda^{(m),(j_{1},j_{2})}=0\\}.$

For any $(i_{1},i_{2})\in I^{(m)}$ such that $u^{(m),(i_{1},i_{2})}<0$, move
$(i_{1},i_{2})$ to $J^{(m)}$.

Step 3

If $J^{(m)}=J^{(m-1)}$, the solution $u=u^{(m)}$. Stop.

Step 4

Solve for $u^{(m+1)}$ and $\lambda$ using

$\displaystyle Au^{(m+1)}+\lambda=\frac{u_{n}}{k}-1,$ $\displaystyle\lambda=0\
\mbox{on}\ I^{(m)},$ $\displaystyle u^{(m+1)}=0\ \mbox{on}\ J^{(m)}.$

This is equivalent to solving sequentially for $(u^{(m+1)},\lambda)$ that
satisfy

$\displaystyle A_{II}u^{(m+1)}_{I}=(\frac{u_{n}}{k}-1)_{I},$ $\displaystyle
u^{(m+1)}_{J}=0,$ $\displaystyle\lambda=\frac{u_{n}}{k}-1-Au^{(m+1)}.$

Here vector subscripts $I$ and $J$ give the sub-vectors with those components
and $A_{II}$ is the block of the matrix A corresponding to the $I$ components.

Step 5

Update $\lambda^{(m+1)}=\min\\{0,\lambda\\}$. Increment $m$.

######  Theorem 3.8.

Let $0\leq u^{(0)}\leq u$ (component-wise). The algorithm above converges in
finitely many steps.

###### Proof.

A proof is found following closely the ideas from [14] for a similar approach
to the elliptic obstacle problem. Monotone behaviour in the index sets
$I^{(m)}$ is shown and since $N$ is finite, the algorithm converges in finite
steps. Use is made of the properties that the sub-matrix $A_{II}^{-1}$ has
positive entries ($A_{II}$ is monotone) and $A_{IJ}$ has non-positive entries
(values zero or $-1/h^{2}$) for any index sets $I$ and $J$. ∎

######  Remark 3.9.

While the proof of iteration convergence above is limited to starting
conditions $0\leq u^{(0)}\leq u$, we implement the method with $u^{(0)}=u_{n}$
and starting index sets from the converged iterations at time step $n$. This
initialization falls out of the scope of the analysis but works well (no
failures, few iterations) in practice.

######  Remark 3.10.

Similar index (active set) iteration methods have been used in capturing
methods for other implicit boundary value problems. Two of these are discussed
in Section 4. A general theory for the convergence of these iteration
strategies is not known, but they can perform well in practice.

####  3.2.2 Numerical Results

We show results in 1D with topological change in Figure 3. Initial conditions
are

$u_{0}(x)=\left\\{\begin{array}[]{ll}((1/3-x).^{2}+0.05)/((5/12)^{2}+0.05)/16&\mbox{\
$x\in[0,3/4]$}\\\ (1-x)^{2}&\mbox{\ $x\in[3/4,1]$}\\\ 0&\mbox{\ $x\geq
1$}\end{array}\right.$

A 2D example is shown in Figure 4. This example has more complicated
topological changes described in the figure caption. Based on evidence from
other computations, the limiting circular shape is generic.

Figure 3:  1D solution of the OD problem with topological changes with the
gradient flow method. Initial conditions are given in Section 3.2.2.

Figure 4: A 2D solution of the OD problem with topological change computed
with the gradient flow method. Solution contours are shown with time
increasing from top left to bottom right. Overall solution levels decrease in
time with depletion. Two topological changes occur. The first is the merger of
the two disjoint sets of $u>0$, the second is the disappearance of the centre
$u=0$ set.

##  4 Other Implicit Free Boundary Value Problems

###  4.1 A biharmonic problem

The OD problem is the simplest second order implicit free boundary problem.
The simplest fourth order problem is the following biharmonic problem shown in
1D for $u(x,t)$:

$u_{t}=-u_{xxxx}-1$

with conditions $u=0$, $u_{x}=0$, and $u_{xxx}=0$ at the implicitly defined
free boundary $x=s(t)$ and $u\equiv 0$ for $x>s(t)$. This can be derived from
the scaled, linear, viscoelastic motion of a beam above a flat, rigid surface.
Note that another boundary value problem occurs if $u_{xxx}=0$ is replaced by
$u_{xx}=0$. However, the third order condition is correct for this application
[25] and also gives the gradient flow structure described below.

We consider the time discretization of this problem as in Section 2.2.2 and
see that it is a discrete $L_{2}$ gradient flow on the energy

$\mathcal{E}^{n}\coloneqq\int\frac{1}{2}|\Delta u_{n}|^{2}+u_{n}$

with $u_{n}\in H_{+}^{2}$. We form a fully discrete scheme as was done in
Section 3.2 and compute the discrete optimization at each time step using
index iterations as described in Section 3.2.1. The convergence of the method
follows the same ideas as presented for the OD problem. Some computational
results are shown in Figure 5.

Figure 5: Two computations at three times each for the biharmonic free
boundary value problem with physical boundary conditions $u(0)=1$, $u_{xx}=0$
approaching the analytic steady state solution shown in dark blue.

######  Remark 4.1.

There has been considerable mathematical interest in the elliptic obstacle
problem as discussed in the introduction. This is the steady state of the OD
problem with nonzero physical boundary conditions. The steady state of the
biharmonic problem (in higher dimensions) described in this section would also
be mathematically interesting. Its analysis would be complicated by the lack
of a maximum principle.

###  4.2 Vector problems

The free boundary in complex fluids with yield stress is of implicit type and
is well studied [11]. Numerical approaches include regularization (increased
viscosity in the unyielded region) and an Augmented Lagrangian approach to the
non-smooth optimization problem that comes from a discretization of a
variational inequality formulation. The literature on this problem is focussed
on capturing the unyielded region rather than considering the free boundary
directly.

Implicit free boundaries in porous media flow can occur when phase change is
present. Boundaries between dry and two-phase (where there is liquid and
vapour present) regions were studied in [12, 4]. The work in [4] had important
implications to simulations of water management in fuel cells. However, many
theoretical questions were left unanswered and this became the motivation of
the corresponding author to attempt the current work.

We present below a class of implicit free boundary value problems that
generalizes the OD problem. The problems are presented in 1D with a single
free boundary at $x=s(t)$ with ${\bf u}^{l}(x,t)$ having $n$ components for
$x<s(t)$ and ${\bf u}^{r}(x,t)$ having $m$ components for $x>s(t)$. Near the
interface we take

${\bf u}^{*}_{t}=D^{*}{\bf u}^{*}_{xx}+{\bf a}^{*}$

for $*\in\\{l,r\\}$, $D^{*}$ positive diagonal matrices, and ${\bf a}^{*}$
constant vectors. At the boundary, we take

$B\left[\begin{array}[]{c}{\bf u}^{l}\\\ {\bf u}^{l}_{x}\\\ {\bf u}^{r}\\\
{\bf u}^{r}_{x}\end{array}\right]={\bf 0}$

where $B$ is an $(m+n+1)\times(2m+2n)$ matrix of full rank. This class can be
reached from a wider class by taking affine combinations of solution
components and $x$, and as an approximation of some nonlinear problems. A
problem statement can be made by adding far field conditions, $n$ on the left
and $m$ on the right. With these far field conditions we label the class as
n+m implicit free boundary value problems. The OD problem is the only well
defined example of the 1+0 class. The model in [4] is of class 2+2, although
one of the components has degenerate diffusion at the free boundary.

There are several open questions related to problems of this type motivated by
the current work on the OD problem. Which lead to well defined problems? (this
could depend on the sign of entries of $\bf a$ as discussed in Remark 2.2).
Which have gradient flow or variational inequality structure? Which allow a
capturing formulation with index iteration similar to that described in
Section 3.2.1? (true of the model in [4]).

##  5 Summary

This work summarizes the ways the Oxygen Depletion problem has been considered
in the literature: with interfaces to be tracked, captured, or found as a
limit of regularized problems. We fill in a gap in the list of formulations,
showing that the OD problem can be considered as a gradient flow with
constraint. A new numerical capturing method based on the gradient flow
formulation is proposed and a convergence proof given. The equivalence of all
formulations is shown. A biharmonic implicit free boundary value problem and a
class of vector problems are introduced.

Several open problems have been presented in the work and are summarized here:

  * •

The regularity of boundary point positions in 1D (Conjecture 2.14) and higher
dimensions.

  * •

A direct analysis of the mapped domain formulation discussed in Remark 2.2
that would possibly extend to a convergence proof of its numerical
approximation (Section 3.1) and an understanding of the general class of
vector problems in Section 4.2.

  * •

A convergence proof for the index iteration (active set) approach that has
been successful in discretizations of implicit free boundary value problems,
see Remarks 3.9 and 3.10.

  * •

A study of the biharmonic obstacle problem discussed in Remark 4.1.

  * •

An understanding of the general class of vector problems introduced in Section
4.2.

We hope the reader will find some of these problems of interest.

## Acknowledgement

Cheng is partially supported by the International Doctoral Fellowship (IDF)
provided by the University of British Columbia and Shanghai ”Super Postdoc”
Incentive Plan. Fu is supported by a PhD Fellowship from the South University
of Science and Technology. Wetton is supported by an NSERC Canada research
grant. We thank Gwynn Elfring for pointing out the literature for the boundary
conditions for the viscoelastic beam problem.

##  Appendix A Proof of Proposition 2.8

With help of the minimality of $u$, we consider a competing function
$u+\varepsilon\phi$ where $\phi$ is an arbitrary smooth function that is
compactly supported inside $\Omega$. By the definition of $\tilde{E}[u]$, it
follows that

$\displaystyle\tilde{E}[u+\varepsilon\phi]\geq\tilde{E}[u],$

that is

$\displaystyle\varepsilon\int_{\Omega}\nabla
u\cdot\nabla\phi+\frac{\varepsilon^{2}}{2}\int_{\Omega}|\nabla\phi|^{2}+\frac{\varepsilon}{k}\int_{\Omega}u\phi+\frac{\varepsilon^{2}}{2k}\int_{\Omega}\phi^{2}\geq-\int_{\Omega}(1-\frac{u_{n}}{k})[(u+\varepsilon\phi)^{+}-u].$
(A.1)

Note that

$\displaystyle\int_{\Omega}(1-\frac{u_{n}}{k})[(u+\varepsilon\phi)^{+}-u]=\varepsilon\int_{\\{u+\varepsilon\phi\geq
0\\}}(1-\frac{u_{n}}{k})\phi-\int_{\\{u+\varepsilon\phi<0\\}}(1-\frac{u_{n}}{k})u,$

ignoring the $O(\varepsilon^{2})$ terms in (A.1), we have

$\displaystyle\varepsilon\int_{\Omega}\nabla
u\cdot\nabla\phi+\frac{\varepsilon}{k}\int_{\Omega}u\phi-\varepsilon\int_{\\{u+\varepsilon\phi\geq
0\\}}(1-\frac{u_{n}}{k})^{-}\phi+\int_{\\{u+\varepsilon\phi<0\\}}(1-\frac{u_{n}}{k})^{-}u$
(A.2) $\displaystyle\geq-\varepsilon\int_{\\{u+\varepsilon\phi\geq
0\\}}(1-\frac{u_{n}}{k})^{+}\phi+\int_{\\{u+\varepsilon\phi<0\\}}(1-\frac{u_{n}}{k})^{+}u.$

In fact we have

$0\leq\int_{\\{u+\varepsilon\phi<0\\}}(1-\frac{u_{n}}{k})^{\pm}u<-\varepsilon\int_{\\{u+\varepsilon\phi<0\\}}(1-\frac{u_{n}}{k})^{\pm}\phi,$

hence (A.2) turns out to be

$\displaystyle\int_{\Omega}\nabla
u\cdot\nabla\phi+\frac{1}{k}\int_{\Omega}u\phi-\int_{\\{u+\varepsilon\phi\geq
0\\}}(1-\frac{u_{n}}{k})^{-}\phi-\int_{\\{u+\varepsilon\phi<0\\}}(1-\frac{u_{n}}{k})^{-}\phi\geq-\int_{\\{u+\varepsilon\phi\geq
0\\}}(1-\frac{u_{n}}{k})^{+}\phi.$

Moreover, we also recall that $u\geq 0$, then in $L^{1}$ sense as
$\varepsilon\to 0$,

$\left\\{\begin{aligned} &\chi_{\\{u+\varepsilon\phi\geq
0\\}}\to\chi_{A_{\phi}\cup\\{u>0\\}}\\\
&\chi_{\\{u+\varepsilon\phi<0\\}}\to\chi_{\\{u=0\\}\cap\\{\phi<0\\}},\end{aligned}\right.$

where $A_{\phi}\coloneqq\\{u=0\\}\cap\\{\phi\geq 0\\}$. Clearly, $A_{\phi}$
and $\\{u>0\\}$ are disjoint. This leads to

$\displaystyle\int_{\Omega}\nabla
u\cdot\nabla\phi+\frac{1}{k}\int_{\Omega}u\phi-\int_{\Omega}\chi_{A_{\phi}\cup\\{u>0\\}}(1-\frac{u_{n}}{k})^{-}\phi-\int_{\Omega}\chi_{\\{u=0\\}\cap\\{\phi<0\\}}(1-\frac{u_{n}}{k})^{-}\phi$
$\displaystyle\geq-\int_{\Omega}\chi_{A_{\phi}\cup\\{u>0\\}}(1-\frac{u_{n}}{k})^{+}\phi,$

or equivalently,

$\displaystyle\int_{\Omega}\nabla
u\cdot\nabla\phi+\frac{1}{k}\int_{\Omega}u\phi+\int_{\Omega}\chi_{A_{\phi}\cup\\{u>0\\}}(1-\frac{u_{n}}{k})\phi-\int_{\Omega}\chi_{\\{u=0\\}\cap\\{\phi<0\\}}(1-\frac{u_{n}}{k})^{-}\phi\geq
0.$ (A.3)

Define a distribution

$T(\phi)\coloneqq\int_{\Omega}\nabla
u\cdot\nabla\phi+\frac{1}{k}\int_{\Omega}u\phi+\int_{\Omega}\chi_{\\{u>0\\}}(1-\frac{u_{n}}{k})\phi,$

then by (A.3),

$\displaystyle
T(\phi)\geq-\int_{A_{\phi}}(1-\frac{u_{n}}{k})\phi+\int_{\\{u=0\\}\cap\\{\phi<0\\}}(1-\frac{u_{n}}{k})^{-}\phi.$

Since $\phi$ is arbitrary, we may replace it with $-\phi$ and as a result,

$\left\\{\begin{aligned}
&T(\phi)\geq-\int_{A_{\phi}}(1-\frac{u_{n}}{k})\phi+\int_{\\{u=0\\}\cap\\{\phi<0\\}}(1-\frac{u_{n}}{k})^{-}\phi\\\
&T(\phi)\leq-\int_{\\{u=0\\}\cap\\{\phi\leq
0\\}}(1-\frac{u_{n}}{k})\phi+\int_{\\{u=0\\}\cap\\{\phi>0\\}}(1-\frac{u_{n}}{k})^{-}\phi.\end{aligned}\right.$
(A.4)

Therefore, $|T(\phi)|\leq C||\phi||_{\infty}$ for some positive constant $C$,
thus by a density argument we derive that $T$ is a radon measure, i.e. there
exists a density function $\rho(x)$ such that

$T(\phi)=\int_{\Omega}\rho\phi\ dx.$

However, by (A.4), we get $\rho=0$ a.e. in $\\{u>0\\}$; moreover, by
definition of $T$ we get $\rho=0$ a.e. in $\\{u=0\\}$. This shows that
$T(\phi)=0$, or

$-\Delta u+\frac{1}{k}u+\chi_{\\{u>0\\}}(1-\frac{u_{n}}{k})=0$

in the weak sense. Equivalently,

$\frac{u-u_{n}\cdot\chi_{\\{u>0\\}}}{k}=\Delta u-\chi_{\\{u>0\\}}.$

##  Appendix B Proof of Theorem 2.13

As the discussion in Section 2.2.3 above showed, $u=\lim_{\epsilon\to
0}u_{\epsilon}$ exists pointwisely by monotonicity. It remains to show $u$ is
the solution to (2.3), that is

$\displaystyle\int_{0}^{t}\int_{\Omega}\partial_{t}u\cdot(v-u)+\int_{0}^{t}\int_{\Omega}\nabla
u\cdot\nabla(v-u)\geq\int_{0}^{t}\int_{\Omega}u-v;\ \mbox{for all}\
v\in\mathcal{J},\ \mbox{a.e.}\ t\in(0,T).$

Intuitively, suppose that $f$ is a smooth approximation, then by maximum
principle $|\nabla u_{\epsilon}|\leq\sup_{\Omega}|\nabla u_{0}|$ for any
$x\in\Omega$ and $\epsilon>0$. Thus $|\nabla u|\leq\sup|\nabla u_{0}|$,
therefore by Dini’s Theorem, such convergence is uniform and as a result,
$u\in\mathcal{J}$ because $u$ also satisfies the boundary condition and
initial condition. Once we have such uniform boundedness of $\nabla
u_{\epsilon}$, $\nabla u_{\epsilon}$ converges to $\nabla u$ weakly and as a
result,

$\lim_{\epsilon}\int_{0}^{t}\int_{\Omega}\nabla
u_{\epsilon}\cdot\nabla(v-u)=\int_{0}^{t}\int_{\Omega}\nabla
u\cdot\nabla(v-u)$

and

$\lim_{\epsilon}\int_{0}^{t}\int_{\Omega}-f_{\epsilon}(u_{\epsilon})\cdot(v-u)=-\int_{0}^{t}\int_{\Omega}\chi_{\\{u>0\\}}\cdot(v-u)=-\int_{0}^{t}\int_{\Omega}v-u+\int_{0}^{t}\int_{\Omega}\chi_{\\{u=0\\}}\cdot
v.$

Indeed we have weak convergence of $\partial_{t}u_{\epsilon}$ thanks to the
equation:

$\lim_{\epsilon}\int_{0}^{t}\int_{\Omega}\partial_{t}u_{\epsilon}\cdot(v-u)=\lim_{\epsilon}\int_{0}^{t}\int_{\Omega}-f(u_{\epsilon})\cdot(v-u)-\nabla
u_{\epsilon}\cdot\nabla(v-u).$

Since $u_{\epsilon}$ converges to $u$ pointwisely and strongly in
$L^{2}((0,T);L^{2}(\Omega))$, then up to a subsequence

$\lim_{\epsilon}\int_{0}^{t}\int_{\Omega}\partial_{t}u_{\epsilon}\cdot(v-u)=\int_{0}^{t}\int_{\Omega}\partial_{t}u\cdot(v-u).$

Note that $v\geq 0$,

$-\int_{0}^{t}\int_{\Omega}\chi_{\\{u>0\\}}\cdot(v-u)\geq-\int_{0}^{t}\int_{\Omega}v-u.$

therefore

$\displaystyle\int_{0}^{t}\int_{\Omega}\partial_{t}u\cdot(v-u)+\int_{0}^{t}\int_{\Omega}\nabla
u\cdot\nabla(v-u)\geq\int_{0}^{t}\int_{\Omega}u-v;\ \mbox{for all}\
v\in\mathcal{J},\ \mbox{a.e.}\ t\in(0,T).$

Indeed, we only require the $H^{1}$ uniform boundedness of $u_{\epsilon}$. To
see this without using smooth $f(u_{\epsilon})$ we write down $u_{\epsilon}$
in the mild form:

$u_{\epsilon}(t)=e^{t\Delta}u_{0}+\int_{0}^{t}e^{(t-s)\Delta}(f_{\epsilon}(u_{\epsilon}))\
ds\ ,$

where $e^{t\Delta}$ represents convolution with heat kernel. As a result, for
any first order differential operator $D$ we have

$\displaystyle
Du_{\epsilon}=De^{t\Delta}u_{0}+\int_{0}^{t}De^{(t-s)\Delta}(f_{\epsilon}(u_{\epsilon}))\
ds$

and hence

$\displaystyle\|Du_{\epsilon}\|_{2}\leq\|De^{t\Delta}u_{0}\|_{2}+\int_{0}^{t}\|De^{(t-s)\Delta}f(u_{\epsilon})\|_{2}\
ds.$

Note that $e^{t\Delta}u_{0}$ solves the standard heat equation with initial
data $u_{0}$, we have

$\|De^{t\Delta}u_{0}\|_{2}=\|e^{t\Delta}Du_{0}\|_{2}\lesssim\|Du_{0}\|_{2}\lesssim
1,$

for any $t\in(0,T)$. On the other hand,

$\displaystyle\|De^{(t-s)\Delta}f(u_{\epsilon})\|_{2}\lesssim\|De^{(t-s)\Delta}f(u_{\epsilon})\|_{\infty}=|K*f(u_{\epsilon})|\
,$

where $K$ is the kernel corresponding to $De^{(t-s)\Delta}$. Since $|f|\leq
1$,

$\displaystyle|K*f(u_{\epsilon})|$
$\displaystyle\leq\|K\|_{2}\cdot\|f(u_{\epsilon})\|_{2}$
$\displaystyle\lesssim\|K\|_{2}.$

We see that from the Fourier side

$\displaystyle\|K\|_{2}^{2}$
$\displaystyle\lesssim\sum_{k\in\mathbb{Z}^{d}}|k|^{2}e^{-2(t-s)|k|^{2}}$
$\displaystyle=\sum_{|k|\geq 1}|k|^{2}e^{-2(t-s)|k|^{2}}$
$\displaystyle\lesssim\int_{1}^{\infty}e^{-2(t-s)r^{2}}r^{1+d}\ dr\ .$

For the 1D case, first we observe that

$\displaystyle\int_{1}^{\infty}e^{-2(t-s)r^{2}}r^{2}\ dr=$
$\displaystyle\frac{\sqrt{2\pi}[1-\operatorname{erf}(\sqrt{2(t-s)})]+4\sqrt{t-s}e^{-2(t-s)}}{16(t-s)^{3/2}}$
$\displaystyle\lesssim$
$\displaystyle\frac{1-\operatorname{erf}(\sqrt{2(t-s)})}{(t-s)^{3/2}}+\frac{e^{-2(t-s)}}{t-s},$

where $\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}\ dt$,
the Gauss error function. Therefore,

$\|De^{\gamma\Delta}f(u_{\epsilon})\|_{2}\lesssim\frac{\left(1-\operatorname{erf}(\sqrt{2(t-s)})\right)^{1/2}}{(t-s)^{3/4}}+\frac{e^{-(t-s)}}{(t-s)^{1/2}}.$

Now we would assume $t\geq 1$, as the other case $t<1$ is easier. Let
$\gamma=t-s$, we split the following integral into 2 parts:

$\displaystyle\int_{0}^{t}\|De^{\gamma\Delta}f(u_{\epsilon})\|_{2}\
d\gamma=\int_{0}^{1}\|De^{\gamma\Delta}f(u_{\epsilon})\|_{2}\
d\gamma+\int_{1}^{t}\|De^{\gamma\Delta}f(u_{\epsilon})\|_{2}\ d\gamma.$

(i) $\gamma>1$:

Then we have

$\frac{\left(1-\operatorname{erf}(\sqrt{2\gamma})\right)^{1/2}}{\gamma^{3/4}}\lesssim\frac{e^{-\gamma}}{\gamma^{5/4}},$

thus

$\displaystyle\int_{1}^{t}\|De^{\gamma\Delta}f(u_{\epsilon})\|_{2}\ d\gamma$
$\displaystyle\lesssim\int_{1}^{t}\frac{e^{-\gamma}}{\gamma^{3/4}}+\frac{e^{-\gamma}}{\gamma^{1/2}}\
d\gamma$ $\displaystyle\lesssim\int_{1}^{t}\frac{e^{-\gamma}}{\gamma^{1/2}}\
d\gamma$
$\displaystyle\lesssim\int_{1}^{\infty}\frac{e^{-\gamma}}{\gamma^{1/2}}\
d\gamma$ $\displaystyle\lesssim 1\ .$

(ii) $\gamma\leq 1$:

We use another estimate for $\|K*f(u_{\epsilon})\|_{2}$. We compute from the
Fourier side:

$\displaystyle\|K*f(u_{\epsilon})\|_{2}^{2}$ $\displaystyle=\sum_{|k|\geq
1}|k|^{2}e^{-2\gamma|k|^{2}}|\widehat{f(u_{\epsilon})}(k)|^{2}$
$\displaystyle\leq\max_{|k|\geq
1}\left\\{|k|^{2}e^{-2\gamma|k|^{2}}\right\\}\cdot\sum_{|k|\geq
1}|\widehat{f(u_{\epsilon})}(k)|^{2}$ $\displaystyle\lesssim\max_{|k|\geq
1}\left\\{|k|^{2}e^{-2\gamma|k|^{2}}\right\\}\cdot\|f(u_{\epsilon})\|_{2}^{2}$
$\displaystyle\lesssim\max_{|k|\geq
1}\left\\{|k|^{2}e^{-2\gamma|k|^{2}}\right\\}.$

Define $g(x)=x^{2}e^{-2\gamma x^{2}}$, where $x\geq 0$. Then,

$g^{\prime}(x)=xe^{-2\gamma x^{2}}\left(1-2\gamma x^{2}\right)\ ,$

this shows the maximum achieves at $x=\frac{1}{\sqrt{2\gamma}}$ and hence

$g(x)\leq g(\frac{1}{\sqrt{2\gamma}})\leq\frac{1}{\gamma}\,$

thus

$\|De^{\gamma\Delta}f(u_{\epsilon})\|_{2}\lesssim\frac{1}{\sqrt{\gamma}}\ ,$

As a result,

$\displaystyle\int_{0}^{1}\|De^{\gamma\Delta}f(u_{\epsilon})\|_{2}\
d\gamma\lesssim\int_{0}^{1}\frac{1}{\sqrt{\gamma}}\
d\gamma\cdot\|f(u_{\epsilon})\|_{2}\lesssim 1\ .$

Similar arguments can be applied to the 2D and 3D cases. In what follows,

$\|Du_{\epsilon}\|_{2}\lesssim 1,$

for any $t\in(0,T)$ and the bound is independent of $\epsilon$.

##  Appendix C Another Proof of the Regularization Result

We recall the variational inequality setting (2.3), that is to solve $u\in
H^{1}_{+}(\Omega)$

$\int_{0}^{t}\langle\partial_{t}u-\Delta u+1,v-u\rangle\geq 0,\ \forall\
v\in\mathcal{J}.$

As in [13], it then has an equivalent formulation, that is to solve $u(t)$ and
$\lambda^{*}(t)$:

$\left\\{\begin{aligned} &\partial_{t}u-\Delta u+1=-\lambda^{*}(t)\geq 0\\\
&u\geq 0,\ \langle u(t),\lambda^{*}(t)\rangle=0,\ \forall\
t>0.\end{aligned}\right.$ (C.1)

To approach this, we introduce a regularized approximation family: we aim to
find $u_{c}$ for any $c>0$ such that the following holds weakly:

$\displaystyle\partial_{t}u_{c}-\Delta
u^{c}+1+\min\left(0,-1+cu^{c}\right)=0.$

By defining $\lambda^{c}=\min\left(0,-1+cu^{c}\right)$, we can rewrite the
above scheme as

$\displaystyle\partial_{t}u^{c}-\Delta u^{c}+1+\lambda^{c}=0.$

It is typical to write the regularization term in this way in some literature,
but the approach is the same as the regularization in Section 2.2.3 with
$c=1/\epsilon$. We then discretize it in time: for any $\phi\in H^{1}$, the
following holds

$\left\langle\frac{u_{n+1}^{c}-u_{n}^{c}}{k},\phi\right\rangle+\left\langle\nabla
u_{n+1}^{c},\nabla\phi\right\rangle+\langle
1,\phi\rangle+\left\langle\min(0,-1+cu_{n+1}^{c}),\phi\right\rangle=0,$ (C.2)

where $u^{c}_{0}$ is chosen to be $u_{0}$. We write $u_{n}$ instead of
$u_{n}^{c}$ for simplicity. Note that the operator
$A(u)\coloneqq\frac{u}{k}-\Delta u+\min(0,-1+cu)$ is coercive and monotone. As
a result, there exists a unique solution $u_{n+1}\in H^{1}$ for sufficiently
small $k>0$ independent of $c>0$. To show $u_{n+1}\in H^{1}_{+}(\Omega)$, we
prove by induction. Assuming $u_{n}\in H^{1}_{+}(\Omega)$, we test the (C.2)
with $(u_{n+1})^{-}$. Therefore we derive that

$\displaystyle\frac{1}{k}\langle
u_{n+1},(u_{n+1})^{-}\rangle+\left\langle\nabla
u_{n+1},\nabla(u_{n+1})^{-}\right\rangle+\langle
1,(u_{n+1})^{-}\rangle+\left\langle\min(0,-1+cu_{n+1}),(u_{n+1})^{-}\right\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{k}\langle
u_{n},(u_{n+1})^{-}\rangle\leq 0.$

We observe that $\left\langle\nabla
u_{n+1},\nabla(u_{n+1})^{-}\right\rangle=\left\langle\nabla(u_{n+1})^{-},\nabla(u_{n+1})^{-}\right\rangle\geq
0$. Moreover, $\langle
1,(u_{n+1})^{-}\rangle+\langle\min(0,-1+cu_{n+1}),(u_{n+1})^{-}\rangle=c\langle(u_{n+1})^{-},(u_{n+1})^{-}\rangle\geq
0$. We thus obtain that $\langle u_{n+1},(u_{n+1})^{-}\rangle\leq 0$ and hence
$u_{n+1}\in H^{1}_{+}(\Omega)$. We then define

$u^{c}_{M}(x,t)=u_{n}+\frac{t-nk}{k}(u_{n+1}-u_{n}),\ \mbox{for}\
t\in[nk,(n+1)k),$

where $M=T/k$. By the same argument in Lemma 2.12, we have $u^{c}_{M}$
converges to function $u^{c}$ in $L^{2}(0,T;H^{1})$ as $M\to\infty$ up to a
subsequence. In fact, it is easy to see that $u_{c}$ is the solution to (C.2).
On the other hand, we show that $u^{c}$ converges to $u^{*}$ as $c\to\infty$.

######  Theorem C.1 (Monotonicity).

Let $u_{n+1}^{c}$ and $u^{c}$ be defined as above. If $0<c\leq b$, then
$u_{n+1}^{c}\geq u_{n+1}^{b}$ for all $n=0,1,2,\cdots$. Therefore
$u^{c}(t)\geq u^{b}(t)$ as a direct application.

###### Proof.

The proof is given by induction. Suppose $u_{n}^{c}\geq u_{n}^{b}$ and for
each $n$ define $\lambda_{n}^{c}$ by

$\lambda_{n+1}^{c}=\min(0,-1+cu_{n+1}^{c}).$

Then the proof is similar to the one showing $u_{n+1}^{c}\geq 0$, we have that

$\displaystyle\frac{1}{k}\langle
u_{n+1}^{c}-u_{n+1}^{b},(u_{n+1}^{c}-u_{n+1}^{b})^{-}\rangle+\left\langle\nabla(u_{n+1}^{c}-u_{n+1}^{b}),\nabla(u_{n+1}^{c}-u_{n+1}^{b})^{-}\right\rangle$
$\displaystyle+\left\langle\lambda^{c}_{n+1}-\lambda^{b}_{n+1},(u_{n+1}^{c}-u_{n+1}^{b})^{-}\right\rangle=\frac{1}{k}\langle
u^{n}_{c}-u_{n}^{b},(u_{n+1}^{c}-u_{n+1}^{b})^{-}\rangle\leq 0.$

Note that $cu^{c}_{n+1}-bu^{b}_{n+1}\leq cu^{c}_{n+1}-cu^{b}_{n+1}$ for $c\leq
b$ and hence
$\langle\lambda^{c}_{n+1}-\lambda^{b}_{n+1},(u_{n+1}^{c}-u_{n+1}^{b})^{-}\rangle\geq
0$. We thus obtain that $u^{c}_{n+1}\geq u_{b}$.

∎

As a corollary of the monotonicity, we obtain the existence of $u(t)$ and it
solves (C.1). Uniqueness can be proved similarly as in [13].

## References

  * [1] U. M. Ascher and L. R. Petzold.  Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations.  Society for Industrial and Applied Mathematics, USA, 1st edition, 1998\.
  * [2] A. Berger, M. Ciment, and J. Rogers.  Numerical Solution of a Diffusion Consumption Problem with a Free Boundary.  SIAM Journal on Numerical Analysis., 12(4):646–672, 1975.
  * [3] H. Brezis and D. Kinderlehrer.  The smoothness of solutions to nonlinear variational inequalities.  Indiana University Mathematics Journal., 23(9):831–844, 1974.
  * [4] L. Bridge and B. Wetton.  A mixture formulation for numerical capturing of a two-phase/vapour interface in a porous medium.  Journal of Computational Physics, 225:2043–2068, 2007.
  * [5] L. A. Caffarelli.  The Obstacle Problem Revisited.  Journal of Fourier Analysis and Applications., 4(4):383–402, 1998\.
  * [6] J. Crank.  Free and Moving Boundary Problems.  Oxford University Press., Walton Street, Oxford., 1984.
  * [7] G. Dal Maso.  An Introduction to $\Gamma$-Convergence.  Progress in Nonlinear Differential Equations and their Applications, 8\. Birkhäuser Boston, Inc., Boston, MA, 1993.
  * [8] A. Fasano and M. Primicerio.  New results on some classical parabolic free-boundary problems.  Quarterly of Applied Mathematics., 38(4):439–460, 1981.
  * [9] A. Figalli, X. Ros-Oton, and J. Serra.  The singular set in the stefan problem.  arXiv, 2103.13379, 2021.
  * [10] M. Focardi, M. S. Gelli, and E. Spadaro.  Monotonicity formulas for obstacle problems with Lipschitz coefficients.  Calculus of Variations and Partial Differential Equations., 54(2):1547–1573, 2015.
  * [11] R. Glowinski and A. Wachs.  On the numerical simulation of viscoplastic fluid flow.  In R. Glowinski and J. Xu, editors, Numerical Methods for Non-Newtonian Fluids, volume 16 of Handbook of Numerical Analysis, pages 483–717. Elsevier, 2011.
  * [12] H. Huang, P. Lin, and W. Zhou.  Moisture transport and diffusive instability during bread baking.  SIAM Journal on Applied Mathematics, 68(1):222–238, 2007.
  * [13] K. Ito and K. Kunisch.  Parabolic variational inequalities: The Lagrange multiplier approach.  Journal de Mathématiques Pures et Appliquées., 85(3):415–449, 2006.
  * [14] T. Kärkkäinen, K. Kunisch, and P. Tarvainen.  Augmented lagrangian active set methods for obstacle problems.  Journal of Optimization Theory and Applications., 119(3):499–533, 2003.
  * [15] J. Kačur.  Method of Rothe in Evolution Equations. Teubner Verlagsgesellschaft., Leipzig, 1985.
  * [16] E. Lindgren and R. Monneau.  Pointwise regularity of the free boundary for the parabolic obstacle problem.  Calculus of Variations and Partial Differential Equations., 54(1):299–347, 2015.
  * [17] J.-L. Lions and G. Stampacchia.  Variational inequalities.  Comm. Pure Appl. Math., 20:439–519, 1967.
  * [18] E. Magenes.  Topics in Parabolic Equations: Some Typical Free Boundary Problems.  Boundary Value Problems for Linear Evolution Partial Differential Equations., 29:239–312, 1977.
  * [19] F. Miranda, J. F. Rodrigues, and L. Santos.  Evolutionary quasi-variational and variational inequalities with constraints on the derivatives.  Advances in Nonlinear Analysis., 9(1):250–277, 2020.
  * [20] S. Mitchell and M. Vynnycky.  The oxygen diffusion problem: Analysis and numerical solution.  Applied Mathematical Modelling., 39(9):2763–2776, 2015.
  * [21] R. Monneau.  On the Number of Singularities for the Obstacle Problem in Two Dimensions.  The Journal of Geometric Analysis., 13(2):359–389, 2003.
  * [22] L. Rubinstein.  The Stefan Problem.  Translations of. Math. Monog. Vol. 27, American Math. Society, Providence R.I., U.S.A., 1971.
  * [23] M. Rudd and K. Schmitt.  Variational Inequalities of Elliptic and Parabolic Type.  Taiwanese Journal of Mathematics., 6(3):287–322, 2002.
  * [24] D. G. Schaeffer.  A New Proof of the Infinite Differentiability of the Free Boundary in the Stefan Problem .  Journal of Differential Equations, 20:266–269, 1976.
  * [25] H. Singh and J. Hanna.  Pick-up and impact of flexible bodies.  Journal of the Mechanics and Physics of Solids, 106:46–59, 2017\.
  * [26] J. Stefan.  Ueber die Theorie der Eisbildung, insbesondere über die Eisbildung im Polarmeere.  Annalen der Physik., 278(2):269–286, 1890.
  * [27] M. Struwe.  Variational methods, volume 34 of A Series of Modern Surveys in Mathematics.  Springer-Verlag Berlin Heidelberg, 2008.
  * [28] R. K. Sundaram.  A First Course in Optimization Theory.  Cambridge University Press, 1996.
  * [29] J. W. Thomas.  Numerical Partial Differential Equations: Finite Difference Methods, volume 22 of Texts in Applied Mathematics.  Springer-Verlag, New York, 1995.
  * [30] C. Vuik.  Some historical notes about the stefan problem.  Nieuw Archief voor Wiskunde 4e serie, 11:157–167, 1993.
  * [31] G. S. Weiss.  A homogeneity improvement approach to the obstacle problem.  Inventiones Mathematicae., 138(1):23–50, 1999.