The Hele-Shaw problem with surface tension in the case of subdiffusion

In this paper we analyze anomalous diffusion version of the 
multidimensional Hele-Shaw problem with a nonzero surface tension 
of a free boundary. We prove the one-valued solvability of this 
moving boundary problem in the Holder classes. In the 
two-dimensional case some numerical solutions are constructed.

Here p and Υ(t) are a desired function and an unknown interface, correspondingly; ∆ y = ∇ 2 y , ∇ y = ∂ ∂y1 , ..., ∂ ∂yn ; µ, γ are positive constants; n t and n are unit normals to Υ(t) and Γ 0 , correspondingly, directed outside Ω(t); the quantity κ(Υ(t)) defined on Υ(t) denotes the mean curvature of this surface. We use the sign convention that convex hypersurfaces have positive mean curvature. V ν nt is the fractional velocity of the boundary Υ(t) in the direction of the normal n t and is represented by (see, e.g. [39]) where ·, · is the notation of the inner product. D ν t Υ(t) denotes the Caputo fractional derivative with respect to t and is defined by (see (2.4.6.) in [21]) In this paper we focuses on two main questions. The first one is the one-to-one solvability of problem (1)-(4) locally in time in the Hölder classes if γ > 0. To this end we adapt the classical approach which is used for a free boundary problem in the case of normal diffusion (see, e.g. [2]) to the subdiffusion case. We reduce a moving boundary problem (1)-(4) to a nonlinear problem in a fixed domain, then linearize the problem and obtain one-valued solvability of the linearized problem. After that we use this result for the reduction of the nonlinear problem to a fixed point theorem.
In this route we have to solve a lot of technically difficulties which are related to the nonlocal behavior of the free boundary velocity (i.e. with fractional derivative in time). The main one of them deals with the investigation of a nonclassical boundary value problem with the Caputo derivative in a boundary condition Note that this is the principal model problem such that nonlinear problem (1)- (4) will inherit the main features of (7). Using Fourier and Laplace transformations, we obtain the solution of (7) as the convolutions: U = K D 2 ℘, ℘ = G f (see (74) and (75)). The kernels G and K can be represented only as integrals which contain the Wright functions. Thus nonlocal form of kernels is the distinguishing feature of the fractional case. As usual in the potential theory, to evaluate the functions U and ℘ it is necessary to describe well the properties of K and G. Unfortunately, in virtue of nonlocal representations for these kernels, it is impossible to get the suitable local estimates as well as in the case of integer order derivative (ν = 1) for G and K and their derivatives. We can obtain just the integral estimates of these kernels which are represented in Lemma 4.2. In this route we essential use the main properties of the Wright functions represented in Proposition 2. Generally speaking, Lemma 4.2, Propositions 2 and 4 are significant under investigation of problem (7) and their proofs contain the main difficulties caused by the presence of the fractional derivative in time.
The second issue considered in the paper is the construction of a numerical solution to problem (1)-(4) in the two-dimensional case. Here we implement like [41] an explicit boundary tracking method, which is related with nonlocal behavior of the free boundary velocity. We recall that in the case of the classical Stefan problem there exists a lot of approaches which allow to avoid explicit boundary tracking. For instance, the level set method [9] and enthalpy formulation [38] can reduce the problem with the moving boundary to the problem in the fixed domain. As it was shown in [40], the enthalpy formulation for a problem with the fractional derivative in general describes a various physical process. Level set method, in turn, demands to find the velocity V nt from the first condition in (2) which is possible only for ν = 1. Thus, these methods can not be applied for the free boundary problems governed by subdiffusion.
Note that this work is a continuation of the investigation of the one-or twophase fractional Hele-Shaw problems reported in [41]- [43]. Thus we will drop here some technical calculations that are insignificantly different from similar one in the papers [41]- [43].
The paper is organized as follows. In Section 2 we describe some useful auxiliary properties related to the fractional derivatives in time: Propositions 1-3. These results are essentially used throughout the paper. In Section 3, we reduce problem (1)-(4) with an unknown boundary Υ(t) to a nonlinear problem in a fixed domain Ω T = Ω × (0, T ) and formulate the main result, Theorem 3.1. In Subsection 3.2, we rewrite this nonlinear problem as Az = F(z), where z = (u, σ) and A is a linear operator, F is a nonlinear one. Section 4 is dedicated to the investigation of model problem like (7). Then in Section 5, we prove the main result of this paper, Theorem 3.1. In Subsection 5.1, adapting the technique of the regularizer for parabolic systems [24] and Schauder technique into the case of a fractional derivative, we prove the one-to-one solvability to the linear problem: Az = F (x, t) for small period of time, Theorem 5.1. Finally, in Subsection 5.2, based on results from Theorem 5.1, Proposition 1 and a fixed point theorem, we show that the operator A −1 F is a contraction one. Section 6 contains numerical study of free boundary problem (1)- (4) if Ω(t) ⊂ R 2 , t ∈ [0, T ]. Proofs of some auxiliary assertions which are applied in Section 4 are given in Appendix.
2. Function spaces and preliminaries. In order to investigate the free boundary problem (1)-(4) we need some definitions and auxiliary results.
Let D be a given domain in R n , D T = D × (0, T );x, x be any points inD, x =x; t, τ ∈ [0, T ], t = τ ; α, β ∈ (0, 1). In this paper we will use the following function spaces: . The spaces C([0, T ], C l+α (D)) are used by many authors and their definitions and properties can be found, for instance, in [31].
Denote by Next we define function spaces C l+α,β,α (D T ).
Definition 2.1. The Banach space C l+α,β,α (D T ) is the set of functions w(x, t) with the finite norm: x,D T and w (β) t,D T are Hölder constants of a function w(x, t) in x and t respectively.
Let d is a smooth surface and Definition 2.3. We will say w ∈ P 4+α .
Moreover, we will use the usual Hölder spaces C l+α (D) and C l+α (∂D), their definitions can be found, for instance, in [25].
Throughout the paper we will need in interpolation inequality (see Corollary 1.2.18 [31]): where ∂D ∈ C l2 , a ∈ (0, 1), 0 ≤ l 1 < l 2 . Next, we define the fractional Riemann-Liouville integral and derivative of a function w(·, t) with respect to t as (see, e.g., (2.1.1) and (2.1.8) [21]): The following results which are well known in the case of an integer order derivative will be essentially applied in Sections 4 and 7.1. ii: The results of Proposition 1 have been obtained in [43] (see Propositions 2.1 and 2.2 there).
Note that the Wright functions (see, e.g. [13] or [32]): play fundamental roles in various applications of the fractional calculus. Here, we formulate the main properties of this function: asymptotic representations, estimates, formula for fractional differentiation and integration, which will be essentially used to investigate the model problem in Section 4 and to obtain the corresponding coercive estimates of its solution. The proofs of the following properties for the Wright functions are represented in [22,34,35] and [23].
Proposition 2. Let C 1 , C 2 , A, b, c be positive constants, d ∈ R 1 , z ∈ C, then there are the following: i: ii: vi: where where either l > 0 and b > 0, or l > −1 and b = 0.
Remark 1. The statements of Proposition 3 are true if Ω ≡ R n + , and Γ 0 = ∅, g 1 ≡ 0, g 0 and g 2 are finite functions which satisfy conditions of Proposition 3.
3. The nonlinear mapping. Main results.
3.1. Reduction of problem (1)-(4) to a problem in the fixed domain. We assume that To prove the solvability of problem (1)-(4), it is convenient to reduce the one to a problem in a fixed domain. To this end, we use the Hanzawa method [17]. Let ω = (ω 1 , ..., ω n−1 ) be some coordinates on Υ. We represent Υ in the form y =m(ω) and denote byn(ω) the normal to Υ directed outside Ω.
We represent the free boundary in problem (1)-(4) as where ρ(ω, t) is an unknown function, and That means the surface Υ(t) in the local variables is given by Using (5), (31) and (33), we can rewrite boundary conditions in (2) as Let such that the transform e −1 ρ maps Ω(t) onto Ω × (0, T ) = Ω T , and Υ(t) onto Υ × [0, T ] = Υ T ; the free boundary is given by and ω(x), λ(x) are the coordinates in X T similar to the coordinates ω(y), η(y) in Y T . After the change of variables (36), we have a new desired function v(x 1 , ..., x n , t) = p(y 1 , ..., y n , t) • e ρ (x, t).
Denote by ∇ ρ = (E * ρ ) −1 ∇ x where E ρ is the Jacobi matrix of the mapping y = e ρ (x, t) so that ∇ y = ∇ ρ and ∆ y = ∇ 2 ρ . Taking into account that y = x near Γ 0T , we can reduce free boundary problem (1)-(4) to the nonlinear problem in the fixed domain: A 0 (ω, ρ, ∇ ω ρ) are some specific smooth functions (their representations and properties can be found in [2,41] and [10]) such that Moreover, one can easily check that Let us define the function v 0 (x) as a solution of the boundary-value problem: We assume that condition (30) hold, and According to the results of Chapter 3 in [25], there exists a unique solution v 0 (x) of problem (46), and where C is a positive constant.
Henceforward the letter C will be used to denote different constants encountered in our formulae. The main result of our paper is the following.
Thus system (39)- (43) can be written briefly in the form where a linear operator A is defined by the left-hand side of (53), and F is a nonlinear operator, (44), (48) and Corollary 1, we can deduce the following results.
Corollary 2. The functions F i (u, σ), i = 0, 2, F 1 (σ) contain the higher derivatives of u(x, t) and σ(ω, t) with the coefficients that tend to zero as t → 0, the "quadratic" terms of minor differential orders of unknown functions. Moreover, Note that equality (55) together with conditions (53) provide In a next step of our investigation we prove the boundedness of the linear operator A in the corresponding function spaces. To this end, we freeze the functional arguments in the functions F i (u, σ),i = 0, 2, F 1 (σ). Then system (54) will be a linear system with variable coefficients, which will be studied in detail in Subsection 5.1. This investigation will be based on the properties of solutions to the corresponding model problems.
4. Model problem with fractional temporal derivative in a boundary condition. In this section we formulate and study the principal model problem such that nonlinear problem (54) will inherit the main feature of this model problem.
Let a 0 and a ij , i, j = 1, n − 1, be some given positive constants, We look for a solution (U (x, t), ℘(x , t)) bounded at infinity and where f 0 , f 1 , f are given functions: for a positive number R 0 . Due to the quadratic form connected with the Laplace equation is positive, we can choose the coordinate x such that the form , and the following estimate holds: . Then there exists a unique solution (U, ℘) of problem (59)-(62): and then there is a unique solution (U, ℘) : , and the following estimate holds: where C i , i = 1, 3, are positive constant independent of the right-hand sides.
Note that statement (iii) of Theorem 4.1 can not be extended into arbitrary right-hand sides f 0 , f 1 , due to the spaces C k+α, k+α 3 ν are related with subdiffusion equation. Nevertheless, this result is essentially used to prove solvability of the homogenous linear problem (109)-(112) (see Section 5.1). As for statements (i) and (ii) of Theorem 4.1, they are necessary to obtain one-valued solvability of the linear problem (109)-(112) in the case of arbitrary right-hand sides.
We remark that conditions of Theorem 4.1 together with (63) lead to In virtue of Proposition 3 and Remark 1, it is enough to prove Theorem 4.1 in the case of (66).
After some simple calculations, we have: One can easily check that where W (z; b, d) is the Wright function (see (14)). Then, applying the inverse Laplace and Fourier transformations in (69), (70), we obtain Remark that by virtue of the smoothness properties of the function f (x , t) and its behavior at infinity all above performed operations are justified.

