{'fraction_non_alphanumeric': 0.06516290726817042, 'fraction_numerical': 0.022877192982456142, 'mean_word_length': 3.5808229243203527, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the logarithmic type boundary modulus of continuity for the Stefan problem

A logarithmic type modulus of continuity is established for weak solutions to a two-phase Stefan problem, up to the parabolic boundary of a cylindrical space-time domain. For the Dirichlet problem, we merely assume that the spatial domain satisfies a measure density property, and the boundary datum has a logarithmic type modulus of continuity. For the Neumann problem, we assume that the lateral boundary is smooth, and the boundary datum is bounded. The proofs are measure theoretical in nature, exploiting De Giorgi's iteration and refining DiBenedetto's approach. Based on the sharp quantitative estimates, construction of continuous weak (physical) solutions is also indicated. The logarithmic type modulus of continuity has been conjectured to be optimal as a structural property for weak solutions to such partial differential equations.


INTRODUCTION
The classical Stefan problem aims to describe the evolution of the moving boundary between two phases of a material undergoing a phase change, for instance the melting of ice to water. The temperature of the material in both phases is governed by the heat equation, subject to the usual initial and boundary conditions. However, the interface of the two phases is free, and additional conditions need to be imposed on that part of the boundary (an unknown hyper-surface), which can be viewed as a kind of energy balance law. As such the classical Stefan problem is genuinely a parabolic free boundary problem, cf. [6]. The model accounts for a large class of physical phenomena, cf. [29].
The Stefan problem also admits a variational perspective, cf. [18,21,24] and [22, Chapter V, Section 9]. Under this point of view, we are led to formulate the initial-boundary value problem for the following nonlinear parabolic equation in a fixed space-time domain E T := E × (0, T ] ⊂ R N +1 , that is, Here β(·) represents the energy (more precisely, the enthalpy) of the two phases, i.e.
[u > 0] and [u < 0], and it is a maximal monotone graph that permits a jump when the temperature u is zero. The equation (1.1) is understood in the sense of differential inclusions. No explicit reference is made to any free boundary in this formulation. Instead, a mushy region, i.e. [u = 0], which consists of a mixture of the two states of the material, is allowed. The solutions are sought in Sobolev spaces and thus are certain weak ones. As a result, the question concerning the regularity of the weak solutions is naturally raised. This is the problem we try to tackle in this note. See [29, § 1.3] for a comparison between the classical and the weak formulations. See [19] for another formulation in terms of variational inequalities. Pioneering works of Caffarelli & Evans [4], DiBenedetto [9,10], Sacks [25], and Ziemer [31], have shown that while the energy admits a jump, the temperature is always continuous, even across the mushy region [u = 0]. Among these contributions, we single out [9,10] whose approach is closely followed in this work. In [9] the interior and up to the boundary continuity of the temperature was obtained, given general Neumann data or homogeneous Dirichlet data on the lateral boundary. However, the case of nonhomogeneous Dirichlet data turned out to be more involved and was solved later in [10]. An explicit modulus of continuity of the temperature, both in the interior and at the boundary, can be derived from the method used in [9,10]. Roughly speaking, it is of the following type (cf. [11]): (I) (0, 1) ∋ r → ω(r) = ln | ln(cr)| −σ for some c, σ > 0.
The method employed in [9,10] is flexible enough to deal with the situation when the Laplacian ∆ in (1.1) is replaced by a general quasilinear diffusion part and when a lower order term appears. Subsequently, the equation (1.1) replacing the Laplacian ∆ by the p-Laplacian ∆ p for p ≥ 2 has been studied in [1,2,26,27,28] by Urbano et al. In particular, we single out the recent work [1] where the interior modulus of continuity has been improved to be of the following type (II) (0, 1) ∋ r → ω(r) = | ln(cr)| −σ for some c, σ > 0.
Discarding a logarithm, this type of modulus of continuity represents an improvement even in the case of Laplacian, let alone the complication brought by the p-Laplacian. The modulus (II) has been conjectured to be optimal in [1]. (See also [5] for the one-phase problem.) Although a boundary modulus of continuity is derived in [2] for Dirichlet data, it is still of type (I).
How to achieve a boundary modulus of type (II) with Dirichlet data, and how to deal with the same issue for Neumann data, remain elusive. We will answer these questions in this note. The significance of continuity for weak solutions (temperatures) stems from their physical bearings; an explicit and sharp modulus has further mathematical implications, cf. [5]. Our quantitative estimates up to the boundary allow us to construct physical solutions with the type (II) modulus of continuity.
1.1. Statement of the results. Our main goal is to establish a modulus of continuity of type (II) for weak solutions to the Stefan problem, up to the parabolic boundary, given either Dirichlet or Neumann conditions. Moreover, for the Dirichlet problem, the p-Laplacian with p ≥ 2 is considered. Evidently, our results are new even for the Laplacian.
More precisely, denoting an open set in R N (N ≥ 1) by E and setting E T := E ×(0, T ], we are concerned with the following nonlinear parabolic equation (1.2) ∂ t β(u) − div A(x, t, u, Du) ∋ 0 weakly in E T .
Here β is a maximal monotone graph in R × R defined by for some constant ν > 0 that represents the exchange of latent heat.
The function A(x, t, u, ξ) : E T ×R N +1 → R N is assumed to be measurable with respect to (x, t) ∈ E T for all (u, ξ) ∈ R×R N , and continuous with respect to (u, ξ) for a.e. (x, t) ∈ E T . Moreover, we assume the structure conditions where C o and C 1 are given positive constants, and we take p ≥ 2. The set of parameters {ν, p, N, C o , C 1 } will be referred to as the (structural) data in the sequel. Before considering the initial-boundary value problems for (1.2), let us recall the parabolic p-Laplace type equation, which is pertinent to (1.2): The boundary datum satisfies (D) g ∈ L p 0, T ; W 1,p (E) , g continuous on E T with modulus of continuity ω g (·).
In addition, the initial datum satisfies Regarding the geometry of the boundary ∂E, we assume a measure density condition, that is, there exist α * ∈ (0, 1) and̺ ∈ (0, 1), such that for all x o ∈ ∂E, for every cube K ̺ (x o ) and 0 < ̺ ≤̺, there holds Here we have denoted by K ̺ (x o ) the cube of side length 2̺ and center x o , with faces parallel with the coordinate planes of R N . Intuitively, the condition (G) means that one can place a cone at x o exterior to E, with an angle quantified by α * .
Throughout the rest of this note, we will use the symbols to denote (backward) cylinders with the indicated positive parameters; we omit the vertex (x o , t o ) from the notations for simplicity. We will also denote the lateral boundary by S T := ∂E × (0, T ] and the parabolic boundary by ∂ P E T : The formal definition of solution to the Dirichlet problem (1.5) will be postponed to Section 1.2. Now we first present the regularity theorem concerning the Dirichlet problem (1.5) up to the parabolic boundary ∂ P E T . Theorem 1.1. Let u be a bounded weak solution to the Dirichlet problem (1.5) under the condition (1.3) with p ≥ 2. Assume that (D), (U o ) and (G) hold. Then u is continuous in E T . More precisely, there is a modulus of continuity ω(·), determined by the data, α * ,̺, u ∞,E T , ω o (·) and ω g (·), such that for every pair of points (x 1 , t 1 ), (x 2 , t 2 ) ∈ E T . In particular, there exists σ ∈ (0, 1) depending only on the data and α * , such that if where C g , C uo > 0 and λ > σ, then the modulus of continuity is with some C > 0 depending on the data, λ, α * , u ∞,E T , C uo and C g .
Now we consider the Neumann problem. In order to deal with possible variational data on S T , we assume ∂E is of class C 1 , such that the outward unit normal, which we denote by n, is defined on ∂E pointwise. Let us consider the initial-boundary value problem of Neumann type: where the structure conditions (1.3) and the initial condition (U o ) are retained. On the Neumann datum ψ we assume for simplicity that, for some absolute constant C 2 , there holds More general conditions should also work (cf. Section 2, Chapter II, [7] for every pair of points (x 1 , t 1 ), (x 2 , t 2 ) ∈ E T . In particular, there exists σ ∈ (0, 1) depending only on the data, C 2 and the structure of ∂E, such that if where C uo > 0 and λ > σ, then the modulus of continuity is for all r ∈ (0, 1), with some C > 0 depending on the data, the structure of ∂E, C 2 , u ∞,E T and C uo . 2) with lower order terms: the convection resulting from the heat transfer is reflected by a lower order term. Also the function β(u) could present more general forms and reflect thermal properties of the material that may slightly change according to the temperature, as long as it permits only one jump. (For the case of multiple jumps, see [17].) The modifications of the proofs can be modeled on the arguments in [7,9,10]. However, we will not pursue generality in this direction, instead focus will be made on actual novelties.
which do not assume any knowledge on the time derivative. In contrast to (1.7), a special character for the notion of solution to the Stefan problem (1.2) is that we require solutions to possess a time derivative in the Sobolev sense. This is necessary to justify the calculations in the proofs of the theorems. On the other hand, we will explain in Section 7 how to construct continuous weak solutions to (1.2) without any knowledge on the time derivative.
for all non-negative test functions and hence all the integrals are well-defined. Note that the above integral formulation itself does not involve the time derivative of u. On the other hand, if we use the time derivative of u, the integral formulation may be written as where ζ is as in (1.8) and ν(x, t) ≥ 0 is given by and χ is the characteristic function of the indicated set. A function u that is both a local weak sub-solution and a local weak super-solution to (1.2) is a local weak solution.

Notion of Solution to the Dirichlet Problem. A function
for all non-negative test functions ζ satisfying (1.8). An equivalent form can be given as in (1.9) involving ∂ t u. Clearly, now the compact set K in (1.8) could touch the boundary ∂E, i.e. K ⊂ E. Moreover, the initial datum is taken in the sense that there exists a selection v o ⊂ β(u o ) and for any compact set K ⊂ E, The Dirichlet datum g is attained under u ≤ (≥)g on ∂E in the sense that the traces of (u− g) ± vanish as functions in W 1,p (E) for a.e. t ∈ (0, T ], i.e. (u−g) ± ∈ L p (0, T ; W 1,p o (E)). Notice that no a priori information is needed on the smoothness of ∂E for the time being.
A function u that is both a weak sub-solution and a weak super-solution to (1.5) is a weak solution.

Notion of Solution to the Neumann Problem. A function
is a weak sub(super)-solution to (1.6), if there is a selection v ⊂ β(u), for every compact set K ⊂ R N and every sub- ψ(x, t, u)ζ dσdt for all non-negative test functions ζ satisfying (1.8), where K is an arbitrary compact set in R N now, and dσ denotes the surface measure on ∂E. The Neumann datum ψ is reflected in the boundary integral on the right-hand side. Moreover, the initial datum is taken as in the Dirichlet problem. An equivalent form can be given as in (1.9) involving ∂ t u.
A function u that is both a weak sub-solution and a weak super-solution to (1.6) is a weak solution.
1.3. DiBenedetto's approach: revisit and refinement. Any known approach to the local continuity issue for the equation (1.1) relies on, in one way or another, De Giorgi's method. Since there are a lot of technicalities to follow that might obscure the main ideas behind them, we briefly review the original approach of DiBenedetto in [9,10] for p = 2, and highlight the main improvement to the method in order to achieve a type (II) modulus of continuity.
For a cylinder we introduce the numbers µ ± and ω satisfying The goal is to reduce the oscillation of u over a cylinder smaller than Q ̺ with the same vertex (x o , t o ). To this end, we first observe that one of following must hold: then u ≥ µ − + 1 8 ω a.e. in Q 1 2 ̺ , which in turn yields a reduction of oscillation ess osc The singularity of β(·) at [u = 0] is reflected by the dependence on ω of α in (1.11). The second step is to consider the case when (1.11) does not hold. Since µ + − 1 4 ω ≥ µ − + 1 4 ω always holds, the reverse of the measure information (1.11) implies that The key to proceeding now is to observe that because of (1.10) 1 , the function v defined by v := 1 4 ω − (u − k) + with k = µ + − 1 4 ω > 0, is actually a non-negative, weak super-solution to the parabolic equation (1.4) with p = 2 in Q ̺ , as the set [u > k] excludes potential singularity of β(u). The information in (1.12) can be rephrased in terms of v: As such the machinery of De Giorgi can be applied to v like in [22,Chapter II,Section 7] and the measure information (1.12) translates into pointwise positivity for v (thus also for µ + − u) up to the top of the cylinder Q ̺ . This kind of property for non-negative, weak super-solutions is called expansion of positivity. More precisely, there holds µ + − u ≥ η(α)ω a.e. in Q 1 2 ̺ (α), which gives the reduction of oscillation ess osc The singularity of β(·), reflected in the dependence of α on ω, has now been passed from the measure information (1.12) to this reduction of oscillation, through the dependence of η on α and hence on ω. Such dependence of η on α is traced by η ≈ 2 − 1 α q , for a generic positive constant q determined by the data; see also [15,Chapter 4,Section 2]. Apart from technical complications, these are the main steps in [9].
Strikingly, the dependence of η on α can be ameliorated in the sense that η ≈ α q , for a generic positive constant q determined by the data. This was first discovered by DiBenedetto & Trudinger in [16] for the elliptic De Giorgi class via a Krylov-Safonov type covering argument. Various parabolic versions have been developed since then; see for instance [15,23,30].
It is exactly the improvement of the dependence of η on α in the expansion of positivity that leads to the refinement of the modulus of continuity from (I) to (II). Indeed, iterating the arguments presented above yields the reduction of oscillation ess osc Qn u ≤ ω n , n = 0, 1, · · · , along a nested family of cylinders {Q n } with their common vertex at (x o , t o ). The sequence {ω n } obeys two types of recurrences, contingent upon the dependence of η on α, recalling also the dependence of α on ω in (1.11): n ω n , (Type I), ω n+1 = 1 − ω q n ω n , (Type II).
Here we have again used q > 0 as a generic constant determined by the data. An inspection on the sequence {ω n } reveals that for n large ω n (ln n) −σ , (Type I), for some proper σ ∈ (0, 1 q ). If the size of Q r ≈ Q n is quantified in a geometric fashion, r ≈ ( 1 2 ) n for instance, then we have n ≈ ln 1 r . Consequently, the modulus of continuity can be estimated by The type (II) interior modulus of continuity has been achieved in [1]. Instead of the expansion of positivity, a weak Harnack inequality was employed in [1]. One purpose of this note is to demonstrate that the expansion of positivity is more flexible when the type (II) modulus of continuity is requested not only in the interior, but also at the boundary, given Neumann boundary data. Needless to say, the expansion of positivity lies at the heart of any kind of Harnack's estimates, cf. [13,15].
In the above outline, the singularity of β(·) at [u = 0] is accommodated by choosing proper levels k for u, either close to µ − or to µ + . However, when the boundary regularity is considered, given general Dirichlet data, the level k has to obey extra restrictions, such that (u − k) ± vanish on the lateral boundary, cf. (2.1). As a result, the analysis becomes more involved. The most crucial part of DiBenedetto's approach in [10] consists in quantifying the largeness of |[u ≤ 0]| where the singularity appears. If it is small enough, then the usual parabolic scaling prevails. Otherwise, the largeness of |[u ≤ 0]| is incorporated into the singular term on left-hand side of the energy estimate, and then balances off the relative largeness brought by the singularity on the right-hand side; in this case, the usual parabolic scaling fails and intrinsically scaled cylinders have to be introduced. Retrospectively, this idea lies at the origin of the method of intrinsic scaling, cf. [7,15,28].
We will adapt this idea, in further balancing the singularity of β(·) and the degeneracy brought by ∆ p , and in refining the approach by an argument from [14,Section 4.4] and [12, Proposition 5.1], which allows us to shrink the measure of the set where u ≈ 0 in a faster fashion than the original De Giorgi method, cf. Lemma 3.2. We will not try to recapitulate all the steps here beforehand, as they are more involved and delicate than the interior case. Instead we feel it is more appropriate to add remarks, which may assist understanding of the technicalities, along the course of the proof.

Acknowledgement.
This research has been funded by the FWF-Project P31956-N32 "Doubly nonlinear evolution equations".

THE DIRICHLET PROBLEM: PRELIMINARY BOUNDARY ESTIMATES
We will deal with the oscillation decay of u in the vicinity of the lateral boundary in Sections 2 -3; in the interior in Section 4; at the initial level in Section 5. Afterwards the proof of Theorem 1.1 follows from standard covering arguments based on these three cases in Section 5. In the present section, we derive energy estimates for sub(super)-solutions to the Dirichlet problem (1.5), near the lateral boundary. Consequently, two De Giorgi type lemmas are implied by the energy estimates.
2.1. Energy estimates. In this section, we derive some energy estimates near S T , for solutions to the Dirichlet problem (1.5). We will work within the cylinder Q R,S := In order to employ the truncated functions (u − k) ± across the lateral boundary, we need to impose certain restrictions on the level k, i.e., Let ζ be a non-negative, piecewise smooth cutoff function vanishing on In such a way, the test functions In what follows, when we deal with integrals of (u − k) ± across the lateral boundary S T , the functions (u − k) ± will be understood as zero outside E T . Likewise, when we deal with a sub(super)-solution u across the lateral boundary S T , we tacitly understand it as for k satisfying (2.1). Upon using the above test function ϕ in the boundary version of (1.9), we have the following energy estimate; the detailed calculation is analogous to the one for [9, (2.7)].
Based on the energy estimate in Proposition 2.1, we may derive two tailored versions. Note that the functions Φ ± carry information near the singularity of β(·). The first tailored version examines the case when Φ ± = 0, and hence singularity disappears and it recovers the energy estimate for the parabolic p-Laplacian (1.4).
Proposition 2.2. Let the hypotheses in Proposition 2.1 hold. If, in addition, k ≥ 0 in the case of sub-solutions or k ≤ 0 in the case of super-solutions, then for every non-negative, piecewise smooth cutoff function ζ vanishing on ∂ P Q R,S , there holds ess sup Proof. We need only to notice that under the current restrictions on k, the functions Φ ± vanish. Also now the cutoff function ζ is required to vanish at the bottom level t o − S.
The second one examines the effect of the singularity of β(·) more carefully.
Proof. Let us omit the reference to x o for simplicity. We only need to calculate Then substituting it back to the energy estimate in Proposition 2.1, we conclude.
as in Section 2.1, we introduce the numbers µ ± and ω satisfying Let θ ≥ θ > 0 be parameters to be determined. The cylinders Q ̺ (θ) and Q ̺ ( θ) are coaxial with Q R,S and with the same vertex (x o , t o ); moreover, we will impose that We now present the first De Giorgi type lemma.
There exists a constant c o ∈ (0, 1) depending only on the data, such that if Proof. Let us assume (x o , t o ) = (0, 0). In order to use the energy estimates in Proposi- Recall that Q ̺n (θ) = K n ×(−θ̺ p n , 0) and Q̺ n (θ) = K n ×(−θ̺ p n , 0). Introduce the cutoff function 0 ≤ ζ ≤ 1 vanishing on the parabolic boundary of Q n and equal to identity in Q n , such that |Dζ| ≤ γ 2 n ̺ and |∂ t ζ| ≤ γ 2 pn θ̺ p . In this setting, the energy estimates in Propositions 2.2 -2.3 combined will yield that Now setting 0 ≤ φ ≤ 1 to be a cutoff function which vanishes on the parabolic boundary of Q n and equals the identity in Q n+1 , an application of the Hölder inequality and the Sobolev imbedding [7, Chapter I, Proposition 3.1] gives that In the second-to-last line we used the above energy estimate. In terms of Y n = |A n |/|Q n |, this can be rewritten as N+p p , which amounts to Since Y n → 0 as n → ∞, we have This concludes the proof of the lemma.
The second De Giorgi type lemma is as follows.
Lemma 2.2. Let the hypotheses in Lemma 2.1 hold. There exists ξ ∈ (0, 1) depending only on the data and α * , such that if To this end, let us first estimate the last integral on the right-hand side of the energy estimate via the given measure theoretical information: This means that the integral can be combined with a similar term involving ∂ t ζ p on the right-hand side of the energy estimate in Proposition 2.3. Consequently, we obtain an energy estimate of (u − k) − from which the theory of parabolic p-Laplacian applies, cf. Lemmas 3.4 -3.5 and Appendix A. Therefore we obtain a constant ξ determined only by the data and α * of (G), such that . The proof is now completed.
Remark 2.1. Before heading to the reduction of boundary oscillation, we try to shed some light on the roles of Lemmas 2.1 -2.2 as a preparation for the technically involved Section 3. Both lemmas are examining the measure information needed in order for the degeneracy of ∆ p to dominate the singularity of β(·). As such the second term on the left-hand side of the energy estimate in Proposition 2.3 is discarded and the time scaling θ = (ξω) 2−p is employed to accommodate the degeneracy of ∆ p ; the reduction of oscillation in this case will be demonstrated in Section 3.1. Attention is called to the difference between Section 3.1 and Section 3.4, where the same time scaling is evoked. Indeed, pointwise information is imposed in Section 3.4 to avoid the singularity of β(·) completely, such that the case falls into the scope of the theory for the parabolic p-Laplacian in [7]. On the other hand, violation of the measure information in Lemmas 2.1 -2.2 means that the singularity of β(·) overtakes the degeneracy of ∆ p and will yield an estimate from below on |[u ≤ 0]| (see (3.13)), which bears the singularity of β(·). This lower bound of |[u ≤ 0]| in turn gives a lower bound for the second term on the left-hand side of the energy estimate in Proposition 2.3; in this case the first term on the left-hand side of the energy estimate in Proposition 2.3 is discarded (cf. Lemma 3.1) and consequently, the "longer" time scaling θ = (δξω) 1−p will be needed to accommodate the singularity, cf. Remark 3.2. The reduction of oscillation in this case will be taken on in Section 3.2.

THE DIRICHLET PROBLEM: REDUCTION OF BOUNDARY OSCILLATION
Let u be a solution to the Dirichlet problem (1.5). Fix (x o , t o ) ∈ S T and define the numbers µ ± and ω as in Section 2.2 with Q R,S replaced by , where̺ is defined in (G), which we may assume to be so small that t o ≥ ( 1 8̺ ) p−1 . We now state the main result of this section, which describes the oscillation decay of u, along shrinking cylinders with the same vertex (x o , t o ) fixed on the lateral boundary S T . g ≤ C g | ln r| λ for all r ∈ (0,̺), with C g > 0 and λ > σ, then for some γ > 1 depending only on the data, α * and C g , ess osc The proof of Proposition 3.1 will be expanded on in Sections 3.1 -3.5. To start with, we introduce the following intrinsic cylinders where δ =ξω q for some positiveξ and q, which together with the positive ξ are to be determined in terms of the data and α * . We may assume with no loss of generality that since otherwise we would have Moreover, one of the following must hold: Let us take for instance the first one, i.e. (3.3) 1 , as the other one can be treated analogously.
Next, we observe that one of the following two cases must hold: Otherwise, the violation of both inequalities implies that If (3.4) 1 holds, then thanks to (3.3) 1 , the function (u − k) + with k = µ + − 2 −n ω will satisfy the energy estimate in Proposition 2.2 for all n ≥ 2. In this case, the reduction of oscillation is achieved in Section 3.4, cf. Lemma 3.4. Consequently we need only to examine the case when (3.4) 2 holds. This case splits into two sub-cases. In Sections 3.1 -3.3 we deal with the situation when |µ − | is near zero; whereas the opposite situation is treated in Section 3.4, cf. Lemma 3.5.
In summary, Sections 3.1 -3.3 will be devoted to reducing the oscillation of u near the set [u ≈ 0], where the singularity of β(u) appears; Section 3.4 treats the cases when the singularity of β(u) is avoided, and hence the behavior of u resembles the one to the parabolic p-Laplacian (1.4) and the theory in [7] applies.
3.1. Reduction of oscillation near zero. Let us deal with the case when (3.4) 2 holds. Throughout Sections 3.1 -3.3, we always assume that for some positive parameters {ξ,ξ, q} to be determined in terms of the data and α * only. When this restriction (3.6) does not hold, the case will be examined in Section 3.4. In Sections 3.1 -3.3, we will work with u as a super-solution and reduce its oscillation near its infimum µ − in a smaller, quantified cylinder with the same vertex (x o , t o ).

3.2.
Reduction of oscillation near zero continued. In this section, we continue to examine the situation when the measure condition in Lemma 2.1 is violated: N+p p |Q 1 2 ̺ (θ)|, and when the condition in Lemma 2.2 is violated: Combining (3.11) and (3.12), we may obtain that for all r ∈ [2̺, 8̺], To proceed, letδ ∈ (δ, 1) be a free parameter and setθ := (δξω) 1−p . Recall that θ = (δξω) 1−p , θ = (ξω) 2−p , and we have assumed θ(8̺) p < ̺ p−1 in (3.1). Therefore Observe also that for any t ∈ [t o −θr p , t o ] and anyδ ∈ (δ, 1), there holds Combining (3.14) and (3.15) gives that ess sup The above analysis together with Proposition 2.3 yields the following energy estimate. and setting k = µ − +δξω with δ ∈ (δ, 1), there exists a positive constant γ depending only on the data, such that for all σ ∈ (0, 1) and all r ∈ [2̺, 8̺] we have Remark 3.1. Note that the restriction of |µ − | in the last line of Lemma 3.1 has been imposed to ensure that the level k ≥ 0 as required in Proposition 2.3. Recall that such a restriction has been always assumed in this section, cf. (3.6). Note carefully that the numberδ ∈ (δ, 1) in Lemma 3.1 is a free parameter, whereas δ is yet to be determined.
Remark 3.2. When we generate the energy estimate of Lemma 3.1, the first term of the energy estimate in Proposition 2.3 has been discarded. Hence the singular effect of β(·) overtakes the degeneracy brought by the p-Laplacian in Lemma 3.1. The machinery of De Giorgi has to be reproduced to reflect this change. This is the reason why the time scaling will be changed from θ = (ξω) 2−p to θ = (δξω) 1−p in the following two lemmas.
Based on the energy estimate in Lemma 3.1, we then show that the measure density of the set [u ≈ µ − ] can be made arbitrarily small in a way that favors a later application for our purpose. The geometric condition (G) of the boundary ∂E comes into play.
Lemma 3.2. Suppose the hypotheses in Lemma 3.1 hold. There exists a positive constant j * depending only on the data and α * , such that for any m ∈ N we have either Proof. We employ the energy estimate of Lemma 3.1 in Q 2̺ ( θ) with σ = 1 2 and the levels k j := µ − + ξω 2 j , j = 0, 1, · · · , j * . Restricting |µ − | ≤ 2 −j * ξω implies that k j ≥ 0 for all j = 0, 1, · · · , j * . Therefore, assuming j * has been chosen and m is fixed for now, and taking into account the definition of θ(j * , m), the energy estimate in Lemma 3.1 yields thaẗ where A j,2̺ := u < k j ∩ Q 2̺ ( θ). Next, we apply [7, Chapter I, Lemma 2.2] slicewise to u(·, t) for t ∈ (− θ̺ p , 0] over the cube K ̺ , for levels k j+1 < k j . Taking into account the measure theoretical information derived from (G):

N. LIAO
Here we used in the last line the notation A j,̺ (t) := u(·, t) < k j ∩ K ̺ . We now integrate the last inequality with respect to t over (− θ̺ p , 0] and apply Hölder's inequality in time. With the abbreviation A j,̺ = [u < k j ] ∩ Q ̺ ( θ) this procedure leads to Here the constant γ in the last line has absorbed α * . Now take the power p p−1 on both sides of the above inequality to obtain Add these inequalities from 0 to j * − 1 to obtain Similarly, replacing ̺ by 2̺, we have The constant γ(data, α * ) appearing in the above two inequalities may be different, but we can take the larger one. Suppose j * has been chosen for the moment. We may repeat the above arguments for the same k j with j = j * , · · · , 2j * . Restricting |µ − | ≤ 2 −2j * ξω implies k j ≥ 0. Hence the energy estimate can be written in the same form with such choices of k j , and the geometry condition (G) will allow us to apply [7, Chapter I, Lemma 2.2] just like above. Consequently this will lead us to Here the constant γ is the same as above.
Combining the above estimates would yield that Iterating the procedure m times would give us that Finally we may choose j * so large that The proof is completed.
Next we apply the energy estimate in Lemma 3.1 to show a De Giorgi type lemma. Notice that the time scaling used here is different from the one in Lemmas 2.1 -2.2; recall Remark 2.1 and Remark 3.2 for an explanation.  ∈ (0, 1). There exists a constant c 1 ∈ (0, 1) depending only on the data, such that if Proof. In order to use the energy estimate in Lemma 3.1, we set as in (2.2) for n = 0, 1, · · · , We will use the energy estimate in Lemma 3.1 with the pair of cylinders Q n ⊂ Q n . Now the constantδ in Lemma 3.1 is replaced by (1 + 2 −n )δ here as indicated in the definition of k n . In this setting, enforcing |µ − | ≤ 1 2 δξω, the energy estimate in Lemma 3.1 may be rephrased as A n = u < k n ∩ Q n . Now setting 0 ≤ φ ≤ 1 to be a cutoff function which vanishes on the parabolic boundary of Q n and equals the identity in Q n+1 , an application of the Hölder inequality and the Sobolev imbedding [7, Chapter I, Proposition 3.1] gives that δξω In the second-to-last line we used the above energy estimate. In terms of Y n = |A n |/|Q n |, this can be rewritten as for a constant γ depending only on the data and with C = C (N, p). Hence, by [7, Chapter I, Lemma 4.1], there exists a positive constant c 1 depending only on the data, such that Y n → 0 if we require that Y o ≤ c 1 (ξω) b , which amounts to Since Y n → 0 as n → ∞, we have This concludes the proof of the lemma.

3.3.
Reduction of oscillation near zero concluded. Under the conditions (3.3) 1 , (3.4) 2 , (3.11) and (3.12) 1 , we may reduce the oscillation in the following way. First of all, let ξ be determined in Lemma 2.2 and j * be determined in Lemma 3.2 in terms of the data and α * . We may assume without loss of generality that Lemma 3.2 holds true for Q 4̺ ( θ) instead of Q ̺ ( θ). Then we choose the integer m so large to satisfy that where c 1 is the constant appearing in Lemma 3.3. Next, we can fix 2δ = 2 −mj * in Lemma 3.3. Consequently, by the choice of m in (3.16), Lemma 3.3 can be applied assuming that |µ − | ≤ 1 2 δξω, and we arrive at All the constants {ξ,ξ, η, q} have been determined in terms of the data and α * .

3.4.
Reduction of oscillation away from zero. In this section, we handle the case when the equation (1.2) resembles the parabolic p- Laplacian (1.4). The reduction of oscillation is summarized in the following two lemmas; the proofs rely on the theory for the parabolic p-Laplacian (1.4) in [7, Chapter III, Section 12], which we collect in Appendix A for completeness.
Under the pointwise conditions of Lemmas 3.4 -3.5, the singularity of β(·) is avoided, and the energy estimates of Proposition A.1 are verified for (u − k) + and for (u − k) − , with ω there replaced by 1 4 ω and 1 2 δξω respectively.

3.5.
Derivation of the modulus of continuity. This is the final part of the proof of Proposition 3.1. Recall that we have defined the quantities µ ± and ω over for some ̺ ∈ (0, 1 8̺ ), and we have performed the arguments in Sections 3.1 -3.4 starting from the intrinsic relation (3.1) 2 , that is, The positive constants {ξ,ξ, q} have been selected in terms of the data and α * , in the course of proof in Sections 3.1 -3.4. Therefore, combining all the arguments in Sections 3.1 -3.4, we arrive at the reduction of oscillation (3.18) ess osc where θ = (ξω) 2−p and {ξ,ξ, q, α o , η, A} have been properly determined in terms of the data and α * . In order to iterate the arguments of Sections 3.1 -3.4, we set We need to choose ̺ 1 > 0, such that the following two set inclusions hold: Just like (3.1), they will validate the intrinsic relation needed to repeat the arguments of Sections 3.1 -3.4. To this end, we may assume that ηω q ≤ 1 2 and then observe that Hence to satisfy the first set inclusion, it suffices to choose ̺ 1 by setting i.e. ̺ p 1 = ξωq̺ p , whereq := 1 + (p − 1)q and ξ depends only on the data and α * . It is not hard to verify the second inclusion also with such a choice of ̺ 1 . As a result of the set inclusions and (3.18), we obtain the following intrinsic relation ess osc which takes the place of (3.1) 2 in the next stage. Therefore the arguments of Sections 3.1 -3.4 can be repeated within Q 8̺ 1 ( θ 1 ), to obtain that either ess osc . Here the constants {ξ,ξ, q, α o , η, A} are the same as in the first step. Now we may proceed by induction and introduce the following notations: . Upon using the induction we have that for all n = 0, 1, · · · , (3.19) Q ′ n+1 ⊂ Q n and ess osc Qn u ≤ ω n+1 .
Next we derive an explicit modulus of continuity of u inherent in (3.19). To this end, we first define a sequence {r n } and cylinders { Q n } by Notice that ̺ n ≤ r n and Q n ⊂ Q n for all n ≥ 0. Observe that if there exists n o ∈ N, such that a sequence {a n } satisfies a n+1 ≥ max a n (1 − ηa q n ), Ar αo n , 2 osc Qn∩S T g for all n ≥ n o and a no ≥ ω no , then a n ≥ ω n for all n ≥ n o . As a matter of fact, we may take a n = (1 + n) −σ max{1, ω} with σ ∈ (0, 1 q ). In such a case, the number n o is determined by the parameters {η, q, σ, α o , A, ξ, C g , λ}. Indeed, the following two assertions are verified by elementary computations: a n+1 ≥ a n (1 − ηa q n ), a n+1 ≥ Ar β n , for all n ≥ n o and n o is determined by {η, q, σ, α o , A}. Using the modulus of continuity of g, there exists γ = γ( ξ, C g , λ), such that Recall that λ > σ. Hence n o here can be determined in terms of { ξ, C g , λ} to ensure that for any n ≥ n o , a n+1 ≥ osc Once the number n o has been determined by the parameters {η, q, σ, α o , A, ξ, C g , λ} and independent of ω, we may assume that n o = 0 for simplicity. As a result of a o ≥ ω o , we arrive at ess osc Qn u ≤ ω n+1 ≤ max{1, ω} (n + 2) σ for all n = 0, 1, · · · .
Let us take a number 4r ∈ (0, ̺). There must be some n such that The right-hand side inequality will yield that Q r (ω 2−p ) ⊂ Q 1 4 ̺n (θ n ) = Q n , and consequently (3.20) ess osc Next we analyze the left-hand side inequality. First notice that by the definition of ω n and assuming ηω q n ≤ 1 2 , we can estimate Use this to estimate ̺ p n+1 from below by Then taking logarithim on both sides, a simple calculation will give us a lower bound of n by n 2 ≥ γ 1 − ln min{1, ω} ln ̺ r , for some γ = γ(p, q, ξ).
Substituting this into (3.20), we obtain the modulus of continuity ess osc for all r ∈ (0, ̺).
This concludes the proof of Proposition 3.1

REDUCTION OF INTERIOR OSCILLATION
For a cylinder 1), we introduce the numbers µ ± and ω satisfying The main result of this section is the following.

