## Abstract

Anomalous diffusion can be characterized by a mean-squared displacement 〈*x*^{2}(*t*)〉 that is proportional to *t*^{α} where *α*≠1. A class of one-dimensional moving boundary problems is investigated that involves one or more regions governed by anomalous diffusion, specifically subdiffusion (*α*<1). A novel numerical method is developed to handle the moving interface as well as the singular history kernel of subdiffusion. Two moving boundary problems are solved: the first involves a subdiffusion region to the one side of an interface and a classical diffusion region to the other. The interface will display non-monotone behaviour. The subdiffusion region will always initially advance until a given time, after which it will always recede. The second problem involves subdiffusion regions to both sides of an interface. The interface here also reverses direction after a given time, with the more subdiffusive region initially advancing and then receding.

## 1. Introduction

A system undergoing classical diffusion will display a mean-squared displacement 〈*x*^{2}(*t*)〉 that is linearly proportional to time. However, there are a variety of processes that will instead display mean-squared behaviour of 〈*x*^{2}(*t*)〉∼*t*^{α} where *α*≠1. These processes are said to be undergoing anomalous diffusion and can be modelled by an integro-differential equation known as the fractional diffusion equation. Here, a class of one-dimensional moving boundary problems is investigated where one or both regions are governed by anomalous diffusion. In some instances, the dynamics will parallel the analogous classical diffusion result, whereas in others the effect of anomalous diffusion on the interface is substantial, resulting in novel non-classical behaviour.

Metzler & Klafter (2000, 2004) present a variety of physical situations where anomalous diffusion occurs, including certain flows in porous media and conduction in amorphous semiconductors. Recently, Gusev *et al.* (1993) and Gusev & Suter (1993) observed anomalous behaviour in simulations of gas in polymer glasses. Eaves & Reichman (2008) explored subdiffusive exit-time behaviour within biomolecular enzyme kinetics. Additional biological studies include that of anomalous electrodiffusion by Langlands *et al.* (2009), anomalous diffusion of morphogens by Yuste *et al.* (2010) and anomalous diffusion of proteins through trehalose-water solutions by Schlichter *et al.* (2001). In addition to proteins, Conrad & de Pablo (1998) find that water molecules follow anomalous diffusion in trehalose-water solutions, albeit only for short periods of time. Hence, subdiffusion is associated with numerous physical phenomena. While there has been considerable work directed at the dynamics within a single, fixed anomalous region, less studied is the moving boundary problem where an anomalous region is in contact with a classical region or another anomalous region. This type of moving boundary problem is investigated here. The governing equations are determined by making general assumptions about wait-time and jump-length distributions in each region. The boundary conditions at the interface are chosen by considering a glass formation front; however, other boundary conditions can be substituted depending on the specific physical application.

Moving boundary problems involving classical diffusion have been studied extensively over the last couple of centuries. Wilson (1978) provides a review of the analyses of such Stefan moving boundary problems. However, anomalous diffusion requires consideration of the fractional diffusion equation that governs subdiffusion. An extensive summary of how fractional calculus is used to formulate the fractional diffusion equation is provided by Metzler & Klafter (2000). Various Green's function approaches to solving the fractional diffusion equation on infinite and semi-infinite domains are also suggested. However, none of these approaches addresses the presence of boundaries. Kosztolowicz (2004) introduces a fixed boundary and solves the two-sided subdiffusion system by mapping the corresponding Green's function for the two-sided classical diffusion system. This approach is then extended to solve systems with both subdiffusion and classical diffusion. Grebenkov (2010) studies a bounded domain with a partially absorbing–reflecting boundary. He is able to make estimates on exit-time behaviour using spectral decomposition. However, both of these approaches do not readily apply to moving boundaries as they use time-independent eigenfunctions. Liu & Xu (2004, 2009) consider a moving boundary in problems involving subdiffusion only, where a similarity variable is available, but their approach is not applicable for the system with both subdiffusion and classical diffusion present. Additionally, with both subdiffusion and classical diffusion present, the fractional derivative must be adjusted so that the history extends back only within the subdiffusion region. It is shown later that this situation will change the form of the anomalous flux depending on the direction of propagation of the interface. Thus, a numerical method is developed here in lieu of the aforementioned analytic methods.

