On the two-phase fractional Stefan problem

The classical Stefan problem is one of the most studied free boundary problems of evolution type. Recently, there has been interest in treating the corresponding free boundary problem with nonlocal diffusion. We start the paper by reviewing the main properties of the classical problem that are of interest for us. Then we introduce the fractional Stefan problem and develop the basic theory. After that we center our attention on selfsimilar solutions, their properties and consequences. We first discuss the results of the one-phase fractional Stefan problem which have recently been studied by the authors. Finally, we address the theory of the two-phase fractional Stefan problem which contains the main original contributions of this paper. Rigorous numerical studies support our results and claims.


Introduction
In this paper we will discuss the existence and properties of solutions for the well-known classical Stefan problem and the recently introduced fractional Stefan problem. A main feature of such problems is the existence of a moving free boundary, which has important physical meaning and centers many of the mathematical difficulties of such problems. For a general presentation of Free Boundary problems a classical reference is [16]. For the classical Stefan problem see Section 2 below. The Stefan problems considered here can be encoded in the following general formulation where the diffusion operator L is chosen as follows: If L = −∆, then (1.1) is called the classical/local Stefan problem.
There is a further choice consisting of considering both types of problems with one and two phases. More precisely, given a constant L > 0 (latent heat) and k 0 , k 1 , k 2 > 0 (thermal conductivities), we take  The first work on the fractional Stefan problem that we know of is due to Athanasopoulos and Caffarelli in [2] where it is proved that the temperature u is a continuous function in a general setting that includes both the classical and the fractional cases. This is followed up by [15], where detailed properties of the selfsimilar solutions, propagation results for the enthalpy and the temperature, rigorous numerical studies as well as other interesting phenomena are established. Other nonlocal Stefan type models with degeneracies like (1-Ph) and (2-Ph) have also been studied. We mention the recent works [4,7,6,8], where it is always assumed that L is a zero order integro-differential operator. There are also some models involving fractional derivatives in time, see e.g. [26].
Organization of the paper. In Section 2 introduce the classical Stefan Problem from a physical point of view, mainly to fix ideas and notations and also to serve as comparison with results on the fractional case. We will also give classical references to the topic and discuss how to deduce the global formulation (1.1).
We address the basic theory of fractional filtration equations in the form of (1.1) in Section 3. We discuss first the existence, uniqueness and properties of bounded very weak solutions. Later we address the basic properties of a class of bounded selfsimilar solutions. Finally, we present the theory of finite difference numerical schemes.
We devote Section 4 to the one-phase fractional Stefan problem, which has been studied in great detail in our paper [15]. We describe there the main results obtained in that article.
Section 5 contains the main original contribution of this article, which regards the two-phase fractional Stefan problem. First, we establish useful comparison properties between the onephase and the two-phase problems. Later, we move to the study of a selfsimilar solution of particular interest. Thus, in Theorem 5.4 and Theorem 5. 6 we construct a solution of the two-phase problem which has a stationary free boundary, a phenomenon that cannot occur in the one-phase problem. Finally, we move to the study of more general bounded selfsimilar solutions. Theorem 5.7 establishes the existence of strictly positive interface points bounding the water region {h ≥ L} and the ice region {h ≤ 0}. In particular, this shows the existence of a free boundary. Theorem 5.10 establishes the existence of a nonempty mushy region {0 < h < L} in the case s = 1/2. The last mentioned results are highly nontrivial and require original nonlocal techniques.
Rigorous numerical studies support also the existence of a mushy region in cases when s = 1/2. We comment on this fact in Section 6.
Section 7 is devoted to the study of some propagation properties of general solutions. We also support these results with numerical simulations that present interesting phenomena that were not present in the one-phase problem.
Finally, we close the paper with some comments and open problems.