4.2.
Estimates of the constructed solution (U (x , x n , t), ℘(x , t)). As usual in the potential theory, to estimate functions ℘(x , t) and U (x , x n , t) we need to describe well the properties of kernels G(y , t), L(y , z) and K(y , x n ). Note that, the behavior of function K(y , x n ) was described in Lemma 3.1 [41]. Thus, it is necessary to evaluate the functions G(y , t) and L(y , z).
Based on the following relation from Lemma 5.3 [3]: where δ(x) is the Dirac delta function, it is easy to get equalities (78), (79) and (81). To obtain (80), we use equality (78) and statement (i) from Proposition 2 and deduce After that it is enough to apply (22) to the integral in the right-hand side in (89).
To prove statement (iv), we use (76) and conclude: Then, performing the consecutive change of variables: and (90) in the right-hand sides of (91), we deduce Note that the last inequality in (93) follows from statements (i) and (vii) of Proposition 2. As for estimate (83), we apply again inequality (76) and change of variables (92) with k = i and (90). Thus, we get The change of variable in the last integral leads to At last, applying equality (22) allows us to obtain estimate (83). As for estimate (85), it is obtained in the same way.
We will use the following results to evaluate the functions D ν t ℘(x , t) and U (x , x n , t) given with (75) and (74).
Since the proof of this result is technically tedious, we provide it in the Appendix. Then Propositions 1 and 4, Lemmas 4.2 and 4.3 together with results of Chapter 3 [25] and Lemma 3.1 [41] allow us to deduce the following results.
In addition, if f ∈ C 1+α,να/3,α (R n−1 Next step of our investigation is a proof of the corresponding estimates to the functions ℘(x , t), D ν t ℘(x , t) and U (x , x n , t) with respect to time. ). Then the following estimates hold: Proof. As for estimate (106), it follows from interpolation inequality (8), Propositions 1 and 4, and Lemma 4.3. Note that, boundary condition (60) together with inequality (106) lead to estimate of U , it is enough to use representation (102) and estimates (3.30), (3.31) from [41]. After that interpolation inequality (8) together with corresponding estimates of U The uniqueness of the constructed solution (U (x , x n , t), ℘(x , t)) follows from (64). Moreover, representation (75) and estimate (96) mean that U = K D 2 x ℘, where ℘ xixj ∈ C 2+α (R n−1 ) for ∀t ∈ [0, T ]. Thus arguments from Chapter 3 [25] allow us to show that U (x , x n , t) satisfies equation (59) with f 0 ≡ 0. Inequality (97) ensures condition (62). Hence, all the written above prove statements (ii) and (iii) of Theorem 4.1. As for validity of statement (i) of Theorem 4.1, it is necessary to show continuity of the function ℘(x , t) and U (x , x n , t) together with their derivatives with respect to time t. To this end, we repeat the arguments from Subsection 3.3 in [41] and apply statement (ii) from Theorem 4.1. Thus, the proof of Theorem 4.1 is finished.