The fractional diffusion equation, specifically the fractional derivative, provides an interesting numerical challenge. Yuste & Acedo (2005) focus on the Grunwald–Letnikov (G-L) fractional derivative that is defined as a limit similar to the definition of the integer-order derivative. As suggested by Podlubny (1999), the G-L derivative is approximated by taking the limit variable close to zero. Although this discretization is first order, it is not straightforward to extend this to higher orders if needed. Oldham & Spanier (1974) consider the Riemann–Liouville form of the fractional derivative that is defined in terms of a convolution integral. They propose approximating the fractional derivative using a combination of integration quadratures and explicit integration. Unlike the approach by Yuste & Acedo, this approximation can be easily extended by using higher-order integration quadratures. Langlands & Henry (2005) use the approximation by Oldham & Spanier (1974), along with a forward Euler finite-difference scheme, to discretize the fractional diffusion equation. An alternative method is suggested by Diethelm *et al.* (2002). They convert a generalized version of the fractional diffusion equation to a Volterra integral equation and then use the Oldham & Spanier (1974) approximation to evaluate the integral. Increasing the order here requires only higher-order integration quadratures, which is an advantage over the approach by Langlands and Henry that requires a higher-order replacement for the forward Euler finite difference as well.

The numerical method that is developed here to solve the moving boundary problems is based on the method of Diethelm *et al.* (2002). The fractional diffusion equation is converted into the corresponding Volterra integral equation and discretized analogous to the Oldham & Spanier (1974) approximation. The moving boundary requires the addition of a front-tracking method in order to enforce the interface conditions. Delayed grid points are necessary to allow for the movement of the interface across these grid points. The resulting system of equations contains mainly linear equations. In order to take advantage of this almost-linear system, a predictor–corrector iteration is developed in order to separate the linear and nonlinear equations, avoiding the need to solve a large nonlinear system at each time-step.

With the numerical method in place, the moving boundary problems are solved on an infinite domain with Heaviside initial data. As mentioned previously, the boundary conditions at the interface are formulated from a consideration of a vitrification front that is postulated to occur during the anhydrobiosis of certain organisms by Aksan *et al.* (2006). With vitrification, continuity of water concentration and flux are expected at the glass formation front. Additionally, a fixed glass transition concentration is expected in an isothermal environment. With these conditions, two problems will be considered. The first problem has an interface bounded by an anomalous region to the one side and a classical diffusion region to the other. For notational convenience, this problem is denoted as the anomalous–classical problem. Instead of moving monotonically, the interface is found to reverse direction after a given time, which depends on material parameters. The anomalous–classical problem is then generalized to involve subdiffusion on both sides of the interface. This new problem is denoted as the anomalous–anomalous problem. This interface here can also reverse direction after a given time, but only if the orders of subdiffusion are distinct. The initial direction also depends on physical parameters.

## 2. Formulation

Consider a one-dimensional moving boundary problem on an infinite domain *Ω*, on which a diffusion process is occurring with concentration *c*(*x*,*t*) (0≤*c*(*x*,*t*) ≤ 1). Suppose there is a critical value *c** so that the process is classically diffusive when *c*>*c**, and subdiffusive when *c*<*c**. Define *h*(*t*) as the location of the interface between the classical diffusion and subdiffusion regions. Denote the subdiffusion region as *Ω*_{1} and the classical diffusion region as *Ω*_{2}. In order to avoid unnecessary complication in future notation, it is assumed that the interface will change direction no more than one time. Consider the process to have a constant diffusivity coefficient in *Ω*_{1}, as well as a constant diffusivity coefficient in *Ω*_{2}. As it is a classical diffusion region, *Ω*_{2} is governed by the classical diffusion equation:
2.1
Likewise, *Ω*_{1} is a subdiffusion region and thus governed by the fractional diffusion equation
2.2
where 0<*α*<1 and is the time that the location *x* has entered the subdiffusion region *Ω*_{1} (which may be the initial time *t*=0). Note that *s*(*x*) could depend on time if the interface reverses direction so that *x* enters *Ω*_{1} twice. However, the dependence on time is suppressed as it will not impact the results here. The operator is the Caputo fractional derivative operator of the form
2.3
where *α* is the order of the derivative and is the gamma function. The initial condition will be *c*(*x*,0)=*H*(*x*). Thus, *Ω*_{2} is initially the positive *x*-axis with *c*=1, while *Ω*_{1} is the negative *x*-axis with *c*=0.

The form of the fractional diffusion equation given in (2.2) is novel in its incorporation of the moving interface through *b*=*s*(*x*). Although (2.2) with *b*=*s*(*x*) is reasonable and could be studied for its novelty alone, more fundamental justification is included in the appendix A. The continuous-time random walk approach used in the appendix A has also been used for other physical problems to justify different forms of the fractional diffusion equation (Yuste *et al.* 2010).

Two boundary conditions need to be given at the interface *x*=*h*(*t*). The first of which is that *c*(*x*,*t*) is continuous and equal to the critical value at the interface:
2.4
Note that other boundary conditions are possible here, including the standard instantaneous phase-change problems discussed by Wilson (1978). In those problems, *c*(*x*,*t*) is a concentration that differs at the interface due to a non-unity segregation coefficient. Here, the underlying material is considered to be in the same phase throughout *Ω*. The second boundary condition at the interface is the balance of flux across the interface:
2.5
*J*_{1}(*x*,*t*) is the subdiffusion flux in *Ω*_{1} and *J*_{2}(*x*,*t*) is the classical diffusion flux in *Ω*_{2}. The classical diffusion equation (2.1) can be derived using Fick's first law and balancing the flux locally (∂*c*/∂*t*)+(∂/∂*x*)*J*=0. However, such an expression for the flux is not known *a priori* for subdiffusion, as subdiffusion requires only that 〈*x*^{2}(*t*)〉∝*t*^{α} (0<*α*<1). Thus, the fractional diffusion equation (2.2) must be manipulated into the form (∂*c*/∂*t*)+(∂/∂*x*)*J* to obtain the subdiffusion flux *J*_{1}(*x*,*t*).