The classical Stefan problem
The Classical Stefan Problem (CSP) is one of the most famous problems in present-day Applied Mathematics, and with no doubt the best-known free boundary problem of evolution type. The mathematical formulation is based on the standard idealization of heat transport in continuous media plus a careful analysis of heat transmission across the change of phase region, the typical example being the melting of ice in water. More generally, by a phase we mean a differentiated state of the substance under consideration, characterized by separate values of the relevant parameters. Actually, any number of phases can be present, but there are no main ideas to compensate for the extra complication, so we will always think about two phases, or even one for simplicity (plus the vacuum state).
The CSP is of interest for mathematicians, because it is a simple free boundary problem, easy to solve today when N = 1, but still quite basic problems are open for N > 1. It is always interesting for physicists, since there exist several processes of change of phase which can be reduced to the CSP. Finally, it is of interest for Engineers, since many applied problems can be formulated as CSPs, like the problem of continuous casting of steel, crystal growth, and others.
Though understanding change of phase has been and is still a basic concern, the mathematical problem combines PDEs in the phases plus a complicated geometrical movement of the interphase, and as a consequence, the rigorous theory took a long time to develop. A classical origin of the mathematical story are the papers by J. Stefan who around 1890 proposed the mathematical formulation of the later on called Stefan Problem also in dimension N = 1 when modelling a freezing ground problem in polar regions, [23]. He was motivated by a previous work of Lamé and Clapeyron in 1831 [20] in a problem about solidification. The existence and uniqueness of a solution was published by Kamin as late as 1961 [19] using the concept of weak solution. Progress was then quick and the theory is now very well documented in papers, surveys, conference proceedings, and in a number of books like [22], [21], and the very recent monograph by S. C. Gupta [18].
2.1. The Classical Formulation. A further assumption which we take for granted in the classical setting is that the transition region between the two phases reduces to an (infinitely thin) surface. It is called the Free Boundary and it is also to be determined as part of the study.
• With this in mind let us write the basic equations. First, it is useful to have some notation. We assume that both phases occupy together a fixed spatial domain D ⊂ R N , and consider the problem in a time interval 0 ≤ t ≤ T , for some finite or infinite T . On the other hand, the regions occupied by each of the phases evolve with time, so the liquid (water in the standard application) will occupy Ω 1 (t) and the solid (ice) Ω 2 (t) at time t. Clearly, for all t Ω 1 (t) ∪ Ω 2 (t) = D. The initial location of the two phases, Ω 1 (0) and Ω 2 (0), is also known. Let us introduce some domains in space-time: Q T := D × (0, T ), and let The energy balance in the liquid takes the form of the usual linear heat equation while for the solid we have (E2) c 2 ρ∂ t u 2 = k 2 ∆u 2 in Ω 2 .
Here u 1 and u 2 are the resp. temperatures in the liquid and solid regions, while k 1 and k 2 are the resp. thermal conductivities, c 1 and c 2 the specific heats, and ρ is the density. All of these parameters are usually supposed to be constant (just for the sake of mathematical simplicity). It is however quite natural to assume that they depend on the temperature but then we have to write ∂ t (ρc i u i ) and ∇ · (k i ∇u i ). In this paper, we will always consider ρ = c i = 1.
• Next we have to describe what happens at the surface separating Ω 1 and Ω 2 , i.e., the free boundary, Γ. A first condition is the equality of temperatures, This is also an idealization, other conditions have been proposed to describe more accurately the transition dynamics and are currently considered in the mathematical research.
We need a further condition to locate the free boundary separating the phases. In CSP this extra condition on the free boundary is a kinematic condition, describing the movement of the free boundary based on the energy balance taking place on it, in which we have to average the microscopic processes of change of phase. The relevant physical concepts are heat flux and latent heat. The result is as follows: if Ψ = k 2 ∇u 2 − k 1 ∇u 1 is the heat flux across Γ, then (FB2) Ψ is parallel to the space normal n to Γ and Φ = L v , v being the velocity with which the free boundary moves. The constant L > 0 is called the latent heat of the phase transition. In the ice/water model it accounts for the work needed to break down the crystalline structure of the ice. Relation (FB2), called the Stefan condition, is not immediate. It is derived from the global physical formulation in the literature. Equivalently, if G(x, t) = 0 is the implicit equation for the free boundary Γ in (x, t)-variables, (FB2) can be written as Ψ · ∇ x G + ρL∂ t Φ = 0.
All things considered, we have the complete problem as follows: Problem about classical solutions. Given a smooth domain D ⊂ R N and a T > 0, we have to: (i) Find a smooth surface Γ ⊂ Q T = D × (0, T ) separating two domains in space-time Ω 1 , Ω 2 . (ii) Find a function u 1 that solves (E1) in Ω 1 and a function u 2 that solves (E2) in Ω 2 in a classical sense, Typically we require u i ∈ C 2,1 x,t inside its domain Ω i , i = 1, 2. (iii) On Γ the free boundary conditions (FB1) and (FB2) hold. (iv) In order to obtain a well-posed problem we add in the standard way initial conditions (v) Boundary conditions on the exterior boundary of the whole domain ∂D for the time interval under consideration. These conditions may be Dirichlet, Neumann, or other type.
The precise details and results can be found in the mentioned literature. Let us remark at this point that it is the Stefan condition (FB2) with L = 0 that mainly characterizes the Stefan problem, and not the possibly different values of c, k on both phases.