Preliminaries.
We first collect a set of lemmas from which the proof of Proposition 4.1 will follow. Let us remark that the energy estimate in Proposition 2.1 still holds true in the interior with arbitrary k ∈ R. Keeping this remark in mind, we present the following De Giorgi type lemma similar to Lemma 2.1. Let θ > 0 be a parameter to be determined. The cylinder Q ̺ (θ) is coaxial with Q o and with the same vertex (x o , t o ); moreover, we will impose that Then we have the following.
Proof. This is exactly Lemma 2.1 in the interior. The proof runs almost verbatim.
Therefore, in any case the term Φ − vanishes in the energy estimate of Proposition 2.1, whereas the last spatial integral at t o also vanishes due to the initial information. As a result, we obtain energy estimates similar to that of the parabolic p- Laplacian (1.4), and the proof can be completed as in [15,Chapter 3,Section 4]. To this end, we only need to impose the additional condition that such that the levels k n become admissible in the application of Proposition 2.1.
The next lemma asserts that if we truncate a solution by a positive number, the resultant function is a sub-solution to the parabolic p- Laplacian (1.4).
A careful inspection of the proof there actually reveals that δ and b can be selected independent of α. Hence we have absorbed them into η here. Moreover, Lemma 4.4 is stable as p → 2.