### (a) Extracting the subdiffusion flux

First, the fractional diffusion equation (2.2) is converted into the Riemann–Liouville fractional diffusion equation
2.6
where the operator is the Riemann–Liouville fractional derivative operator of the form
The conversion assumes sufficient smoothness of *c*(*x*,*t*) and is accomplished by applying to both sides of (2.2). Note that (2.6) emerges from the CTRW derivation by Metzler & Klafter (2000) for homogeneous domains, where asymptotic expansions are used in Fourier–Laplace space instead of real-space.

When the interface is to the left (*h*(*t*)≤0) with the subdiffusion region receding, it is clear from figure 1 that *s*(*x*)=0, ∀*x*∈*Ω*_{1}.

Thus, the spatial derivative can be interchanged with the fractional derivative in (2.6), and the equation takes on the form . Now, the flux can be written as and the interface condition (2.5) becomes
2.7
Note that the history here is taken at *x*=*h*(*t*) only.

When the interface is to the right (*h*(*t*)>*h*(0)), figure 2 shows that *s*(*x*) ≠ 0, ∀*x*∈*Ω*_{1}.

Thus, the spatial derivative cannot simply be interchanged with the fractional derivative in (2.6). In order to extract the flux, several steps are needed. First, Leibniz's rule is applied to compute
This result is now substituted into the right-hand side of (2.6) to find
For the second integral on the right-hand side, introduce *v*=*s*(*u*). The temporal derivative is now applied within the square bracket to get
After rearranging to the form (∂*c*/∂*t*)+(∂/∂*x*)(*J*_{1}(*x*,*t*))=0, the flux *J*_{1}(*x*,*t*) is easily identified and the interface condition (2.5) becomes
2.8
Note that (2.8) includes history taken along the moving interface as well as history at *x*=*h*(*t*), as in (2.7).

Finally, when the interface is to the left (*h*(*t*)<*h*(0)) and the subdiffusion region is advancing, figure 3 shows again that *s*(*x*)≠0, ∀*x*∈*Ω*_{1}. Note also that *s*(*x*) has a jump discontinuity at *x*=*h**.

Similar steps are used to extract the flux as the previous two cases. However, in order to keep the flux continuous at *x*=*h**, an extra term is added during the application of Leibniz's rule:
Proceeding as before, the Riemann–Liouville fractional diffusion equation (2.6) becomes
Note that the additional integral added does keep *J*_{1}(*x*,*t*) continuous at *x*=*h** although *s*(*x*) has a jump discontinuity there. Thus, the interface condition (2.5) is
2.9
Here, the history is taken partially along the moving interface, partially at *x*=*h*(*t*), and partially at *x*=*h**.

## 3. Numerical method

A numerical method is developed to solve the anomalous–classical problem defined by (2.1), (2.2), (2.4) and (2.5). The anomalous–classical problem will be solved on a finite computational domain *Ω*=[−*L*,*R*] with Heaviside initial data so that *h*(0)=0, *Ω*_{1}=(−*L*,0), and *Ω*_{2}=(0,*R*). The endpoints *x*=−*L* and *x*=*R* are selected large enough so that the solution is unaffected by the size of the domain. Define Δ*x*=(*R*+*L*)/(*M*−1) and the spatial grid points {*x*_{0}=(−*L*),*x*_{1}=(−*L*+Δ*x*),…,*x*_{m}=(−*L*+*m*Δ*x*),…,*x*_{M}=(*R*)} as shown in figure 4. The temporal domain [0,*T*] is discretized similarly with {*t*_{0}=(0),*t*_{1}=(Δ*t*),…,*t*_{n}=(*n*Δ*t*),…,*t*_{N}=(*T*)}. Furthermore, define *I*_{n} to be the index of the grid point just to the left of the interface at *t*=*t*_{n}, so that *x*_{In}<*h*(*t*_{n})<*x*_{In+1}. If *h*(*t*) falls on a grid point *x*_{m}, the time-step is adjusted so that *h*(*t*) lands to the side of *x*_{m}.

### (a) Discretizing the diffusion equations