2.2.
The one-phase problem. Special attention is paid to the simpler case where one of the phases, say the second, is kept at the critical temperature (e.g., in the water-ice example, the ice is at 0 • C). Then the classical problem simplifies to: Finding a subset Ω of Q T bounded by a internal surface Γ = ∂Ω∩Q T and a function u(x, t) ≥ 0 such that where v is the normal speed on the advancing free boundary. The theory for the one-phase problem is much more developed, and essentially simpler.
2.3. The global formulation. In order to get a global formulation we re-derive the model from the general energy balance plus constitutive relations. In an arbitrary volume Ω ⊂ D of material we have d dt where Q(Ω) =´Ω e(x, t)ρ(x, t) dx is the energy contained in Ω at time t, Ψ(∂Ω) =´∂ Ω ϕ(x, t) · n dS is the outcoming energy flux through the boundary ∂Ω, and F(Ω) =´Ω f (x, t) dx is the energy created (or spent) inside Ω per unit of time. Therefore, e represents an energy density per unit of mass (actually an enthalpy). We need to further describe these quantities by means of constitutive relations. One of them is Fourier's law, according to which where u(x, t) is the temperature and k > 0 is the heat conductivity, in principle a positive constant. Thus, we get the global balance law It is useful at this stage to include ρ into the function e by defining a new enthalpy per unit volume, h = ρe. Equivalently, we may assume that ρ = 1. Using Gauss' formula for the first integral in the second member we arrive at the equation This is the differential form of the global energy balance, usually called enthalpy-temperature formulation. We have now two options: (i) Either assuming the usual structural hypothesis on the relations between h, u, and k, and performing a partial analysis in each phase, deriving the equations (E1), (E2), plus a free boundary analysis leading to the free boundary conditions (FB1), (FB2), or (ii) trying to continue at the global level, avoiding the splitting into cases. If we take the latter option which allows for a greater generality and conceptual simplicity.
We will take this latter option, which allows us to keep a greater generality and conceptual simplicity. We only need to add a structural relation linking h and u. This is given by the two statements: (ii) At u = 0 we have a discontinuity. More precisely, h jumps from 0 to L > 0 at u = 0.
After some easy manipulations contained in the literature we get the relations stated in the Introduction and the integral formulation in Definition 3.2 with −(−∆) s replaced by ∆. If the space domain D is bounded, we need boundary conditions on the fixed external boundary of D. We see immediately that this is an implicit formulation where the free boundary does not appear in the definition of solution.