4.2.
Proof of Proposition 4.1. Assume (x o , t o ) = (0, 0). Let θ = L( 1 4 ω) 2−p with L = L o ω −a for some a, L o > 0 to determined in terms of the data. To proceed we will assume that otherwise we would have As in (3.3) we may also assume that µ + − 1 4 ω > 1 4 ω > 0, and hence Lemma 4.4 is at our disposal.
Let us suppose L is a large integer with no loss of generality, and slice Q ̺ (θ) into a vertical stack of L coaxial cylinders all congruent with Q ̺ (θ). Here we have setθ := ( 1 4 ω) 2−p for ease of notation. Let us focus on the bottom cylinder Q L−1 . The proof unfolds along the following two alternatives: Let us first suppose (4.2) 1 holds true. An application of Lemma 4.1 would give us that u ≥ µ − + 1 8 ω a.e. in 1 2 Q L−1 . Here the notation 1 2 Q L−1 should be self-explanatory in view of Lemma 4.1. In particular, the above pointwise lower bound of u holds at the time level t o = −(L − 1)θ̺ p and it serves as the initial datum for an application of Lemma 4.2. As a result, by choosing ξ(ω) so small that γ o (ξω) 2−p ≥ L( 1 4 ω) 2−p , we obtain that (4.3) u ≥ µ − + 1 2 ξω a.e. in K 1 2 ̺ × (t o , 0). Note that the choice of ξ(ω) is made out of the data and L = L o ω −a , that is, where a, L o > 0 will be determined later in terms of the data.
Next, we turn our attention to the second alternative (4.2) 2 . Since µ + − 1 4 ω ≥ µ − + 1 4 ω, we may rephrase (4.2) 2 as Then it is not hard to see that there exists Based on this measure theoretical information at t * , an application of Lemma 4.4 yields that (4.4) µ + − u(·, t) ≥ η o α κ ω a.e. in K 1 2 ̺ for all times Now we may determine L by setting Hence we have chosen a := N +p p (p − 2)κ and L o := (4η o c κ o 4 − N+p p κ ) 2−p . Combining (4.3) and (4.4), and taking (4.1) into account, we obtain the reduction of oscillation ess osc where q = N +p p κ and η depends only on the data. In order to iterate the above arguments, we define we need to choose ̺ 1 > 0, such that Since we may assume ω 1 ≥ 1 2 ω and estimate , the first set inclusion amounts to requiring that i.e. ̺ 1 = c̺ for some c ∈ (0, 1) determined by the data only. The second set inclusion is verified similarly for such a choice of ̺ 1 . As a result, we arrive at ess osc Now we may proceed by induction and introduce the following notations: Upon using the induction we have that for all n = 0, 1, · · · , ess osc Qn u ≤ ω n .
From now on the deduction of a logarithmic type modulus of continuity is performed as in Section 3.5.

