BDF Schemes in Stable Generalized Finite Element Methods for Parabolic
Interface Problems with Moving Interfaces

There are several difficulties in generalized/extended finite element methods (GFEM/XFEM) for moving interface problems. First, the GFEM/XFEM may be unstable in a sense that condition numbers of system matrices could be much bigger than those of standard FEM. Second, they may not be robust in that the condition numbers increase rapidly as interface curves approach edges of meshes. Furthermore, time stepping schemes need carrying out carefully since both enrichment functions and enriched nodes in the GFEM/XFEM vary in time. This paper is devoted to proposing the stable and robust GFEM/XFEM with efficient time stepping schemes for the parabolic interface problems with moving interfaces. A so-called stable GFEM (SGFEM) developed for elliptical interface problems is extended to the parabolic interface problems for spatial discretizations; while backward difference formulae (BDF) are used for the time stepping. Numerical studies demonstrate that the SGFEM with the first and second order BDF (also known as backward Euler method and BDF2) is stable, robust, and achieves optimal convergence rates. Comparisons of the proposed SGFEM with various commonly-used GFEM/XFEM are made, which show advantages of the SGFEM over the other GFEM/XFEM in aspects of stability, robustness, and convergence.


Introduction
Numerical simulations of multi-material physical phenomena have become increasing important aspects in computational science and engineering. The multi-material problems can be modelled by various partial differential equations with interfaces. The interfaces are moving in time in the real world and engineering applications, such as fluid-structure interaction, multiphase flows, free boundary problems, flows in porous media, 3D fiber composites. See [1][2][3][4] and references therein. Conventional numerical methods [5][6][7], such as finite element methods (FEM), have to update or refine meshes to match the interface movement at every time step to get approximation accuracy. However, the mesh updating or refining fitted to the interfaces can be extremely difficult, especially in 3D time-dependent problems.
Much attention has been paid to develop unfitted FEM for the moving interface problems in the last decades. The meshes in the unfitted FEM are simple, fixed, and independent of the interface movements. The typical unfitted FEM includes immersed FEMs [8][9][10], discontinuous Galerkin methods [11,12], cutting FEMs [13][14][15][16], Nitsche extended FEMs [17,18]. These methods are nonconforming, and various stabilized or penalty techniques were proposed to control conforming errors. The penalty parameters need determining with care and cannot be treated using a unified approach. We focus on the conforming unfitted FEM for the moving interface problems in this paper. Generalized or extended FEMs (GFEM/ XFEM) turn to be prominent approaches to construct conforming unfitted shape functions. The GFEM/ XFEM augments the standard FEM with special functions that mimic local features of exact solutions to solve complicated non-smooth engineering problems. These special functions are "pasted" together using techniques of partition of unity (PU) [4,19]. Meshes in the GFEM/XFEM are typically simple, fixed, and independent of the non-smooth features of problems. The GFEM/XFEM has been extensively applied to a wide range of engineering problems, e.g., crack growth, material modelling, multiphase flows, and fluid-structure interaction. We refer to the review articles [20][21][22][23][24] for various aspects of GFEM/XFEM. The GFEM/XFEM was proposed to simulate elliptical and moving interface problems quite early; see [25][26][27][28][29][30][31][32][33][34][35][36][37] for examples. The GFEM/XFEM can be viewed as an instance of the partition of unity methods (PUM) [4,19]. In the rest of paper, we will use the term GFEM instead of GFEM/XFEM.
There are several essential difficulties in applications of the GFEM to the moving interface problems. In aspect of spatial approximations, almost linear dependence between the FEM shape functions and enriched special function causes stability issues, i.e., scaled condition numbers (SCN) of stiffness matrices of the GFEM are much higher than those of the standard FEM. On the other hand, the SCN rapidly increases as the interface curves get close to boundaries of elements. It gives rise to robustness issues of the GFEM. The bad stability and robustness may cause disastrous round-off errors in elimination methods or the slow convergence of iterative schemes to solve underlying linear systems. Many stability and robustness techniques have been proposed for the GFEM, such as, locally adapting positions of either nodes or interface curves to balance the ratios of the two volumes created by the intersections of the interfaces with the elements [35,38], preconditioning the stiffness matrices or orthogonalization [39][40][41][42][43], and correcting the enrichments by interpolant [25,26,31,[44][45][46][47][48][49].
In aspect of temporal discretization, both basis functions and numbers of enriched degrees of freedom of the GFEM/XFEM for the moving interface problems change in time so that conventional time integration schemes for the standard FEM can not be directly applied to the GFEM. Time-space element methods for the GFEM were studied in [18,32,36,50,51] for the moving interface problems. The optimal convergence was reported in [36], and a theoretical analysis on the convergence was executed in [18]. The space-time GFEM suffers from a substantial increase of computational cost since the time-space element methods with one dimension higher than the spatial dimensions are used. Therefore, efficient time-stepping schemes are very important for the GFEM of the moving interface problems [23,33,36]. We briefly analyze major difficulties of time stepping in the GFEM. Let u(t) and u h ðtÞ ¼ P j2J c j ðtÞn j ðtÞ be the exact and approximate solutions of moving interface problems, respectively, where h is a discretization parameter, ξ j 's are the GFEM shape functions. We note that ξ j (t) changes with time. The derivative of u h (t) with respect to t is When stepping schemes for (u h ) t are employed, appearance of (ξ j ) t causes extra matrices composed of ξ j and (ξ j ) t , in addition to stiffness and mass matrices, see [10] for instance. Moreover, in the GFEM, we have ðP h uÞ t 6 ¼ P h u t ; (2) where P h is an elliptical projection of u to the GFEM subspace, see [17,52]. It is known that the equality of (2) is crucial for theoretical analysis on discretization errors in space and time [52]. The features (1) and (2) of the GFEM cause complexities in employment of time stepping [23,36]. This is very similar with situations of adaptive FEM for evolving partial differential equations [53,54]. The time stepping schemes in the GFEM can be found out in [29,33,35,36,55,56]. The optimal convergence rate of a GFEM with backward Euler methods was proven in [17]. To our best knowledge [17] is the first theoretical analysis on time stepping of the GFEM, when applied to the moving interface problems.
Recently, a stable GFEM (SGFEM) was proposed to address the stability and robustness of GFEM for elliptical interface and crack problems in [25,26,[44][45][46][47][48][49][57][58][59]. The SGFEM (i) yields the optimal convergence rates, (ii) is stable in that its SCNs are of same order as those of FEM, and (iii) is robust, i.e., the conditioning does not worsen as the interface approaches the boundaries of elements. In this paper we carry out a numerical study on time stepping in the SGFEM for parabolic interface problems with moving interfaces. The SGFEM developed for the elliptical problems is extended to the parabolic problems for spatial discretizations. The time stepping schemes are derived by first discretizing time other than space because the shape functions and enrichment schemes in SGFEM are all time-dependent. The backward difference formulae (BDF) [60,61] are proposed for the time stepping in this paper. The first and second order BDF are referred to as the backward Euler method and BDF2 in the literature, which typically produce the first and second accuracies in time. It is shown in the paper that the SGFEM with the backward Euler method and BDF2 reach optimal convergence orders for mild and large contrast parameters. We compare the proposed SGFEM with other commonly-used GFEMs, including topological, geometric, and corrected GFEMs, and the stability and robustness advantage of SGFEM over these methods are shown in the paper. The optimal convergence of SGFEM based on the backward Euler method and BDF2 will be proven in a forthcoming study. This paper is organized as follows. The model problems and their conventional GFEM and SGFEM are presented in Sections 2 and 3, respectively. The BDF schemes are proposed in Section 4. The stability, robustness, and convergence of the various GFEMs and SGFEM are analyzed and compared in Section 5. The concluding remarks are presented in Section 6.

Moving Interface Problem
Let & R 2 be a bounded and convex domain with piecewise smooth boundary ∂Ω, and we consider the parabolic moving interface problem with a Neumann boundary condition @u @t À r Á ðaruÞ ¼ f ðP; tÞ; P 2 ; t 2 ½0; T ; a @u @ñ ðP; tÞ ¼ gðP; tÞ; P 2 @; t 2 ½0; T ; uðP; 0Þ ¼ u 0 ðPÞ; P 2 ; subject to jump conditions on a moving interface Γ(t) whereñ andñ À denote the unit outward normal vectors to the boundary ∂Ω and the interface Γ(t), respectively. The interface Γ(t), t ∈ [0, T], divides Ω into two sub-domains Ω 0 (t) and Ω 1 (t) such that ¼ 0 ðtÞ [ 1 ðtÞ, Ω 0 (t) ∩ Ω 1 (t) = [, and À ¼ 0 ðtÞ \ 1 ðtÞ, see Fig. 1 for an illustration. The interface Γ(t) varies with time t, and is smooth for any t ∈ [0, T]. The notation v ½ À is jump of a quantity v on the interface Γ, which is defined as the difference on Γ of v values, limited to 1 ðtÞ and 0 ðtÞ. For a point P ¼ ðx; yÞ 2 R 2 , we assume that the diffusion coefficient a(P) is a piecewise-constant function defined by aðP; tÞ ¼ a 0 ; P 2 0 ðtÞ; a 1 ; P 2 1 ðtÞ; where a 0 , a 1 are two positive constants. The contrast ρ of the diffusion coefficients is defined by It is known from the interface condition (4) that the solution u to (3) is continuous, but its gradients are not continuous on account of the discontinuity of the coefficient a(P) on the interface Γ. Such a property is frequently termed as a weak discontinuity [23,34].
Then we pose the equivalent variational form of the parabolic interface problem (3)

GFEM and SGFEM
The conventional FEM has to update or refine meshes to match the interfaces since the interfaces move in time. However, the mesh updating or refining can be extremely difficult. The GFEM augments standard FEM with non-polynomial shape functions that are derived based on a priori knowledge about the solution of underlying variational problem. The meshes in the GFEM are fixed, simple, and independent of interface movements, such as Cartesian meshes. Hence, the mesh updating or refining can be avoided.
To describe the GFEM, we start with a uniform finite element mesh T h ¼ fe s g, characterized by a mesh-size parameter h. The element e s can be triangular or quadrilateral, and the mesh T h is fixed and independent of the moving interfaces. The element e s is closed. Denote {P i }i ∈ I h to be the set of nodes associated with the mesh T h , where I h is the index set of nodes. For any i ∈ I h , let ω i be a patch with respect to P i , which is the union of all elements sharing the common node P i , namely, It is clear that ω i is a closed set. Let f i be the standard linear (bilinear for the quadrilateral elements) FE hat function, associated with the node P i with supp{f i } = ω i . Subordinate to the cover where C > 0 is independent of i, h.
The standard FEM is used for the spatial approximation of the variational (6), and its approximation space is defined as follows: As discussed above, the FEM (8) cannot approximate the solution of (6) with optimal convergence rates since the mesh does not match the interfaces.
The approximation space of GFEM is generally defined as where S t ENR is called enrichment part of GFEM, Å t i are referred to as enrichment functions that mimic the discontinuities in the solution of underlying variational problem, and I t h;enr denotes the set of enriched nodes. All S t ENR , Å t i , and I t h;enr may change with time t. It is seen from (9) that the GFEM has great potentials to improve the approximation accuracy because of introduction of the enrichment part S t ENR . For the interface problems with the weak discontinuities, the enrichment function Å t i can be taken as the absolute value of level set functions or the following distance functions D(P,t) [23,31,45]: We now describe various forms of GFEM, applied to the interface problem.
Topological GFEM: Define an index set of the enriched nodes and the topological GFEM [23,45] is defined as The topological GFEM has a minimal number of the enriched nodes, i.e., CardfI À;t h;enr g :¼ Oðh À1 Þ; (11) see Fig. 2 right plot and Fig. 9 right plot, where Card{X} is the dimensionality of a set X. It was reported in [45] the topological GFEM only produces a sub-optimal convergence error Oð ffiffi ffi h p Þ other than the optimal O(h).
Geometric GFEM: The optimal convergence can be restored by increasing the number of enriched nodes. The geometric GFEM [23,45,48] enriches the nodes in a fixed neighborhood of the interface Γ(t), namely, h;enr . The approximate subspace of the geometric GFEM is as follows: The geometric GFEM gives rise to the optimal convergence error O(h) for elliptical interface problems [45]. However, the geometric GFEM results in much more enriched degrees of freedom (DOF), in which Furthermore, the introduction of more enriched DOFs causes bad conditioning O(h −4 ) for stiffness matrices, which is much bigger than that of the standard FEM [45].
Corrected GFEM: Effects of blending elements were analyzed in [23,62], and a corrected GFEM was proposed there to improve the approximation accuracy of GFEM. The enrichment strategy of corrected GFEM at each time t is defined as The following approximate subspace is used in the corrected GFEM [23,62]: where η(P) is called a ramp (or cutoff) function defined as h;enr f i ðPÞ: Apparently, the corrected GFEM has a number O(h −1 ) of the enriched nodes that is slightly larger than topological GFEM (see Fig. 2 right plot and Fig. 9 right plot), and we have The corrected GFEM is able to achieve optimal convergence error O(h). However, the corrected GFEM may not be robust with respect to the mesh because its conditioning may blow up when the interface Γ(t) gets close to the mesh line, as shown in numerical experiments.
Stable GFEM: To address the ill conditioning of GFEM and the lack of robustness, a simple local modification of subtracting the interpolant of the original enrichment functions was proposed in [25,26,[44][45][46][47][48], recently. The modified GFEM is referred to as stable GFEM (SGFEM). The approximate subspace of the SGFEM for the interface problems at each time t is defined as where T h f is the FE interpolant of a continuous function f, based the preliminary FE hat functions. Clearly, the enriched nodes in SGFEM are the same as in the topological GFEM. It was proved [45] that the SGFEM (17) (a) yields the optimal convergence order O(h), (b) has an order of SCN around O(h −2 ) that is of same order as that of FEM, and (c) is robust with respect to the mesh. These developments in the SGFEM are achieved for the elliptical interface problems, and extensions to the moving interface problems have not been made yet. In next section we propose time stepping schemes for the SGFEM applied to the moving interface problems, which satisfy the features (a)-(c).
Remark 3.1: We analyze the calculation complexity of SGFEM and make comparisons with the other methods. The mesh in SGFEM (also the other GFEMs) is simple, fixed, and independent of the moving interfaces. Therefore, the SGFEM enjoys significant computational advantages over the FEM because the mesh in FEM has to be updated at every time step to match the evolving interfaces. Compared with the geometric and corrected GFEM, the SGFEM possesses the minimal number of enriched nodes, see (16). Although the topological GFEM has the same number of enriched nodes as that of SGFEM, it cannot derive optimal convergence errors, as illustrated in the numerical experiments. The stiffness matrix of SGFEM is slightly bigger than that of FEM because of introduction of enriched functions. The number of enriched nodes in SGFEM is O(h −1 ) (11), which is one dimension less than the number of FE nodes, O(h −2 ). Therefore, the computational costs in enriched part are negligible for small h. Moreover, this is quite trivial considering the merit of non-matching meshes in the SGFEM.

BDF Schemes in GFEM and SGFEM
Conventional time-space element methods and time stepping schemes have been proposed in fields of the GFEM, see [23] for a review. The time-space GFEMs for the moving interface problems [18,32,36,50,51] suffer from a substantial increase of computational cost since the time-space element methods with one dimension higher than the spatial dimensions are used. Therefore, time-stepping schemes with high accuracy are very important for the GFEM when applied to the moving interface problems. There are two kinds of commonly used time stepping schemes in the standard FEM: (a) if the shape functions are independent of time, spatial discretization is first carried out to derive systems of ordinary differential equations with respect to coefficients of the shape functions, which are then solved by finite difference methods [52]; (b) if the interfaces evolve with time, moving meshes can be considered in arbitrary-Lagrangian-Eulerian (ALE) frameworks, and problems of shape functions at different time levels can be solved using mesh velocities [52,63]. However, these two approaches can not be directly applied to the GFEM for moving interface problems because, typically, the shape functions of GFEM are timedependent, and the meshes are fixed relative to the interface movements. The objective of this paper is to propose efficient time stepping schemes for the SGFEM of moving interface problems, and the resultant GFEMs are stable and robust with regard to the interface movements, and has optimal convergence rates.
We analyze which of spatial and temporal discretizations are first carried out. In the presence of moving discontinuity of the solution u, the approximation subspaces S t h of GFEM and SGFEM are time-dependent because the enrichment functions D(P,t) and the enrichment schemes I t h;enr are all time-dependent. Let fn i ðtÞg & S t h be the basis of GFEM or SGFEM, then ξ i (t) and their coefficients are all time-dependent, see (1). If the space is first discretized, the dependence of ξ i (t) on t yields blending matrices of ξ j and (ξ j ) t , in addition to stiffness and mass matrices, see [10] for instance. This increases computational complexity in assembling and solving linear systems. Therefore, it is feasible to first discretize time other than space in the GFEM and SGFEM, as discussed in [17,36].

Temporal Discretization
In this paper, we mainly focus on backward difference formulae (BDF) [60,61] with the coefficients δ j given by They are used to generate kth-order time semi-discrete schemes, based on the weak formulation (6).

Spacial discretization
The semi-discrete scheme (20) are approximated by the various GFEMs and SGFEM presented in Section 3. In order to guarantee symmetric and positive definite system matrices, the test and trial function space should be chosen at the same time step (the current time step t n ). For any u n h 2 S t n h , 0 ≤ n ≤ N, let d j u nÀj h ; for 1 n N ; and based on the various GFEMs and the semi-discrete scheme (20), we suggest fully discrete formulation for parabolic problem with moving interface as follows, where u 0 h is an approximation to u 0 . Remark 4.1: In this paper, we only consider the first and second order accuracies in time, i.e., k = 1, 2 in (21) since we construct the GFEM or SGFEM using the linear (or bilinear) elements, which produce the first and second order accuracies in space for H 1 and L 2 errors, respectively. The cases k ≥ 3 can be treated similarly, and we do not present them here. Taking k = 1, 2 gives us and Find u n h 2 S t n h ; 8 1 n N; such that respectively. The formulations (22) and (23) are referred to as backward Euler and BDF2 methods in the literature, respectively. We mention that in the BDF2 we take u À1 h ¼ u 0 h so that calculating u 1 h of the first time step t 1 is reduced to the backward Euler method with time step length 2 3 s. The various GFEMs and SGFEM with the backward Euler method (22) and BDF2 method (23) are analyzed numerically in next section, where their convergence, stability, and robustness are examined and compared.
In the end of this section, we define the scaled condition number (SCN) of stiffness and mass matrices of GFEMs, which will be computed below to show the conditioning and robustness of various GFEMs and SGFEM. Denote a stiffness or mass matrix by A, associated with the approximation space of GFEM S t h ¼ S FEM È S t ENR , and it is of the form where A 11 and A 22 are the sub-matrices associated with the FEM part S FEM and enrichment part S t ENR of GFEM respectively, while A 12 is the sub-matrix associated with the cross-terms between S FEM and S t ENR .
Consider a diagonal matrices D with D ii ¼ A À1=2 ii , and defineÂ ¼ DAD. The scaled condition numbers (SCN) of A, A 11 are defined by respectively, where κ 2 (·) is the condition number based on the Euclidean vector norm || · || 2 .
In next section, we will show numerically that the proposed SGFEM with backward Euler method and BDF2 achieves optimal convergence orders, and is stable and robust in that its SCN K is of same order as that of FEM, and K will not blow up as the interfaces are close to the edges of mesh. The geometric and corrected GFEMs are not stable, and the topological GFEM only has sub-optimal convergence. The stability and robustness of GFEM are relevant for the moving interface problems because linear systems are solved frequently according to various possible relative positions of interfaces and edges of meshes.

Numerical Studies and Discussion
In h;enr . Integration for Nonsmooth Enrichments: We describe the numerical integration formulae for the moving interface problems. For an element not cut by the interface, the standard 4 × 4 Gaussian rule is adopted. For an element passed by the interfaces of one or two time layers, we first replace the interface curves with straight lines by connecting intersection points of the interfaces and the boundaries of element. In this way, the element is divided into several sub-polygons, each one of which is decomposed into triangles, and 6 × 6 Gauss quadrature rule in each triangle is employed. Such an integration rule is easy to implement and commonly used in the GFEMs for moving interface problems. We refer to [23,41] for more discussions about the numerical integration of non-smooth enrichments.
We will compute the relative errors at the final time step t N in L 2 and energy norms and SCNs to compare and analyze the various GFEMs and SGFEM with backward Euler and BDF2 methods. The SCN of stiffness or mass matrix is defined in (25), and the relative errors are defined as follows: where u N h and u(·,T) are the numerical and exact solutions at the final time point, respectively. In all the computations below, we test two contrast parameters, mild one ρ = 100 (a 0 = 100, a 1 = 1) and large one ρ = 1000 (a 0 = 1000, a 1 = 1). The uniform M × M meshes with M = 2 i+1 , i = 1, 2, …, 7 are employed.
We consider a manufactured solution of (3) given by u ¼

Convergence Studies
We first take τ = h, and the log-log plots of L 2 and energy errors with respect to h for the various GFEMs and SGFEM with the backward Euler and BDF2 methods are presented in Figs. 3 and 4, respectively. For the backward Euler scheme, a L 2 convergence error O(h 2 + τ) and an energy convergence error O(h + τ) are expected. Since we take τ = h, these two errors are all O(h). Fig. 3 Left shows the L 2 error O(h) for all GFEMs and SGFEM, whereas the energy optimal error O(h) is only obtained by the SGFEM and GFEM Cor in Fig. 3 Right. The GFEM Topo and GFEM Geo based on the backward Euler scheme can not produce the optimal energy convergence errors. The optimal spatial convergence error O(h 2 ) in L 2 norm is not obtained for the various GFEMs and SGFEM because the backward Euler scheme is of first order in time (remembering that τ = h is taken). The optimal L 2 spatial convergence error O(h 2 ) can be observed in the BDF2 method since it produces the time discretization of second order. For the BDF2, the optimal L 2 error O(h 2 + τ 2 ) and energy error O(h + τ 2 ) are expected. Therefore, these error are O(h 2 ) and O(h), respectively, since we take τ = h. It is pointed in Fig. 4 that the SGFEM and GFEM Cor with the BDF2 yield optimal convergence errors O(h 2 ) in L 2 and energy norms O(h) for both mild and large contract parameters. Again, these optimal rates can not attained by the GFEM Topo and GFEM Geo .
In sum, the GFEM Topo and GFEM Geo can not reach the optimal convergence rates for the BDF schemes. The SGFEM and GFEM Cor , based on the backward Euler and BDF2 schemes, converge with the optimal rates in both L 2 and energy norms for mild and large contrast parameters.
We note that even though both the proposed SGFEM and GFEM Cor exhibit excellent convergence, the SGFEM has fewer DOFs than the GFEM Cor . Moreover, we will see that the SGFEM is stable and robust, compared with the GFEM Cor , in the following conditioning studies.

Conditioning of Stiffness/Mass Matrices and Robustness Test
We first lay an emphasis on the SCN of stiffness/mass matrices of various methods, which indicates the stability of GFEMs. The SCNs K A of stiffness matrices with respect to h at time t = 0.25 and t = 0.50 are plotted in Fig. 5. It is clear that growths of K A for GFEM Topo , GFEM Cor , and SGFEM are of order O(h −2 ), which is of the same order as that of the standard FEM. However, K A of GFEM Geo grows with an order about O(h −4 ), which is much bigger than that of FEM. This is caused because of introduction of too many enriched DOFs in GFEM Geo .
The SCNs K M of mass matrices with respect to h at time t = 0.25 and t = 0.50 are shown in Fig. 6. It is noted that K M of GFEM Topo , GFEM Cor , and SGFEM are of order O(1) that is the same as that of the standard FEM, while the growths of K M in GFEM Geo are of order O(h −4 ), which is significantly higher than other methods. We next study the robustness of various GFEMs when the interfaces approach the boundaries of elements. To this end, we consider an interface Γ represented by equation y = 0.5 + δ, and the discontinuous coefficient a(P) is set by aðPÞ ¼ a 0 ; y ! 0:5 þ d; a 1 ; y , 0:5 þ d: The domain Ω = [0,1] 2 is divided with a fixed 32 × 32 uniform mesh, and the enrichment schemes of various GFEMs are displayed in Fig. 7. We investigate the conditioning of various GFEMs and SGFEM as the interface Γ approaches the mesh line y = 0.5, namely, δ decreases to zero. We take The SCNs of the stiffness and mass matrices against δ above with ρ = 100 and ρ = 1000 are plotted in Fig. 8. It shows that the FEM, SGFEM and GFEM Topo are robust in a sense that their SCNs K A and K M do not vary as δ decreases. However, K A and K M of GFEM Cor cause massive increases as δ decreases. Therefore, it indicates that GFEM Cor is not robust, while the SGFEM remains robust with respect to the relative position of interface and boundaries of elements. We emphasize that the robustness is relevant for the GFEMs of parabolic moving interface problems, because various relative positions of interfaces and meshes may occur when the meshes are fixed and the interfaces move. Detailed studies on the robustness of GFEMs are relatively new, we refer to [26,45,47,49] for more details. These numerical results illustrates that among the GFEMs tested in this section, the SGFEM is the only one who both achieves the optimal convergence orders and is stable and robust.

More Relations between τ and h
In the subsection 5.1.1 we verify the optimal L 2 error O(h 2 + τ) and energy error O(h + τ) of backward Euler method, and the optimal L 2 error O(h 2 + τ 2 ) and energy error O(h + τ 2 ) of BDF2 scheme by taking τ = h.  In this subsection we further test these optimal rates for the SGFEM by taking more relations between τ and h. We keep all the settings in the subsection 5.1.1, including the domain, equation, exact solution, mesh, and so on. We take τ = 4h 2 for the backward Euler method and s ¼ ffiffi ffi h p for the BDF2 scheme. The L 2 errors of SGFEM with backward Euler method and the energy errors of SGFEM with BDF2 scheme against h ¼ 1 2 iþ1 ; i ¼ 1; 2; Á Á Á ; 6 are given in Tab. 1. It can be seen in Tab. 1 that the optimal spatial L 2 error of SGFEM with backward Euler method O(h 2 ) is derived. This coincides with the optimal L 2 error O(h 2 + τ) of backward Euler method because of τ = 4h 2 . We mention that the spatial L 2 error O(h 2 ) is not observed in the subsection 5.1.1. In Tab. 1 we also see the optimal spatial energy error O(h) of SGFEM with BDF2 scheme. This again verify the optimal energy error O(h + τ 2 ) of SGFEM with BDF2 scheme since s ¼ ffiffi ffi h p .

An Interface Problem with Curved Interfaces
We next show that the behaviors of convergence and conditioning of SGFEM in the case of straight interfaces are maintained when the interfaces are curved. We consider a domain Ω = [−1,1] 2 with an  = (x − x c ) 2 + (y − y c ) 2 − r(t) 2 = 0 is the circle of center O(x c ,y c ) and time-dependent radius r(t). Let Ω 0 : = {P ∈ Ω: γ(P,t) > 0 } and Ω 1 : = {P ∈ Ω: γ(P,t) < 0} as shown in Fig. 9 Left. We use (x c ,y c ) = (0,0) and rðtÞ ¼ r 0 Á sin tþ3 4 with r 0 ¼ p 6:28 in the computations. Besides, the parameter of enriched region is taken as R ¼ 1 3 rðtÞ at time t for the GFEM Geo . In Fig. 9 Right, the enrichment schemes of various GFEMs are displayed at time t = 1, where R ¼ 1 3 rð1Þ for the GFEM Geo .
We consider a manufactured solution of (3) given by where the right hand side f of (3), the boundary function g, and the initial value u 0 are obtained from the exact solution u using Eq. (3). Fig. 10 Left shows that the SGFEM with the backward Euler scheme produce the optimal convergence errors in L 2 and energy norms for the mild and large contrast parameters (the L 2 error is O(h) because of τ = h). Likewise, Fig. 10 Right demonstrates that the SGFEM based on the BDF2, a time discretization of second order, is of optimal convergence for the mild and large contrast parameters, and the errors in L 2 and energy norms are O(h 2 ) and O(h), respectively. These observations agree with the situations of moving straight interfaces. Therefore, the SGFEM based on the backward Euler and BDF2 schemes show the optimal convergence and robustness with the high contrast parameters, which are advised to applied to the space and time discretizations of parabolic moving interface problems.

More Relations between τ and h
As in the subsection 5.1.3, we test the optimal convergence rates for the case of curved interfaces by taking more relations between τ and h. We take τ = 4h 2 for the backward Euler method and s ¼ ffiffi ffi h p for the BDF2 scheme. The L 2 errors of SGFEM with backward Euler method and the energy errors of SGFEM with BDF2 scheme against h ¼ 1 2 iþ1 ; i ¼ 1; 2; . . . ; 6 are given in Tab. 2. We see the optimal spatial L 2 error O(h 2 ) of SGFEM with backward Euler method and the optimal spatial energy error O(h) of SGFEM with BDF2 scheme in Tab. 2. These again account for the optimal L 2 error O(h 2 + τ) of backward Euler method (τ = 4h 2 ) and the optimal energy error O(h + τ 2 ) of BDF2 scheme (s ¼ ffiffi ffi h p ), respectively.

Conditioning of Stiffness and Mass Matrices
We present the SCNs of SGFEM against h with the contrast parameters ρ = 100, 1000 at time points t = 0.25 and t = 0.50 in Fig. 11. It is clear that the SCNs of stiffness and mass matrices grow with orders around O (h −2 ) and O(1), respectively, which are the same as those of FEM. Moreover, the SCNs are independent of the relative positions between interfaces and boundaries of elements. These results are in accordance with the case of moving straight interfaces. This means that the SGFEM is indeed a stable and robust GFEM for the moving interface problems.

Conclusions
The SGFEM is shown to be stable and robust for the elliptical interface problems in [25,26,44,45], namely, its SCNs are of the same order as those of FEM, and do not blow up as the interfaces approach boundaries of elements. The stability and robust of SGFEM are of significantly importance for moving interface problems where various possible relative positions between interfaces and meshes may occur because the meshes are fixed, and the interfaces are moving. Since the shape functions of SGFEM are time-dependent, the time stepping schemes need carrying our carefully. In this paper, we proposed the BDF schemes for SGFEM when applied to the parabolic moving interface problems. The first and second order BDF (the backward Euler method and BDF2) are employed to establish the first and second order schemes in time, respectively. The numerical experiments indicated that the SGFEM with the backward Euler method and BDF2 produces the optimal convergence rates for mild and large contrast parameters. We compared the SGFEM with other conventional GFEMs, such as topological, geometric, and corrected GFEMs in aspects of convergence and conditioning. The SGFEM exhibited its advantages of convergence, stability, and robustness over these GFEMs for the parabolic moving interface problems indeed. The optimal convergence of SGFEM with the backward Euler method and BDF2 will be proven in future. The SGFEM with certain time stepping schemes for nonlinear moving interface problems will also be an interesting direction. Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.