Common theory for nonlinear fractional problems
The theory of well-posedness and basic properties for fractional Stefan problems can be seen as a part of a more general class of problems that we call Generalized Fractional Filtration Equations (see [24] for the local counterpart). More precisely, one can consider the equation for s ∈ (0, 1), N ≥ 1, and Together with (3.1) one needs to prescribe an initial condition h(·, 0) = h 0 .
Our theory is developed in the context of bounded very weak (or distributional) solutions. More precisely: 3.1. Well-posedness and basic properties. The following result ensures existence and uniqueness (see [17]).
We will also need some extra properties of the solution. For that purpose, we rely on the argument present in Appendix A in [15] where bounded very weak solutions are obtained as a monotone limit of L 1 ∩ L ∞ very weak solutions. The general theory for the latter comes from [14,13]. See also [9,10] for the theory in the context of weak energy solutions.
Additionally, the results of [2] ensure that for fractional Stefan problems, the temperature is a continuous function. We refer to Appendix A in [15] for an explanation of how the result [2] is applied to our concept of solutions.

3.2.
Bounded selfsimilar solutions. The family of equations encoded in (3.1) admits a class of selfsimilar solutions of the form for any initial data satisfying h 0 (ax) = h 0 (x) for all a > 0 and all x ∈ R N . It is standard to check the following result, and we refer the reader to [15] for details.
, where the selfsimilar profile H satisfies the stationary equation When N = 1, we can choose a more specific initial data that will lead to a more specific selfsimilar solution from which we will be able to prove several properties for the general solution of (3.1). Indeed, we have the following Theorem which is new in the general context we are treating.
Theorem 3.8. Under the assumptions of Theorem 3.7, and additionally, that N = 1, and that for some Then the corresponding solution h ∈ L ∞ (Q T ) is selfsimilar as in Theorem 3.7. Moreover, it has the following properties: Proof. Part (a) follows by translation invariance and uniqueness of the equation The bounds in (b) are a consequence of comparison and the fact that any constant is an stationary solution of (3.1). The limits in (b) are obtained by selfsimilarity and the fact that the initial condition is taken in the sense of Remark 3.3 (see Lemma 3.13 in [15] for more details). Finally, (c) follows from Theorem 3.6 and H(ξ) = h(ξ, 1).
Remark 3.9. By translation invariance, one can obtain selfsimilar solutions not centred at x = 0 by just considering h 0,c = h 0 (· + c) for any c ∈ R. In this way, one obtains selfsimilar profiles of the form H c = H(·+c). Moreover, selfsimilar solutions in R also provide a family of selfsimilar solutions in R N by extending the initial data constantly in the remaining directions. See Section 3.1 in [15] for details.
3.3. Numerical schemes. As in [15], we can have a theory of convergent explicit finitedifference schemes (see also [14]). More precisely, we discretize (3.1) by where V is the approximation of the enthalpy defined in the uniform in space and time grid On the other hand, L ∆x is a monotone finite-difference discretization of (−∆) s (see e.g. [13]). It takes the form: where ω γ,∆x = ω −γ,∆x are nonnegative weights chosen such that the following consistency assumption hold: . Together with (3.3) one needs to prescribe an initial condition. Since h 0 is merely L ∞ we need to take [15] (see also [14]), we get the following convergence result.
Theorem 3.10. Assume (A Φ ). Let h ∈ L ∞ (Q T ) be the very weak solution of (3.1) with h 0 ∈ L ∞ (R N ) as initial data such that h 0 − h 0 (· + ξ) ∈ L 1 (R N ) for all ξ > 0, ∆t, ∆x > 0 be such that ∆t ∆x 2s , L ∆x be such that (3.4) and (3.5) hold, and V j β be the solution of (3.3). Then, for all compact sets K ⊂ R N , we have that The above convergence is the discrete version of convergence in C([0, T ]; L 1 loc (R N )). Remark 3.11. We would like to mention that all the results of Section 3 apply also in the local case, i.e., replacing (−∆) s by −∆ in (3.1). More precisely, • The existence part as in Theorem 3.4 is a classical matter (see [24]). We also refer to Appendix A in [15] for a modern reference in a more general local-nonlocal context. • Properties as in Theorem 3.5 follow from the results in Appendix A in [15]. See also [12,11].
• Regularity of Φ(h) as in Theorem 3.6 is the classical result of Caffarelli and Evans in [5].
• Convergence of numerical schemes as in Theorem 3.10 follows from the results of [14] replacing L ∆x in (3.4) by the standard monotone finite-difference discretization of the Laplacian:

The one-phase fractional Stefan problem
Here we list a series of important results regarding the one-phase fractional Stefan problem where Φ 1 is given by (1-Ph) with k 0 = 1. Since such Φ 1 is locally Lipschitz, all the results listed in Section 3 apply for the very weak solution h of (4.1). Moreover, we list (without proofs) a series of interesting results recently obtained in [15].
We start by stating the fine properties of the selfsimilar profile.
Theorem 4.1. Assume that Φ 1 is given by (1-Ph) with L > 0, and let the assumptions of Theorem 3.8 hold with b 1 = L + P 1 and b 2 = L − P 2 for P 1 , P 2 > 0. The profile H has the following additional properties: (a) (Free boundary) There exists a unique finite ξ 0 > 0 such that H(ξ 0 ) = L. This means that the free boundary of the space-time solution h(x, t) at the level L is given by the curve for all t ∈ (0, T ).
Moreover, ξ 0 > 0 depends only on s and the ratio P 2 /P 1 (but not on L).
for some α > 0, and (SSS) is satisfied in the classical sense in R \ {ξ 0 }.
Again, we remind the reader that selfsimilar solutions in R also provide a family of selfsimilar solutions in R N by extending the initial data constantly in the remaining directions. Once the above properties are established in that case as well, one can prove that the temperature u := (h − L) + has the property of finite speed of propagation under very mild assumptions on the initial data.
Theorem 4.2 (Finite speed of propagation for the temperature). Let h ∈ L ∞ (Q T ) be the very weak solution of (4.1) with h 0 ∈ L ∞ (R N ) as initial data and u : Moreover, the temperature not only propagates with finite speed, but it also preserves the positivity sets, an important qualitative aspect of the solution.  (a) If h 0 ≥ L + ε > L in B ρ (x 1 ) for x 1 ∈ R N and ρ, ε > 0, then h(·, t) > 0 for all t ∈ (0, T ). (b) If additionally supp{h 0 } ⊂ B η (x 0 ) for x 0 ∈ R N and η > 0 , then h(x, t) 1/|x| N +2s for all t ∈ (0, T ) and |x| large enough.
The question of asymptotic behaviour is still under study, but we refer to the preliminary results of the one-phase work [15].

The two-phase fractional Stefan problem
In this section we treat the two-phase fractional Stefan problem, i.e., where Φ 2 is given by the graph (2-Ph). Again, we make the choice k 1 , k 2 = 1.

5.1.
Relations between one-phase and two-phase Stefan problems. Here we will see that any solution of the two-phase Stefan problem is essentially bounded from above and from below by solutions of the one-phase Stefan problem.
We need two lemmas to prove this result.
That is, . Finally, the initial data relation follows from Remark 3.3.