PROOF OF THEOREM 1.1
Before starting the proof of Theorem 1.1, we still have to deal with the continuity at the initial level. This is considerably simpler than the lateral boundary case and the interior case. Let us define µ ± and ω to be the supreme/infimum and the oscillation of u over the forward cylinder Q + for some x o ∈ E and ̺ ∈ (0, 1). We will impose that ω 2−p ̺ p ≤ ̺ p−1 , such that Then we have the following oscillation decay estimate.
Proposition 5.1. Suppose the hypotheses in Theorem 1.1 hold true. There exist γ > 1 and α 1 ∈ (0, 1) depending only on the data, such that ess osc Proof. Observe first that one of the following two alternatives must hold: otherwise we obtain With no loss of generality, let us suppose the second alternative holds true. Then we are in the position to apply Lemma 4.2 (recalling Remark 4.1), such that for some γ o ∈ (0, 1) depending only on the data, there holds This in turn yields a reduction of oscillation ess osc Notice that the above reduction of oscillation at the initial level bears no singularity of β(·) and in fact yields a Hölder type decay. As such the final oscillation decay is dominated by that of u o and g near the initial level. The technical realization runs similar to Section 3.5; see also [3, Section 7.1]. We refrain from giving further details, to avoid repetition.
The proof of Theorem 1.1 consists of several cases, which rely on Proposition 3.1 near the lateral boundary, Proposition 4.1 in the interior and Proposition 5.1 near the initial level.
Let us consider two points z i := (x i , t i ) in E T with i = 1, 2, satisfying where̺ is defined in the geometric property (G) of ∂E. Suppose first that dist p (z i , ∂ P E T ) := inf{dist p (z i , Q) : Q ∈ ∂ P E T } > 1 2̺ for i = 1, 2, then according to the interior estimate of Proposition 4.1, we have for some proper σ ∈ (0, 1) that where C > 0 takes into account the data and u ∞,E T . Next suppose that one of {z 1 , z 2 }, z 1 for instance, satisfies dist p (z 1 , ∂ P E T ) ≤ 1 2̺ . Suppose further that t 1 , t 2 > ( 1 2̺ ) p−1 . Then we may apply Proposition 3.1 to obtain the same type of estimate as (5.1). The present constant C involves the data, α * , λ, C g and u ∞,E T .
Finally if one of {t 1 , t 2 }, t 1 for instance, satisfies that t 1 ≤ ( 1 2̺ ) p−1 , then we may apply Proposition 5.1 to obtain the same type of estimate as (5.1). The present constant C involves the data, λ, C g , C uo and u ∞,E T .