5.
Proof of the main result: Theorem 3.1. Our proof splits into two steps.
In the first stage we show the boundedness of the linear operator A (see (54)) in the spaces C k+α,β,α . After that we prove that the nonlinear operator A −1 F is a contraction one in the corresponding classes.
Note that due to statement (i) from Proposition 3, it is is enough to prove Theorem 5.1 in the case of homogenous right-hand sides in (109) and (110).
First of all we assume that and prove existence of solution (109)-(112) in the smoother classes.
For the sake of simplicity, we consider δ = 0 in (112). Note that the case of δ = 1 will be proved in the same way.

Definition 5.2. An operator
we call the right regularizer of the operator L, if R is bounded and there is the following equality for each function F 2 ∈ C 0 1+α, 1+α 3 ν (Υ T ) Here T is a bounded operator in C 0 1+α, 1+α 3 ν (Υ T ): for enough small T .
As it follows from the principle of contracting maps, inequality (122) ensures the existence of (I + T ) −1 where I is the identity operator. Then one can conclude from (121) that the operator L −1 := R(I + T ) −1 is the right inverse operator of L and A nonlocal character of system (120) causes difficulties requiring technical tricks related to a suitable localization. As suggested in either Chapter 4 from [24] or Section 4 from Using the functions ξ m , we define the function For each m ∈ N , we pick out one point x m ∈ θ m ∩ Γ which will be the origin of the local coordinate system. The description of this system can be found in Chapter 4 [24].
We define local coordinate systems connected with each point x m , m ∈ N . Let Υ be described with y n = Ψ m (y 1 , ..., y n−1 ) in a small vicinity of every point x m , m ∈ N , and where B (m) = (β (m) ij ) i,j=1,n is an orthogonal matrix with elements β m ij , and (β m ij ) −1 is an element of the inverse matrix to B (m) . To obtain the local "flatness" of the boundary, we change the variables as z i = y i , z n = y n − Ψ m (y 1 , ..., y n−1 ), i = 1, n − 1. (124) Thus, the variables (x 1 , ..., x n ) are connected with (z 1 , ..., z n ) by the mappings Z m (see (123), (124)), such that It is easily to check that Let us introduce the function F and construct the operator R.
We define the functions u (m) (z, t), σ (m) (z , t), m ∈ N, as a solution of problem Applying statement (iii) from Theorem 4.1 to problem (130), we get the onevalued solvability of (130): Based on (130), we will say that Denote byσ After that we can define the operator R and the function u as Next we introduce the function v(x, t) as a solution of the following Dirichlet problem The definition of the operator L (see (120)) gives Then, taking into account (129)-(134), we rewrite (137) as Proposition 7. Let α, ν and Γ 0 , Υ satisfy conditions of Theorem 5.1; (117) and (118) hold. Then where for small quantities r and T .
Since the proof of Proposition 7 is standard (see, e.g., Section 4 [5]) and technically tedious, we represent it in Appendix 7.2.
Note that results of Lemma 5.3 hold in the case of a 1 i = 0, a 1 i ∈ C 2+α,ν,α (Υ T ). To show this it is enough to apply the contraction mapping theorem to (109)-(112) together with Lemma 5.3 and interpolation inequalities.
Here we used properties (115) of the functions a ij and a 1 i . Based on inequality (11), one can easily check that 6. Numerical solution of (1)- (4). In this section we employ the numerical technique from [41] to get some numerical solution of problem (1)- (4) if Ω(t) ⊂ R 2 , t ∈ [0, T ]. To this end we separate (1)-(4) into two problems. The first problem is the elliptic one defined on Ω(t): while the second problem is the initial value problem on Υ(t): Note that problem (152) is obtained from the first boundary condition in (2) and representations (5) and (32) for the "fractional" velocity of the free boundary and the unknown interface Υ(t).
Problems (151) and (152) are connected through the reciprocal usage of their solutions. That is why, the numerical approach for solving quasistationary fractional free boundary problem (1)-(4) may be constructed as the splitting scheme. Other words, at each time point problems (151) and (152) are solved separately with the mutual exchange of their solutions on the next time step.
Elliptic problem (151) is solved by the finite element method implemented using the finite element library DOLFIN [30]. At every time point t i , we construct a new mesh in accordance with solution of (152) which provides Υ(t i ). The curvature κ(Υ(t i )) is given by where (see (31)) The predictor-corrector scheme from [11] is employed to solve (152) where on the predictor stage we use an analytical approximation like (52): Thus the solving algorithm for the coupled problems (151) and (152) looks like the algorithm from Subsection 6.1 [41] and can be represented as: Algorithm Given Ω, Υ, µ, φ(y),n(ω), δ, γ, the temporal grid t i = i∆t; i ∈ 0, N , ∆t = T /N . First, we solve problem (151) and compute ∂p0 ∂n(ω) Next for each i = 1, N we perform • compute ρ pr (ω, t i ) from (153) • update Υ (and Ω) according to (31) using ρ pr (ω, t i ) • solve problem (151) and compute (155) • update Υ (and Ω) according to (31) using ρ(ω, t i ) • solve for the problem (151) and compute ∂pi ∂nt i Next we represent some example of this simulations with different ν ∈ (0, 1) and the same initial spatial domain Ω. Let Ω be restricted by two elliptic boundaries: Numerical experiments are performed for three different values of ν, namely, ν 1 = 0.25, ν 2 = 0.5, ν 3 = 0.75. We analyze the dynamics of the free boundary for different values of ν. As it follows from Fig. 1, smaller values of ν require longer time to change the boundaries.
Moreover, Fig. 2 demonstrates that ν controls dynamics of the moving boundary Υ(t) in the following way. There is a critical time T * := T * (ν), such that for t < T * the free boundary moves rapid if ν decreases. However, the effect is opposite for t > T * , i.e. the free boundary moves slowly for small ν and much faster for ν approaching 1. As it follows from Fig. 2, T * is located in the neighborhood T = 0.75. Drawing the shape of free boundaries for t near T = 0.75 and ν = ν i , i = 1, 2, 3, (see Fig.3 in the case of ν 2 = 0.5, ν 3 = 0.75) we find critical time: • for l, j = 1, n. Hence, representation (102) is an easy consequent of representation (156) if l = n. As for case l = j = n, one can check with simple calculations that After that we integrate by parts in (156), l = j = n, and taking into account relation (157) together with estimate (3.30) [41] to get representation (102).  Next, we apply equality (88) to representation (73) and get condition (60). Thus, to finish the proof of Proposition 4, it is necessary to show that the functions U (x , x n , t) and ℘(x , t) represented by (74) and (75) satisfy condition (61). :  After that we apply (88) to the second term in (158) and (16) to the first one in (158) and obtain Then, we return to representation (101) and calculate ∂ ν τ G(y , τ ). Proposition 2, equality (81) together with integrating by parts lead to the relation Therefore, condition (61) follows from equalities (101), (159) if r and T are enough small.
Note that the rest terms in the left-hand side of (167) are estimated with the same way. Hence, inequalities (166) and (167) ensures the corresponding estimate of the first term in the right-hand side of (162) n i,j=1 w xixj (α) x,Ω T ≤ c 2 [r + ae + T . (172) To evaluate the second term in the right-hand side of (162), we use interpolation inequality (8) and deduce