5.2.
A selfsimilar solution with antisymmetric temperature data. We continue to address properties of the same type as in Theorem 4.1 for the two-phase Stefan problem. In the case where the initial temperature is an antisymmetric function, the solution has a unique interphase point between the water and the ice region, and it lies at x = 0 for all times (stationary interphase). As a consequence, the enthalpy h is continuous for all x = 0 and discontinuous at x = 0 for all times t > 0. This is the precise result: Theorem 5.4. Assume N = 1, P > 0, and Let h be selfsimilar solution (given by Theorem 3.8) of (5.1) with initial data h 0 . Let also H and U be the corresponding profiles. Then, additionally to the properties given in Theorem 3.8, we have that: To prove this theorem, we need a simple lemma. Proof. ClearlyΦ 2 (h) = Φ 2 (h) in Q T and ∂ t h = ∂ th in D (Q T ). The initial data relation follows from Remark 3.3.
Proof of Theorem 5.4. 1) Antisymmetry. We consider the translated problem as in Lemma 5.5 and prove that u(·, t) and h(·, t) are antisymmetric. We recall that the initial datum is and Φ 2 (h) := max{h − L 2 , 0} + min{h + L 2 , 0}. We avoid the superscript tilde onΦ 2 andh in the rest of the proof for convenience. .
. Note also that h 0 (x) = −h 0 (−x) and thus h 1 is a very weak solution with initial data h 0 . By uniqueness, this implies that h 1 = h, which proves the antisymmetry result. The antisymmetry of u(·, t) follows. Note that the translation did not affect the u.
One might be tempted to try to use the behaviour at the interphase of the above constructed solution to obtain estimates close to the free boundary of a general solution in the spirit of Theorem 4.1(d). However, in this special case, the behaviour at the interphase is quite different, as we show in the following result.
where P s is the fractional heat kernel. In particular, for ξ is close enough to 0 and for some C = C(P, s) > 0.
Proof. According to the previous results, the interphase points coincide and we know that H(ξ) > L for ξ < 0 and H(ξ) < 0 for ξ > 0. Consequently, for all t ∈ (0, T ), we have that We examine the first term of the very weak formulation for h: By using the above relation in the definition of very weak solution for h, we get with u 0 := Φ 2 (h 0 ). To sum up, u = Φ 2 (h) is the unique very weak solution of the fractional heat equation with initial data Consequently, u is given by the convolution formula (5.2) (see e.g. [3]). Moreover, it is wellknown that P s (z, t) is a smooth function for all t > 0 and z ∈ R, and that it has the following properties P s (z, t) = P s (−z, t) and P s (z, t) t (t 1/s + |z| 2 ) (1+2s)/2 . Then, for ξ < 0, we have Since P s is a positive and smooth kernel, where C = 2P P s (0, 1). The bound for ξ > 0 follows by antisymmetry.

5.3.
Analysis of general selfsimilar solutions. Now, we analyze the fine properties of general selfsimilar solutions where the initial temperature is not antisymmetric, i.e. P 1 = P 2 . The main difference will be that the interface is never stationary. We may assume that P 1 > P 2 without loss of generality; the case P 2 > P 1 is obtained by antisymmetry. The solution constructed in Theorem 5.4 will be used in the analysis.
Our running assumptions during this section will be N = 1, P 1 > P 2 > 0, and Denote h as the selfsimilar solution of (5.1) with initial data h 0 . Let also H and U be the corresponding profiles.
Our first main result in the section (see Theorem 5.7 below) establishes the existence of strictly positive interface points bounding the water region and the ice region, as well as the behaviour in the mushy region if it exists.
Our second main result (see Theorem 5.10 below), restricted to the case s = 1/2, establishes the existence of a nonempty mushy region lying in the positive half space.
Theorem 5.7. Under the running assumptions, additionally to the basic properties given in Theorem 3.8, we have that: (a) (Unique interphase points) There exist unique points ξ w and ξ i with 0 < ξ w ≤ ξ i < +∞ such that H(ξ − w ) = L, H(ξ) > L for ξ < ξ w , H(ξ + i ) = 0, H(ξ) < 0 for ξ > ξ i . This means that the free boundaries of the space-time solution h(x, t) at the levels L and 0 are given by for all t ∈ (0, T ).
(b) (Improved monotonicity in the mushy region) If ξ w = ξ i , then H is strictly decreasing and smooth in [ξ w , ξ i ].
The proof of Theorem 5.7 will be divided into two parts. In the first part we will prove everything except the fact that ξ w > 0, obtaining only ξ w ≥ 0. The analysis for the strict inequality requires more refined and elaborate arguments.
Proof of Theorem 5.7 (Part 1). We do not prove the results in the order stated.
1) Interphase points. Recall that U is continuous, nonincreasing, and has limits P 1 and −P 2 at −∞ and +∞. Then, if we define the mushy region as the set this is either a closed finite interval M = [ξ w , ξ i ] or just a point when ξ w = ξ i . If M is just one point, then of course there is a unique interphase point for the ice and water region with no mushy region. Note that in the already studied case P 1 = P 2 , we are back to the setting of Theorem 5.4 where ξ w = ξ i = 0 and there is no mushy region.
with U (ξ) ≥ 0 for ξ ≤ −1 and U (ξ) ≤ 0 for ξ ≥ 1. Since U is continuous and takes the limits P 1 and −P 2 at −∞ and +∞, there exist a, b ≥ 1 such that which is a contradiction. The argument for smoothness inside this region is the same as in the one-phase problem.
3) Unique interphase points. By strict monotonicity in [ξ w , ξ i ], we have that for small enough ρ > 0, H(ξ w + ρ) < H(ξ w ) ≤ H(ξ − w ) = L, which proves that the only interphase point with the water region can be ξ w . A similar argument shows that ξ i is the only interphase point with the ice region.