With *x*_{m}∈*Ω*_{2} for *m*=*I*_{n}+3,…,*M*−1, the classical diffusion equation (2.1) is discretized across these points using the standard backward-time centred-space stencil (denote and ):
With *x*_{m}∈*Ω*_{1} for *m*=1,…,*I*_{n}−2, the fractional diffusion equation (2.2) is discretized across these points. Apply Riemann–Liouville fractional integration to both sides of the fractional diffusion equation (2.2) to obtain the corresponding Volterra integral equation:
Note the relabelling of *T* to *t* in the second line. The general integral is discretized in a manner suggested by Diethelm *et al.* (2002). At location *x*=*x*_{m}, denote *s*_{m} as the index where *t*_{sm−1}≤*s*(*x*_{m})<*t*_{sm}. Then, the discretization is
3.1
where and . A linear approximation can be used to approximate *F*(*x*_{m},*v*) for added accuracy, if needed. Thus, the discretization for (2.2) is

### (b) Including the interface

The remaining governing equations at the grid points *x*_{In−1}, *x*_{In}, *x*_{In+1} and *x*_{In+2} involve the location of the interface. The interface may move through the grid points *x*_{In} or *x*_{In+1} during this iteration, with the time-step restricted so that *h*(*t*) does not go past *x*_{In−1} or *x*_{In+2}. Thus, *x*_{In} and *x*_{In+1} are ignored until *h*(*t*_{n+1}) is known. Thus, the equation at *x*_{In−1} will be the same as *x*_{In−2} except that approximates (∂^{2}*c*/∂*x*^{2})(*x*_{In−1},*t*_{n+1}) using the points *x*_{In−2}, *x*_{In−1} and *h*(*t*_{n+1}). Similarly, the equation at *x*_{In+2} is the same as *x*_{In+3} except that approximates (∂^{2}*c*/∂*x*^{2})(*x*_{In+2},*t*_{n+1}) using the points *h*(*t*_{n+1}), *x*_{In+2} and *x*_{In+3}. The flux balance (2.5) to compute the *h*(*t*_{n+1}) is nonlinear when discretized. Thus, a predictor–corrector iteration is used to separate the linear equations for *c*(*x*,*t*) away from the interface, from the nonlinear equations near the interface and the flux balance for *h*(*t*_{n+1}). Using a guess for *h*(*t*_{n+1}), the predicted concentration values are computed using a linear banded solver at {*x*_{0},…,*x*_{In−1},*x*_{In+2},…,*x*_{M}}. A corrected guess for *h*(*t*_{n+1}) is then computed using Newton's method or a bisection/secant method. Once the iteration for *h*(*t*_{n+1}) converges, the concentration values at *x*_{In} and *x*_{In+1} are computed using the discretization scheme of the region (*Ω*_{1} or *Ω*_{2}) they end up in. In the event that one of these grid points has switched regions, *s*(*x*) is determined using an interpolant between *s*(*h*(*t*_{n}))=*t*_{n} and *s*(*h*(*t*_{n+1}))=*t*_{n+1}. If the new interface position is not within (*x*_{I−1},*x*_{I+2}), a smaller time-step is taken until it is within said region.

In order to discretize the flux balance (2.7), a first-order differencing scheme is used for the Riemann–Liouville derivative:
Now, the same discretization steps from (3.1) are used to obtain
3.2
where *D*_{1L}(*c*^{j}) approximates (∂*c*/∂*x*)(*h*(*t*_{n+1}),*t*_{j}) using *x*_{In−2}, *x*_{In−1} and *h*(*t*_{n+1}). Because the location *h*(*t*_{n+1}) will be in-between grid points for *j*<*n*+1, an interpolated value is used for *c*(*h*(*t*_{n+1}),*t*_{j}). Similar to , . Finally, *D*_{1R}(*c*^{j}) approximates (∂*c*/∂*x*)(*h*(*t*_{j}),*t*_{j}) using *h*(*t*_{j}), *x*_{In+2}, *x*_{In+3}.