THE NEUMANN PROBLEM
In this section, we will focus on the Neumann problem (1.6). The treatment more or less parallels the interior regularity in Section 4. However, we are not able to treat a general p-Laplacian for the moment; see Remark 6.2.
First of all, an energy estimate is in order.  1 and (N) holds. There exists a constant γ > 0 depending on C o , C 1 and the structure of ∂E, such that for all cylinders , every k ∈ R, and every non-negative, piecewise smooth cutoff function ζ vanishing on and ν(x, t) is a selection out of [0, ν] for u(x, t) = 0, and ν(x, t) = ν for u(x, t) < 0.
Proof. We deal with sub-solutions only as the other case is analogous. Upon using the test function (u − k) + ζ 2 in the weak formulation of Section In order to treat the boundary integral, we make use of (N), apply the trace inequality (cf. [ Substituting it back yields the desired energy estimate. Notice also that the structure of ∂E enters the constant γ through the trace inequality. Remark 6.1. The vertex (x o , t o ) of the cylinder Q R,S needs not to be on the lateral boundary S T . If Q R,S ⊂ E T , then the energy estimate in Proposition 6.1 reduces to the interior case in Proposition 2.1 for p = 2 and C 2 = 0. This remark also applies to the following Lemmas 6.1 -6.2.
with ̺ ∈ (0, 1), we introduce the numbers µ ± and ω satisfying The numbers µ ± and ω are used in the following lemmas.
Lemma 6.1. Let the hypotheses in Proposition 6.1 hold true and ξ ∈ (0, 1). There exist constants γ > 0 and c o ∈ (0, 1) depending only on the data and the structure of ∂E, such that if γC 2 ̺ < ξω and The proof runs almost verbatim as in Lemma 2.1. One needs only to intersect all cylinders and cubes with E T and take p = 2. As for the term of C 2 , it is absorbed into other terms on the right by imposing γC 2 ̺ < ξω. Certainly, if Q ̺ ⊂ E T , the term of C 2 does not appear and we are back to Lemma 4.1 for p = 2.
We present the following lemma on expansion of positivity. Lemma 6.2. Let the hypotheses in Proposition 6.1 hold true. Suppose that µ + − 1 4 ω > 0 and that for some α ∈ (0, 1) and some t * ∈ (t o − ̺ 2 , t o ) there holds Then there exist positive constants γ, δ, η and κ depending only on the data and the structure of ∂E, such that either γC 2 ̺ ≥ ω or for all times t * + 1 2 δ̺ 2 ≤ t ≤ t * + δ̺ 2 , Repeating in this fashion, we obtain the following construction for n ≥ 0: By induction, we have for all n ≥ 0 that ess osc Then the derivation of the modulus of continuity inherent in the above estimate can be performed as in the Section 3.5.
Finally, we recall that the vertex of Q r could be any point in E T and therefore the above oscillation decay holds in the interior and up to the lateral boundary. This joint with the continuity up to the initial level presented in Proposition 5.1 would allow us to conclude the proof of Theorem 1.2, just like in Section 5 for the Dirichlet problem.

UNIFORM APPROXIMATIONS
Existence of weak solutions to the Stefan problem (1.2), given proper initial and boundary conditions, is an issue of independent interest. A standard device in the construction of weak solutions to boundary value problems consists in first solving regularized versions, deriving a priori estimates that suggest where to find a solution, and then obtaining a solution in a proper limiting process via the a priori estimates and compactness arguments. In general, such a limiting process requires the function A(x, t, u, ξ) to satisfy more stringent structural conditions than (1.3), such as monotonicity in ξ, cf. [22, Chapter V, Theorem 6.7]. The so-obtained limit function will be in the function space L ∞ 0, T ; L 2 (E) ∩ L p 0, T ; W 1,p (E) , and meanwhile, verifies one of the integral formulations in Section 1.2, cf. [18,26] and [22, Chapter V, Section 9]. Nevertheless, Theorems 1.1 -1.2 do not grant continuity to this kind of solution as it does not possess time derivative in the Sobolev sense. A priori estimates on the time derivative are generally unavailable. On the other hand, the importance of continuous weak solutions (temperatures) lies in their physical bearings.
The purpose of the present section is to exhibit that the arguments presented in the previous sections permit us to identify a limit function, which is continuous with the same kind of moduli as in Theorems 1.1 -1.2, with the aid of the Ascoli-Arzela theorem (cf. [8,Chapter 5,Section 19]). This, joint with the existence results described above, will yield continuous weak solutions without any knowledge on the time derivative.
Consider the regularized Dirichlet problem: Here the boundary datum g and the initial datum u o are as in Theorem 1.1. For each fixed ε, the notion of solution to (7.1) can be defined via a similar integral identity as in Section 1.2.2. A main difference now is that the function space for solutions becomes  , g ∞,S T } and is equicontinuous in E T . More precisely, there exist positive constants γ and q depending only on the data and α * , and a modulus of continuity ω(·), determined by the data, α * ,̺, M , ω o (·) and ω g (·), independent of ε, such that for every pair of points (x 1 , t 1 ), (x 2 , t 2 ) ∈ E T . In particular, there exists σ ∈ (0, 1) depending only on the data and α * , such that if ω g (r) ≤ C g | ln r| λ and ω o (r) ≤ C uo | ln r| λ for all r ∈ (0,̺), where C g , C uo > 0 and λ > σ, then the modulus of continuity is for all r ∈ (0,̺) with some C > 0 depending on the data, λ, α * , M , C uo and C g .
Similarly we may consider a regularized Neumann problem: Here the boundary datum ψ and the initial datum u o are as in Theorem 1.2. The notion of solution to (7.3) can be defined via a similar integral identity as in Section 1.2.3. In particular, the function space in (7.2) is used. The following theorem can be regarded as a "cousin" of Theorem 1.2.
for every pair of points (x 1 , t 1 ), (x 2 , t 2 ) ∈ E T . In particular, there exists σ ∈ (0, 1) depending only on the data, C 2 and the structure of ∂E, such that if ω o (r) ≤ C uo | ln r| λ for all r ∈ (0, 1), where C uo > 0 and λ > σ, then the modulus of continuity is ω(r) = C ln 1 r − σ 2 for all r ∈ (0, 1), with some C > 0 depending on the data, the structure of ∂E, C 2 , M and C uo .
Remark 7.1. Theorems 7.1 -7.2 lay the foundation of obtaining continuous solutions to the Dirichlet problem (1.5) or the Neumann problem (1.6), once we can construct solutions {u ε } to (7.1) or (7.3) and identify the limit function u as a proper weak solution to (1.5) or (1.6). Clearly, the so-obtained solution u will enjoy the logarithmic type modulus of continuity.
The second term containing H ′ ε on the right-hand side is discarded because now ζ = 0 on ∂ P Q R,S . Using these remarks, one can perform the De Giorgi iteration in Lemma 2.1 and reach the same conclusion.
As for Lemma 2.2, now the condition becomes |µ − | ≤ ξω and, letting k = µ − + ξω, The derivation of modulus of continuity runs similar to Section 3.5; here it is affected only by an extra control of order ε or ε 1 1+q . Therefore we may conclude the proof of Theorem 7.1.