4)
The interphase points are nonnegative: ξ w ≥ 0. Consider the solutionh of (5.1) with initial datah 0 given as h 0 with P 1 = P 2 (cf. Theorem 5.4). By comparison, we get that h ≥h or H(ξ) = h(ξ, 1) ≥h(ξ, 1) for a.e. ξ ∈ R. Sinceh(0 − , 1) = L, the continuity of (H(ξ) − L) + gives H(0 − ) ≥ L. Now, if H(0 − ) > L, then monotonicity and continuity gives ξ w > 0, and if H(0 − ) = L, then Step 3) gives that ξ w = 0. Note also that Theorem 5.6 gives that H(ξ) C|ξ| near the origin for ξ < 0, and this will be used below. Now, we improve the information about the interphase points ξ w , ξ i . The only thing left to show in Theorem 5.7 is the following result. Proof. Note that we know that 0 ≤ ξ w ≤ ξ i by Part 1 of the proof of Theorem 5.7. Assume by contradiction that ξ w = ξ i = 0. Since the interphase points are unique, we proceed as in the proof of Theorem 5.6 to show that u = Φ 2 (h) is the unique very weak solution of the fractional heat equation with initial data Finally, we note, by symmetry of the heat kernel P s in the first variable and the fact that P 1 > P 2 we have that which is a contradiction with our assumption.
Then, for ξ ∈ (0, ξ i ) close enough to 0 we have that Finally, taking limits as ξ 1 → 0 + in (5.3) we get which is a contradiction and shows that ξ w > 0.
The only thing left to show is assertion (5.4) holds if ξ w = 0. We will do it in a series of steps.
1) Let us consider g : [0, +∞) → R as an auxiliary function defined as g := U , the known continuous and bounded solution. By equation (SSS) we get that U satisfies the following external boundary value problem for ξ ∈ [0, +∞).

5)
The strict inequality C 1 > C 2 needs an extra argument done by approximation. For convenience, now we will use the notations: U P 1 := U , ξ i,P 1 := ξ i and g P 1 := g since the dependence on the parameter P 1 will be important. We fix P 1 and P 2 and pick a value P * 1 := P 1 − ε > P 2 for some small ε > 0. We then consider U P * 1 , ξ i,P * 1 , and g P * 1 . By comparison, we have U P 1 ≥ U P * 1 , and ξ i,P 1 ≥ ξ i,P * 1 . In particular, we have that g P 1 * (ξ) ≤ g P 1 (ξ) for ξ > 0. On one hand, this estimate implies that On the other hand, from Step 3) we get that small ξ < 0 we havẽ U P 1 * (ξ) ≥ C * 1 |ξ| s , and Step 4) for U P * 1 implies that C * 1 ≥ C * 2 . The final point is to note that due to the initial data of the space-time version,Ũ P 1 (ξ) andŨ P * 1 (ξ) are proportionalŨ P 1 (ξ) = λŨ P * 1 (ξ) with λ = (P 1 /(P 1 − ε)) > 1, so that C 1 > C * 1 . We conclude that and thus, for some c > 0. This ends the proof.
As a further result, we prove the existence of a mushy region for s = 1/2.
Proof. We only need to show that the configuration 0 < ξ w = ξ i is not possible when s = 1/2. Assume that 0 < ξ w = ξ i . We follow first the ideas of the proof Theorem 5.6 and show that u satisfies a fractional heat equation, this time with a forcing term. More precisely (we take T = 1 for simplicity), for ψ ∈ C ∞ c (R × [0, 1)), we have that, Note that, for the last term above, we have that where δ ξwt 1 2s (x) denotes the Dirac delta measure in the variable x at the point ξ w t 1 2s . By using the above relations in the definition of very weak solution for h, we have that u satisfies the following identity, 2s (x) dt. This means that u is the distributional solution of the following heat equation with a forcing term: Note that u =ũ +û where On one hand, note that On the other hand, we can use Duhamel's representation formula to show that Using the above estimate in x = ξ w t 1 2s for some t > 1/2 we have that Of course this is a contradiction with the fact that u is bounded.
6. Numerical evidence on the existence of a mushy region Since we have established a good numerical framework for the fractional Stefan problem, the form of the selfsimilar solutions of both the one-phase and the two-phase problems have been verified numerically. In particular, in examining the two-phase problem we have proved partial results of the existence of a mushy region, it is theoretically established only in the case s = 1/2 and the numerics agrees. Moreover, the existence of such a mushy zone is also observed for other values of s ∈ (0, 1) in view of the numerical simulations presented in Figure  2, which are very clear in the case when P 2 is very close to 0.