In order to discretize the flux balance (2.8), a first-order differencing scheme is again used followed by the discretization steps from (3.1). When the interface is moving to the right, so that *s*(*h*(*t*))=*t*, *D*_{1L}(*c*^{j}) is replaced with *D*_{1Lh}(*c*^{j}) that approximates (∂*c*/∂*x*)(*h*(*t*_{j}),*t*_{j}) using *x*_{Ij−2}, *x*_{Ij−1} and *h*(*t*_{j}). When the interface is moving left, so that *s*(*h*(*t*))<*t*, the scheme requires *m* such that *t*_{m}<*s*(*h*(*t*_{n+1}))<*t*_{m+1}. The scheme then uses *D*_{1L} for *j*>*m*, and *D*_{1Lh} for *j*<*m*. For *j*=*m*, *D*_{1Lh*}(*c*) uses *c*(*h*(*t*_{n+1}),*s*(*h*(*t*_{n+1}))) and interpolated values for *c*(*x*_{In−2},*s*(*h*(*t*_{n+1})) and *c*(*x*_{In−1},*s*(*h*(*t*_{n+1})).

In order to discretize the flux balance (2.9), the same steps are used to discretize as for (2.7) and (2.8). Here, the scheme requires *m* such that *t*_{m}=*t**, and uses *D*_{1Lh} for *j*≥*m*, and *D*_{1L*} for *j*<*m*. *D*_{1L*} approximates (∂*c*/∂*x*)(*h*(*t*_{m}),*t*_{j}) using *x*_{Im−2}, *x*_{Im−1} and *h*(*t*_{m}). The location *h*(*t*_{m}) will be in-between grid points for *j*<*m*, and thus interpolated values are used for *c*(*h*(*t*_{m}),*t*_{j}).

An approximation to the Heaviside function is used for the initial condition along with *c*(−*L*,*t*)=0 and *c*(*R*,*t*)=1 as the boundary conditions. The approximation to the Heaviside function is defined piecewise as
3.3
where *ϵ* is an approximation parameter. This particular initial condition is chosen because it satisfies the flux balance (2.5) uniformly as .

### (c) Handling the subdiffusion history

All previous *c*(*x*,*t*) values in *Ω*_{1} must be stored in order to compute the new *c*(*x*,*t*) values using (3.1). While the computational storage is not an issue for the examples presented here, the corresponding computational time is, as it grows like *O*(*N*^{2}). For this reason, an adaptive time-stepping routine is implemented to decrease *N*. The routine compares the interface locations computed at *t*+Δ*t* and *t*+0.5Δ*t*. If their difference is above a certain threshold (*ϵ*_{max}), the time-step is halved. If their difference is below a certain threshold (*ϵ*_{min}), the time-step is doubled. This is particularly effective here because the time-step restriction depends on the interface speed, which varies over time. Additionally, the *Ω*_{1} region is split up so that the history can be computed in parallel across multiple cores. Although unnecessary for the computations in this paper, higher-order temporal schemes for (3.1) would help decrease *N* further.

### (d) Convergence analysis

In order to examine the global order of the scheme, the values *α*=0.8, , *c**=0.8, *L*=*R*=4 and *ϵ*=0.6 are used. To first study the spatial convergence, results are obtained for *T*=0.01 with Δ*t*=1.0×10^{−6} and *M*∈{51,101,201,401,801}. These results are compared with the solution with the same parameters except *M*=3201:
Note that the presence of the free interface introduces fluctuations in the values for the *L*_{2}(2Δ*x*):*L*_{2}(Δ*x*) and ratios. However, the values still indicate that the method is at least first-order in space. The convergence is not second order because the first-order stencils used near the interface, and the free interface itself, degrade the overall order.

To study the temporal convergence, results are obtained for *T*=0.0001 with *M*=3201 and Δ*t*∈{1×10^{−6},5×10^{−7},2.5×10^{−7},1.25×10^{−7},6.25×10^{−8}. These results are compared with the solution with the same parameters except Δ*t*=1.5625×10^{−8}:
The values for the *L*_{2} and ratios show first-order convergence in time.

The final convergence test is for *T*=50 with *α*=0.8, , *c**=0.8, *L*=*R*=25 and *ϵ*=0.6. For this larger *T*, the adaptive time-stepping routine is essential to manage the *O*(*N*^{2}) computational time. Instead of using a fixed Δ*t*, the time-steps *t*_{1},…,*t*_{M} are generated from the adaptive routine with Δ*t*_{initial}=1×10^{−6}, *M*=16001, *ϵ*_{max}=1×10^{−4} and *ϵ*_{min}=1×10^{−5}. These time-steps are then used to compute results for *M*∈{251,501,1001,2001,4001}, which are compared with the *M*=16001 result:
Thus, even at larger *T*, the convergence is still at least first-order in space.

### (e) Generalizing to the anomalous–anomalous problem

As noted earlier, the more general anomalous–anomalous problem will also be considered. The numerical procedure to solve the anomalous–anomalous problem can be constructed by paralleling the numerical method for the anomalous–classical problem presented earlier. When the interface is moving to the right, the new subdiffusion region *Ω*_{2} will be receding and will use the existing receding stencil. When the interface is moving to the left, *Ω*_{2} will be advancing and will use the existing advancing stencil. The remaining details are straightforward and will not be discussed.

## 4. Numerical results

The anomalous–classical problem defined by (2.1), (2.2), (2.4) and (2.5) is solved using the numerical method presented in §3. Results for the interface position *h*(*t*) with *α*=0.8 in *Ω*_{1} are presented in figure 5.

Figure 5*a* shows how the interaction between these two regions causes the interface in the anomalous–classical problem to evolve on an overall time scale that is slower than that of pure classical diffusion (*α*=1.0). In order to determine that overall time scale, curves of the form *k*_{1}*t*^{k2} are fitted to the interface position. Figure 5*b* shows that the curve for *α*=1.0 is of the expected form *h*(*t*)∼*t*^{1/2}, while the curve for *α*=0.8 implies *h*(*t*)∼*t*^{ν} where *ν*<*α*/2. This suggests that the overall time scale is not only slower than that of pure classical diffusion, but even slower than that of pure subdiffusion. However, when the solution is computed for longer time, it is clear that *h*(*t*)≠*k*_{1}*t*^{k2} for *α*=0.8.

In figures 6 and 7, *h*(*t*) is plotted versus *t* for various values of *α* and *c**. Note that *h*(*t*) increases at first, but then decreases as time moves forward. The time at which the interface reverses direction will be denoted as *t**, with *h** defined as *h**=*h*(*t**).

Figure 6*a* shows that for *c**=0.65, the more anomalous the subdiffusion (smaller *α*), the quicker and shorter the direction change (smaller *t** and *h**). With pure classical diffusion (*α*=1) and *c**>0.5, Gruber *et al.* (2011) show that the interface moves monotonically to the right (). Thus as the subdiffusion tends towards classical diffusion (), the interface returns to being monotone (*t** and *h** tend to ) as expected. Additionally, figure 6*b* shows that for *c**=0.4, the more anomalous the subdiffusion, the longer and further the direction change (larger *t** and *h**). With pure classical diffusion and *c**<0.5, Gruber *et al*. also show that the interface would move monotonically to the left (*t**=0=*h**). Again, as the subdiffusion tends towards classical diffusion (), the interface returns to being monotone (*t** and *h** tend towards 0) as expected. With *α* fixed, figure 7 shows that decreasing *c** causes *t** and *h** to decrease. The dependence of *t** and *h** on both *α* and *c** will be further explored in §4*a*.

The direction reversal is a result of the coupling of the subdiffusion and classical diffusion regions. In order for the *Ω*_{1} region to continue to advance, the diffusion within *Ω*_{1} must be faster than the diffusion into the region from *Ω*_{2}. If *α* is set to 1 (*Ω*_{1} follows classical diffusion), then the interface motion returns to being monotone, as expected by the self-similarity with . This is because the diffusion within *Ω*_{1} is faster than diffusion into the region from *Ω*_{2} for *c**>0.5. However, when the two regions are on different time scales (*α*<1), the diffusion will only be faster for small time. Whereas the *Ω*_{1} region may start off evolving faster than the *Ω*_{2} region (*t*^{α/2}>*t*^{1/2} for small *t*), it will eventually evolve slower (*t*^{α/2}<*t*^{1/2} for large *t*). Thus, the diffusion within the *Ω*_{1} region will eventually be too slow to keep up with the diffusion into the region from *Ω*_{2} and the *Ω*_{1} region will begin to recede.

Evidence of this can be seen in the *c*(*x*,*t*) profiles in figure 8. Figure 8*a*,*b* shows the development of a ‘kink’ in the *c*(*x*,*t*) profiles over time where the subdiffusion and classical diffusion regions meet. This forms as the diffusion within the *Ω*_{1} region falls behind that of the diffusion into the region from *Ω*_{2}. As expected, the *c*(*x*,*t*) profile with the greater disparity in time scales (*α*=0.7) shows greater kink formation than the *c*(*x*,*t*) profile with less disparity (*α*=0.9). Note that similar kink formation will occur with classical diffusion on both sides of the interface (*α*=1) when . There, the kink is also caused because the diffusion in *Ω*_{1} is faster/slower than the diffusion in *Ω*_{2}. However, because both regions involve classical diffusion, and thus evolve on the same time scale, the interface remains monotone , as shown by Gruber *et al.* (2011).

### (a) Quantifying *t**(*α*,*c**) and *h**(*α*,*c**)

In order to explore the dependence of *t** on *α* and *c**, numerical results were obtained for a range of *α* and *c** values. The results suggest the following relation:
4.1
The base of the power law in (4.1) is determined by considering the monotone interface direction when *c**=0 and when *c**=1. The exponent is found from the slope of a linear fit for *y*=*f*(*x*), where and . The values for *β*(*α*) in figure 9*a* are obtained by averaging across various *c** and *α* values. The results are fitted with a power law curve of the form *k*_{1}*α*^{k2}.

Additionally, the results suggest that *h**=*θ*(*α*)(*c**/(1−*c**))^{f[1/(1−α)]}, where *f*(*y*)=1.018*y*−0.7317 after a linear fit is found. Again, the values for *θ*(*α*) in figure 9*b* are then found by averaging across various *c** and *α* values. This time, the results are fitted with a curve of the form *k*_{1}(1−*α*)/*α*+*k*_{2}.

Although the fitted curves in figure 9 will suffice for the purposes here, better fits may be possible for with data for *α*<0.5. However, decreasing *α* causes *h*′(*t*) to increase for small *t*. Thus, the problem becomes very stiff and computationally expensive for *α*<0.5. Regardless, the dynamics for *α*<0.5 are expected to follow (4.1). The numerically computed *t** and *h** values are shown to match very closely with their respective fitted curves in figure 10.

Numerical results were also obtained for larger values of *L* and *R* (*L*=*R* = 400 instead of *L*=*R*=100) and compared with those presented earlier. The difference between the results was negligible, indicating that the values used for *L* and *R* are sufficiently large that the boundary conditions do not affect the solution. Additionally, results were obtained for the smoother initial condition (*x*/*ϵ*)/*π*+0.5. These results also show the subdiffusion region initially advancing until time and then receding. Furthermore, these values have the same form as (4.1), but with different values for *β*(*α*).

### (b) Results for the anomalous–anomalous problem

The *Ω*_{2} region is now considered to be a subdiffusion region with *α** as the order of the fractional derivative in (2.2). When the order of the derivative in *Ω*_{2} is equal to that in *Ω*_{1} (*α**=*α*), the anomalous–anomalous problem can be shown to follow the self-similarity *ξ*=*x*/*t*^{α/2}. Thus, the corresponding interface position is given by *h*(*t*)=*λt*^{α/2} for some constant *λ*. This can be seen in the numerical results obtained for *α**=*α* in figure 11.

The interface location data in figure 11*a* is fitted with curves of the form *k*_{1}*t*^{k2}. The self-similarity of the anomalous–anomalous problem with *α**=*α* is confirmed by *k*_{2}≈*α*/2 for both *α*=0.7 and *α*=1.0. Furthermore, figure 11*b* does not show a kink because the two subdiffusion regions are evolving at the same rate and on the same time scale.

In order to evaluate the behaviour when *α* is distinct from *α**, both *α**>*α* and *α**<*α* are studied. As expected, the disparity in time scales of the subdiffusion in *Ω*_{1} and *Ω*_{2} causes the interface to reverse direction after a specific time *t**. Similar to the anomalous–classical problem, the more subdiffusive (lesser of *α* and *α**) region will initially evolve faster than the less subdiffusive (greater of *α* and *α**) region. Thus, initial direction of the interface is always towards the less subdiffusive region. In figure 12, the *t** values are plotted against the *α** values.

When *α**>*α*, the interface initially moves towards *Ω*_{2}. Thus, corresponds to a monotone interface towards *Ω*_{2}, and *t**=0 corresponds to a monotone interface towards *Ω*_{1}. As , figure 12 shows for *c**=0.6 and for *c**=0.4. When *α**<*α*, the interface initially moves towards *Ω*_{1}. Thus, corresponds to a monotone interface towards *Ω*_{1}, and *t**=0 corresponds to a monotone interface towards *Ω*_{2}. As , figure 12 shows for *c**=0.6 and for *c**=0.4. Note in both cases, the interface returns to the monotone behaviour found by Gruber *et al.* (2011) for *α**=*α* (towards *Ω*_{2} for *c**>0.5 and towards *Ω*_{1} for *c**<0.5).

## 5. Conclusions

In order to study the effect of subdiffusion on moving interfaces, two moving boundary problems are solved. The first moving boundary problem, the anomalous–classical problem, is formulated as a diffusion process that is measured by *c*(*x*,*t*), governed by (2.1) and (2.2), and subject to (2.4) and (2.5) on the interface. Whereas the interface location *h*(*t*) in moving boundary problems involving only classical diffusion follows the power law *h*(*t*)=*λt*^{1/2}, the anomalous–classical problem reveals that the presence of subdiffusion causes the interface to reverse direction. Furthermore, the dependence of the reversal time *t** on the order 1−*α* and critical value *c** is determined to be (4.1). The direction change distance *h** has a similar form.

The second moving boundary problem, the anomalous–anomalous problem, generalizes the anomalous–classical problem to include subdiffusion regions on both sides of the interface. When the two subdiffusion regions are of the same order, the interface location *h*(*t*) follows the power law *h*(*t*)=*λt*^{α/2}. However, the anomalous–anomalous problem confirms that the interface will reverse direction when the orders of the two subdiffusion regions are distinct. Additionally, the interface will initially move towards the less subdiffusive region (greater of *α* and *α**).

The driving force behind the direction reversal is found to be the difference in time scales of the two regions *Ω*_{1} and *Ω*_{2}. This disparity is seen by the formation of a kink in the *c*(*x*,*t*) profile. If *Ω*_{1} is more subdiffusive than *Ω*_{2}, the *Ω*_{1} region evolves faster initially. The diffusion within *Ω*_{1} causes the region to advance initially, driving the interface towards *Ω*_{2}. However, the *Ω*_{1} region will eventually evolve slower, at which point the diffusion from *Ω*_{2} into *Ω*_{1} will overwhelm the diffusion within *Ω*_{1} and cause the interface to reverse direction. If *Ω*_{1} is less subdiffusive than *Ω*_{2}, the opposite interface behaviour will occur.

These results may suggest some interesting dynamics in the desiccation experiments by Aksan *et al.* (2006) and Chakraborty *et al.* (2011). Assuming that trehalose glass forms in these experiments, and that the diffusion of water follows subdiffusion within that trehalose glass, the results here suggest a non-monotone vitrification front. In the experiments by Aksan *et al.* (2006), they attempt to desiccate microchannels filled with trehalose-water solution. They hypothesize that the non-uniform drying observed is caused by a skin of trehalose glass. If the glass is present and the channel is left to equilibrate without the drying mechanism, the anomalous–classical results might suggest that the glass region would eventually recede and disappear. Thus, an intermittent drying mechanism may be more efficient. In the experiments by Chakraborty *et al.* (2011), they attempt to preserve cells after the samples have been vitrified by spinning. While the cells remain viable for very short storage times, the cells do not survive for longer times. Again, the anomalous–classical results might suggest that the presence of any unvitrified region would cause the glass region to recede. Eventually, the sample would fully reconstitute and leave the exposed cells in an environment where they cannot survive. Thus, the samples may need to be fully vitrified in order to store them.

## Acknowledgements

We acknowledge partial financial support from NIH grant 1R01GM086886 and NSF RTG grant DMS-0636574.

## Appendix A. Continuous-time random walk derivation

Consider a random walker whose probability of being in location *x* at time *t* is given by *p*(*x*,*t*). Following the notation of Metzler & Klafter (2000), let *λ*(*x*) be the jump-length probability density function (PDF) for the walker, let *w*_{s}(*t*) be the wait-time PDF while the walker is in the subdiffusion region *Ω*_{s} and let *w*_{c}(*t*) be the wait-time PDF while the walker is in the classical diffusion region *Ω*_{c}. Define as the probability the walker remains at its current location, in *Ω*_{s} or *Ω*_{c}, for a period of time *t*. Finally, let define *η*(*x*,*v*) as the probability the walker arrives at location *x* at time *v*. Although the following derivation can be applied to both leftward and rightward moving interfaces, the focus here is on the rightward moving case (an anomalous region moving into a classical region). The following expressions are now defined:
A1
A2
and
Note that *G*(*x*,*t*) uses a one-sided delta function to handle when the walker enters a region without jumping. Whereas Metzler & Klafter (2000) define *w*_{s}(*t*) in terms of long-time behaviour only, all PDFs are explicitly defined here:
A3

The analysis here is not limited by the use of specific PDFs. The crucial assumptions of the derivation are that the power law tail in *w*_{s}(*t*) matches the Metzler & Klafter long-time behaviour, and that the functions are PDFs. The standard approach to deriving a fractional diffusion equation, as seen in Metzler & Klafter (2000), is to transform (A1)–(A3) into Fourier and Laplace space. The resulting equations are exchanged for their long-time approximations and then solved. The approach here uses a similar limit, while remaining in physical space, by letting *τ* and *σ* tend to zero.

The focus is narrowed to *x*∈(*h*(0),*h*(*t*)), as this is the most interesting region. The *η*(*x*,*v*) in (A1) can be replaced by its definition (A2). The contribution from is exponentially small by Laplace's method. Integration-by-parts techniques show that *ψ*_{s}(*b*−*a*)=*O*(*τ*^{α}) for *b*≠*a*. This, along with careful application of Laplace's method, gives that the contribution from is *O*(*σ*^{γ}*τ*^{α}) where *γ* is a parameter chosen in Laplace's method (0<*γ*<1). The remaining expression for *p*(*x*,*t*) is

By changing the order of integration, noting *ψ*_{s}′(*t*)=−*w*_{s}(*t*), and using integration-by-parts on , the expression for *p*(*x*,*t*) becomes
Again, a careful application of Laplace's method shows the double integral is *O*(*σ*^{γ}*τ*^{α}). Using (A1), the integral equation for *p*(*x*,*t*) is now complete:
A4

Denote and use integration-by-parts to get
In the first integral, substitute *F*(*v*)=*F*(*t*)−(*t*−*v*)*F*′(*μ*(*v*)) where *v*≤*μ*(*v*)≤*t*. The contribution from the *F*′(*μ*(*v*)) term is *O*(*τ*^{2−α}). The result is
The remaining integrals are expanded using Laplace's method, proceeded by changing the order of integration if necessary:
After using that *ψ*_{s}(*t*−*s*(*x*))∼*τ*^{α}/[*Γ*(1−*α*)(*t*−*s*(*x*))^{α}] from integration-by-parts, the parameters *σ* and *τ* are both taken to zero so that . The leading-order result is the fractional diffusion equation (2.2):
A5

- Received March 19, 2012.
- Accepted May 21, 2012.

- This journal is © 2012 The Royal Society