Figure 2.
Existence of mushy regions for s > 1/2 and P 2 small. This is consistent with the results ensuring existence of a mushy region obtained in the one-phase problem, and the L 1 loc continuous dependence result on the initial data as P 2 → 0. However, the numerical simulations are not conclusive for P 2 close to P 1 and s > 1/2, as we can see in Figure 3. In the figures the enthalpy is displayed and L = 1.

Results on the speed of propagation
In the one-phase Stefan problem, knowing the precise behaviour of the selfsimilar solution is enough to conclude that the temperature has finite speed of propagation (see Theorem 4.2). However, in the two-phase problem an analogous result is simply not true.
In this section we present a series of partial results regarding the speed of propagation of the temperature, as well as numerical simulations. They exhibit interesting and very different behaviour. Figure 3. Solutions for s > 1/2 and P 2 close to P 1 . Proposition 7.1 (Control through one-phase selfsimilar solutions). Let h ∈ L ∞ (Q T ) be the very weak solution of (5.1) with h 0 ∈ L ∞ (R N ) as initial data and u := Φ 2 (h). If supp{Φ 2 (h 0 + ε) + } ⊂ B R (x 0 ) for some ε > 0, R > 0, and x 0 ∈ R N , then supp u(·, t) + ⊂ B R+ξ 0 t 1 2s (x 0 ) for some ξ 0 > 0 and all t ∈ (0, T ). However, we cannot at the same time conclude that u − is compactly supported. In fact, this is not true as one can see from Figure 4 where it has finite or infinite speed of propagation depending on the initial amount of entalphy above the level L and below the level 0. Proposition 7.2 (Control through two-phase selfsimilar solutions). Let h ∈ L ∞ (Q T ) be the very weak solution of (5.1) with h 0 ∈ L ∞ (R N ) as initial data and u := Φ 2 (h). If h 0 ≤ −C outside B R (x 0 ) for some C > 0, R > 0, and x 0 ∈ R N , then: (a) supp u(·, t) + ⊂ B   5 exhibits the control obtained by the two-phase selfsimilar solution described in e.g. Theorem 5.7. We see that the water region is first expanding, then contracting, and finally disappearing. Only the expansion was possible in the one-phase Stefan problem (cf. Theorem 4.3). In the end, all we can say is that the water region is not expanding more than what the two-phase selfsimilar solution allows it to do, and we obtain an upper estimate on the support. Similar results are also shown for the mushy region, but with a possibly bigger radius.

Comments and open problems
• The classical Stefan Problem has been thoroughly studied in the last decades with enormous progress, but still many fine details are under development. There are also plenty of studies of variants of the equation or systems which appear in physical or technical applications.
• We have concentrated our interest in presenting our results on the fractional Stefan problem, which are quite recent. A basic theory is ready in the one-phase fractional Stefan problem. The regularity of solutions and free boundaries is still at an elementary stage and needs much further development, even in the one-phase problem.
• The two-phase fractional Stefan problem with k 1 = k 2 has not been considered here and needs attention because different conductivity in the two phases agree with the practical evidence.