A direct sampling method for time-fractional diffusion equation

This paper introduces a direct sampling method tailored for identifying the location of the source term within a time-fractional diffusion equation (TFDE). The key aspect of our approach involves the utilization of a versatile family of index functions, which can be chosen according to the specific characteristics of the source term. Recognizing the key role of the TFDE’s fundamental solution within the index function, we further enhance our method by deriving its asymptotic expansions. This advancement not only enhances the accuracy, but also significantly improves the computational efficiency of our method. To validate the effectiveness and robustness of the proposed sampling method, we conduct a series of comprehensive numerical experiments.


Introduction
Fractional differential equations are equations that involve fractional derivatives or integrals.They are used to describe a wide range of phenomena with non-local mechanisms (see [35]    * Author to whom any correspondence should be addressed. Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
for some examples).In particular, a time-or space-fractional diffusion equation is an effective model for anomalous diffusion processes, including diffusion in fractal media [34], diffusion or material transport in 1D systems [21,26], diffusion of cosmic rays [20,49], and release of drug [48].Additionally, it models dielectric relaxation [33].Metzler and Klafter [29] derived the time-and space-fractional diffusion equations by studying the stochastic movement of particles.
We will focus on the time-fractional diffusion equation (TFDE), which reads as follows: Here, Ω ∈ R n is a bounded region with smooth boundary, γ ∈ (0, 2] is the order of the fractional derivative, κ is the conductivity constant, ∂ n is the normal derivative and g ∈ S ′ (Ω × [0, +∞)) is the source term.The fractional derivative is taken in Caputo's sense, i.e.
For 0 < γ < 1, the TFDE characterizes diffusion processes that are slower than usual, hence it is called a time-fractional sub-diffusion equation [29].Sub-diffusion usually occurs when the medium is porous or has fractal structure, such as soil, gas reservoirs and reinforced concrete [4,12,43,45].Experiments have been conducted to show the validity of sub-diffusion equations.When 1 < γ < 2, the TFDE becomes a time-fractional diffusion-wave equation [27].This equation is primarily of mathematical interest as it shares some similarities with the sub-diffusion equation.The TFDE also finds its applications in signal processing [23] and image smoothing [9].
In this article, we aim to find the support of the source term g from measurements u| (∂Ω)×[0,T] on the boundary.This can be applied, for example, to locate pollutant sources when it is not possible to measure data in the interior of a region.There are already numerous results for the forward problem, both analytical [15,25,39] and numerical [1,10,11,31,37,44].For the inverse problem, Sakamoto and Yamamoto [39] proved the stability of reconstructing the temporal component of the source term.Liu and Zhang [24] lifted an impractical assumption in their result.Li et al [22] established the uniqueness and stability of simultaneous recovery of the space-dependent diffusion coefficient and the fractional order, while Ruan et al [38] proved the counterpart of recovering the fractional order and the space-dependent source term.Numerically, a commonly used method is to minimize a cost functional.Wei and Zhang [46], and Gong and Wei [14] constructed the cost functional by discretizing the boundary integral equations corresponding to the TFDE and applied Tichonoff regularization to address the ill-posedness of the problem.Ruan et al [38] proposed a multi-term cost functional and designed an alternating minimization algorithm to recover the coefficients.Apart from minimization methods, Liu and Zhang [24] proposed an iterative method to recover the temporal component of the source term.The reciprocal gap method can be used to determine the location of point sources on R 2 along with other information in evolution equations; see [2] for its application in the heat equation and [13] in the TFDE.An exhaustion-type method is also proposed by Yang and Deng [47] to recover the location of point sources.
Another significant method is the sampling method.The sampling method primarily seeks an index function characterizing the location of inclusions, obstacles, inhomogeneities, etc.The sampling method was introduced by Colton and Kirsch [8], where they used far-field observation data to determine the position of a scatterer under the Helmholtz equation.The sampling method consists a variety of sub-methods, including the factorization method [18], the linear sampling method [16,19,32], reverse time migration [5], direct (orthogonality) sampling method [3,6,7,36] and the topological sensitivity method [40,41].In the direct sampling method, or the orthogonality sampling method, the index function is constructed by pairing the boundary observations with some probing function under an inner product.Then an orthogonality argument is applied to obtain the desired property of the index function.It offers advantages over the aforementioned three approaches.In those methods, the measurement is an integration kernel of two spatial variables (or its equivalent, such as a Neumann-to-Dirichlet map), which is difficult to obtain and computationally inefficient.On the other hand, the direct sampling method only requires a single set of observation data on the boundary, and the probing function can be computed offline.The topological sensitivity method is based on the idea of minimizing the cost functional, but it only requires one iteration.This approach rivals the direct sampling method in simplicity and the amount of computation.Inspired by the work of Chow et al [7], we propose a family of index functions that can adapt to the properties of the source term g.Given that the probing function relates to the fundamental solution of the TFDE, we would also like to focus on efficient and accurate evaluation of the fundamental solution.
The rest of this article is organized as follows.In section 2, we will describe the sampling method and discuss the choice and the validity of the index function.In section 3, we will derive the asymptotic expansion for the fundamental solution of the TFDE when the dimension n ⩽ 3.In section 4, we describe the detailed algorithm and present the numerical results.Finally, in section 5, we provide some concluding remarks.

Sampling method for TFDEs
Consider the TFDE defined over the entire space R n : The source term g is represented as m j =1 g j (t)δ xj (x), where x j ∈ Ω. Hereafter, for a function ϕ(y 1 , . . ., y N ) with y j ∈ R dj , we denote the partial derivative in the direction of y j by ∂ j ϕ if d j = 1 and denote the L p norm with respect to y j by ||ϕ|| L p j .Denote the fundamental solution of the equation (2) as solves equation ( 2) with g replaced by δ (x,t) .We provide some elementary properties of the fundamental solution.
Proof.For x, y ∈ R n and t, s ⩾ 0, define It is straightforward to check the above relation also holds for γ = 1.Moreover, for any γ ∈ (0, 2) we have ∆H (x,t) (y, s) = ∆G (0,0) (y − x, s − t).Hence, Note that H (x,t) also satisfies zero initial condition.Therefore, we obtain Now the statements in the proposition follow directly from above.
Let f be the boundary data of equation (2), i.e. f = u| (∂Ω)×[0,T] .Using a modification of the index function introduced in [7], we define where α, β ∈ N. Eventually, we find the location of the sources by searching the region where the index function I α,β takes relatively large values.Hence, the index function is constructed so that it approximately reaches its maximum when x is one of the source locations x j .Note that our index function includes two parameters α and β.This gives us the flexibility to fit it into the properties of the source term g.With the index function defined, we can now explore its properties under different conditions.By employing Green's formula for the TFDE (see, for instance, [42]), the boundary value is expressed as It follows that From equation ( 4) and the definition of the index function (3), it is natural to define Note that, for fixed y and s, K α,α attains its maximum 2 (0,T) when (x, t) = (y, s).This is crucial in order to let the index function have the desired property.
In addition, we show a decay property for K α,β (x, t; y, s).By proposition 2.1, for x ∈ ∂Ω, y ∈ Ω and 0 ⩽ s ⩽ T, we have It follows that The right hand side of the inequality above ranges from [0, 1] and is decreasing with respect to s.Now we will discuss the choice of α and β for two typical scenarios.
Assuming g j (0) > 0 is not negligible compared to |g ′ j (t)| near t = 0 for some j.When (x, t) = (x j , 0), the first term inside the bracket reaches its maximum, while the second term remains small relative to the first term.This is because g ′ j is small near s = 0 compared with g(0), and when s is far away from 0, estimate (6) guarantees that K α,α (x, t; x j , s) has low magnitude compared with the maximum value 2 (0,T) of K α,α (x, t; x j , 0).Hence, we can choose I α,α+1 (x, 0) as the index function.
While we could extend this discussion further, doing so would necessitate the use of index functions that incorporate higher-order derivatives.These are known to introduce numerical instability.Instead, we opt to refine our approach by transitioning from point sources to sources that resemble points.Specifically, we consider g(x, t) = m j =1 g j (t)h j (x), where each h j is nonnegative and the supports ω j = supp(h j ) are disjoint, compact subsets of Ω.Given the conditions outlined in Case I, a similar calculation yields: When ω j is small, for any x j belonging to ω j , serves as a good approximation to I α,α+1 (x, t).Thus, we expect that I α,α+1 (x, t) will have relatively large values when x ∈ ω j and t = 0.The analog of Case II remains valid.Finally we consider the TFDE in the bounded domain (1).Let G Ω be the corresponding fundamental solution.If we define then with the same definition of the index function I α,β , the relationship between I α,β and K α,β remains consistent.As G Ω closely approximates G when the domain Ω is large, we expect that the index function to remain valid in bounded domains.
To conclude this section, we aim to present the stability of the index function with respect to the source term.Furthermore, we validate the index function when the source term is nonregular in time, i.e. if g j is sufficiently close to a regular function g j in some L p norm, then the index functions discussed previously remain applicable.
, where x j ∈ Ω or ω j := supp h j ⊂⊂ Ω, and the index functions I α,β and I α,β are defined as in equation ( 3) associated with g j and g j , respectively, with f being the boundary data of either equation (1) or equation ( 2) with source term g, then Proof.We still begin with the case of whole space with point sources.Let f be the boundary data corresponding to the functions g j .Since where 1 p + 1 q = 1.By definition of K α,β (5), we have , by the triangle inequality, we can repeat the proof above with all the norms taken into the integral ´ωj .It follows that For the bounded domain case, the proof above remains valid, with the only modification being the replacement of the fundamental solution G with G Ω .Remark 2.3.One can also recover the sources using partial boundary measurement.Let µ ⊂ ∂Ω be a subboundary, and define , analogous to equations (3) and (5).By the same argument, K α,α still takes its maximum at (x, t) = (y, s), and the decay property (6) still hold.However, this will suffer from lack of accuracy when determining the sources numerically.

Asymptotic expansion of the fundamental solution
To enhance the computational efficiency of the proposed index function, we next derive an asymptotic expansion for both the fundamental solution and its derivatives.While Kim and Lim [17] have provided the leading term of this expansion using Fox-H functions to represent the fundamental solution, this single-term approximation lacks the precision required for subsequent computations.Directly extending Kim and Lim's method to obtain additional terms presents significant challenges.Consequently, we adopt an alternative strategy.
According to [28], the fundamental solution of the TFDE on R n is given by Here, E γ (z) is the Mittag-Leffler function and F n represents the n-dimensional Fourier transform.This formulation has been verified to hold true for 0 < γ ⩽ 1 and n = 1, but it extends naturally to the scenario when 0 < γ < 2 and n ∈ N. When γ = 2, the TFDE transforms into a wave equation, and the fundamental solution incorporates Dirac delta functions.Thus, our subsequent discussions are confined to 0 < γ < 2. Observing that G (x,t) is merely a translated version of G (0,0) , our focus narrows to the latter.By the scaling property of the Fourier transform, we obtain This leads us to focus on the expansion of the reduced Green function . Henceforth we omit the subscript n on Φ n,γ when no confusion arises.
We focus on cases when n ⩽ 3.In theorem 3.1 we will find the asymptotic expansion of the reduced Green function when the dimension n is 1.It is achieved based on an integral representation of the reduced Green function together with the steepest descent method.Subsequently, we show that the three-dimensional reduced Green function has a strong relationship with the one-dimensional one, which is established via a density argument in proposition 3.4.This provides the expansion when n = 3.For the case n = 2, we relate it to the aforementioned scenarios by augmenting the dimension of the underlying equation, and the expansion is detailed in theorem 3.6.

Case of n = 1
Theorem 3.1.When n = 1 and 0 < γ < 2, the reduced Green function Φ γ has the following asymptotic expansion, with coefficients a > 0 and a k : Moreover, termwise differentiation provides the asymptotic expansion for its derivatives The symmetry of Φ γ confirms it is an even function.Hence, we consider x > 0 for the subsequent analysis.
Given n = 1, the reduced Green function is also represented by , where M is the M-Wright function [28].By the integral representation of the M-Wright function [28], we get where Ha is the Hankel path.Repeated changes of variables yield Step 1: our goal is to modify the contour in the integral (11) to be more suited for the steepest descent method.First, write ν = γ 2−γ ∈ (0, ∞), then we can express Let , and let σ 1 be a path crossing z 0 on which η 1 is real.Write z = re iθ , then Im η 1 (z) = 0 implies So σ 1 is parametrized by (r(θ), θ) with r(θ) determined by (13).On the path σ 1 , We denote the expression on the right hand side to be η 2 (θ).
In the following proof, we use b k to denote coefficients of a series, which may vary in different contexts.For θ ∈ (0, θ 0 ), we have Note that all the series involved here are convergent series.It follows that for w ∈ [w 0 , +∞).
To derive the expansion of η 3 , we will make use of equations ( 14) and ( 15) to get The claim is proven by substituting equation ( 18) into the expansion of η 3 .
Combining the claim with equation ( 17) and applying Watson's lemma, we can derive This matches the desired asymptotic expansion.Substituting ν with γ 2−γ yields an expression identical to equation (10), where a = w 0 > 0 and a k = Γ(k+1/2) Step 3: to derive the asymptotic expansion of Φ ( j) γ , we use equation (17) and deduce , where F i (x, w) is of the form b i x pi (w − w 0 ) qi exp(−x ν+1 w) and b i is a constant depending only on w 0 and ν.Then we can repeat the process of equation (19) to derive expansions for each integral involving F i (x, w).Denote these expansions by by equation (19).On the other hand, according to the form of F i (x, w), we have Thus, This shows that the asymptotic expansion of Φ γ can be calculated by differentiating the expansion of Φ γ termwise.Remark 3.2.Let ν = γ 2−γ as above, then the coefficient a in the expansion is ν ν /(ν + 1) ν+1 , which has been shown in the proof.The first few coefficients a k are given in the table 1 below: These coefficients can be computed using mathematical software, such as Mathematica.Note that when γ = 1, the first term of the expansion is 1   2 √ π exp − 1 4 |x| 2 , which is precisely the reduced Green function of the 1D heat equation.
Consequently, the derivative of the first K terms in the asymptotic expansion of Φ γ comprises K + 1 terms, whose first K terms (sorted in descending order by the exponent |x|) constitute the first K terms of the expansion of Φ γ .The same is true when j is other than 1.

Case of n = 2 and 3
In this section, we show that the reduced Green function in two and three dimensions are closely related to that in one dimension.Specifically, for n = 3, we present the following proposition: ).This means that the Fourier transform can be interpreted in the L 2 sense.Let {u N } N be any sequence in ).In the following proof, we fix a sequence {u N } N tending to Φ 3,γ such that u N ∈ C ∞ c (R 3 ) and is rotation-invariant, i.e. u N (x) = u N (|x|, 0, 0).For any rotation-invariant function u, denote ũ(x) = u(|x|, 0, 0), x ∈ R. We aim to first establish the relation: for any rotation-invariant u ∈ C ∞ c (R 3 ).Subsequently, we apply this to u N and show the validity of taking the limit as N → ∞.
Since u N is invariant under rotation, the same is true for its inverse Fourier transform.Denoting v 1 the first component of a vector v, we have The inner integral can be expressed as: It follows that

Note that
Thus, The last equality holds because F −1 1 ũ is an even function on R.This confirms equation (20).Apply equation (20) to u N , then the left-hand side converges to Φ 3,γ in L 2 (R 3 ) by Plancherel's identity.To show that the right hand side converges as well, for any rotationinvariant u ∈ L 2 (R 3 ), we establish These identities confirm the convergence of the right-hand side, completing the proof.
The asymptotic expansion of the reduced Green function in three dimensions is a direct consequence of theorem 3.1 and proposition 3.4.Now let us turn to the case where n = 2.
, where the subscript n on G and G indicates the dependency on the dimension, then we have and (b) The reduced Green function in two dimensions satisfies Proof.We can treat the two-dimensional problem as a specific case of the three-dimensional problem.More precisely, for any bounded function ), consider the classical solution u of the equation Recall equation (9), which tells us that By the exponential decay of Φ 3,γ (theorem 3.1 and proposition 3.4), passing R to infinity, we get the pointwise limit ũ (x, t) = ˆt 0 ˆR3 G 3 (x, t; y, s) F (y, s) dy ds.
It follows by Fubini's theorem and the arbitrary choice of F that In particular, Since the transformation from G n,(0,0) to G n,α is only concerned with the temporal variable, we can replace G n,(0,0) by G n,α in equation ( 21), and part (a) is proved.
To prove part (b), take t = κ −1/γ in equation ( 23), and we get By proposition 3.4, which presents the relationship between the 1D and 3D Green's functions, we obtain and the proof is complete.
Proposition 3.5 allows us to derive the expansion of Φ 2,γ .
Theorem 3.6.When n = 2 and 0 < γ < 2, the reduced Green function has the following asymptotic expansion: for some coefficients a > 0 and a 2,k .Moreover, termwise differentiation yields the asymptotic expansion of its derivatives, Φ Here the derivative is taken in the direction of |x|, as Φ 2,γ is rotation-invariant.
Proof.The proof utilizes Watson's lemma.
Throughout this proof, let b k represent the coefficients of a series, which may vary in different contexts.According to theorem 3.1 and remark 3.3, when y > 0, we have where R K (y) = o y −2K 2−γ .Now we can calculate the asymptotic expansion of Φ 2,γ by proposition 3.5(b), which implies that where and By a change of variables, we can express I k as Since F(r) := (1 + r) (1+p)(1−γ/2)−1 ((1 + r) 2−γ − 1) −1/2 has an expansion of the form l⩾0 b l r l−1/2 , by Watson's lemma we have which leads to , for any ϵ > 0, there exists which implies that Equations ( 25)- (27) show that Since we can differentiate termwise in the expansion of Φ ′ 1,γ (y) by theorem 3.1, for any fixed r > 0 we have Then by proposition 3.5 and the definition of I k , where With a similar proof as shown in theorem 3.1, one can differentiate termwise to obtain the asymptotic expansion of I ( j) k .Besides, J K = o(I ( j) K ) can be inferred through a similar argument used to estimate J K .The above facts together with equations ( 25), ( 26) and ( 28), imply that one can differentiate termwise to obtain the asymptotic expansion of Φ ( j) 2,γ .
Remark 3.7.Let ν = γ 2−γ , then the coefficient a in the expansion is ν ν /(ν + 1) ν+1 as in the one-dimensional case, and the first few coefficients a 2,k are given in the table 2 below.

Detailed algorithmic approach
Broadly, the algorithm is divided into three steps: • compute the reduced Green function by its definition when |x| is small, or by its asymptotic expansion when |x| is large; • recover the fundamental solution (and its derivatives) from the reduced Green's function (and its derivatives); • compute the index function using equation (3).
The details and subtleties of the process will be described below.We start with the case n = 1.
Step 1: when |x| is large, the reduced Green function and its derivatives are computed via the asymptotic expansion in theorem 3.1.When |x| is small, we use the integral formula (11) instead.By differentiating both sides of the first equality in equation ( 11), we get Let the integrand in equation ( 29) be A(r), then the infinite integral is computed as following: Additionally, we truncate the second integral to the interval [ϵ, 1] for some ϵ ∈ (0, 1).The following proposition gives an error bound of the truncation.
Remark 4.2.(a) The assumptions are practical since we will cover the case |x| 1 by the asymptotic expansion.A valid tuple of parameters, for example, can be (γ 0 , x 0 , B, δ) = (1.5, 10, 10 4 , 0.5).(b) The assumption on B is optimal in the sense that as γ 0 tends to 2 and (2−γ0) for some δ ′ > 0, then by the stationary phase theorem, the integral 1 π γ ´∞ B A(r) dr has at most polynomial decay on B. Proof of proposition 4.1.By a simple observation, Let I(j) denote the right hand side of the inequality above.
First, let us consider I(t) for t ⩽ 0. Due to the assumption on B, we have 2 γ r 2/γ−1 − x 0 ⩾ γ 2δ r 2/γ−1 for r ⩾ B. Hence, When t > 0, we can integrate by parts to obtain Again by 2 γ r 2/γ−1 − x 0 ⩾ γ 2δ r 2/γ−1 for r ⩾ B, we have Therefore, the following estimate holds: where and k 0 is the largest k such that t − 2k/γ > 0. Now the estimate for I(t) when t ⩽ 0 yields The proposition is proved by applying the definition of c 0 .
Step 2: the relationship between the fundamental solution and the reduced Green function is given by equation (23).We only need to take derivatives on both sides of this equality.
Step 3: instead of computing the index function directly, we first make an important observation.By the chain rule, ∂ α 2 G (0,0) (x, t) is of the form where c α,j is independent of t.When γ is small, the exponent of t in the denominator is large, which suggests that the function ∂ α 2 G (0,0) may be concentrated near t = 0.The following figure 2 justifies this assumption.
Such a sharp 'spike' will increase the difficulty of integration when computing the index function, so we do a change of variable ) as in proposition 3.5.After some tedious calculation we have Let f be the boundary data.Since in application f is only given on a grid of ∂Ω × [0, T], we can only perform the change of variable on the denominator of the index function.Hence, the index function is computed by a quadrature to the formula where ∂ β 2 f is computed by some finite difference formulas.Now let us consider n = 2 and 3. When n = 3, the reduced Green function is related to the one-dimensional function by proposition 3.4.After that we simply repeat steps 2 and 3 above, and we omit the details.When n = 2 and |x| is large, the reduced Green function is computed by theorem 3.6, and then we still repeat steps 2 and 3.For small |x|, we can obtain ∂ α 2 G 2,(0,0) and G 2,α from ∂ α 2 G 3,(0,0) and G 3,α by proposition 3.5(a).Then we can directly apply equation (30).There are two remaining issues: (a) the number of terms to use in the asymptotic expansion is not determined; (b) the terms 'small' and 'large' are vague, i.e. the dividing line between two methods of computing the reduced Green function is not determined.(When n = 2, we bypass the computation of the reduced Green function in the algorithm.However, it is equivalent to use proposition 3.5(a) to compute ∂ α 2 G (0,0) and G α and to use proposition 3.5(b) to compute Φ γ .)These issues will be resolved in the next section.

Comparative study of reduced Green function computation methods
In this section.we aim to compare the methods to compute the reduced Green function, to demonstrate the accuracy of the asymptotic expansion and address the aforementioned issues.
First, we consider the case n = 1.Let Φ ( j) γ, * denote the result of computation via equation (29).Then, let Φ ( j) γ,k represent the sum of the jth derivative of first k + 1 terms in the asymptotic expansion (10).Also, let δ ( j) k denote the 'relative error' between Φ ( j) γ, * and Φ k under different γ and k are shown in figure 3.
We can infer from figure 3 that: (a) Φ (0) γ, * is accurate when x is small, but becomes numerically unstable as x increases, and the turning point decreases as γ increases; (b) Φ (0) γ,k approximates Φ γ better as k increases.However, we should take into account the complexity of the coefficients in Φ (0) γ,3 and its derivatives, which will make the computation less efficient.This suggests that we can approximate Φ γ by Φ (0) γ,2 when |x| is greater than the turning point x 0 mentioned above.x 0 can be dependent on γ, for example, x 0 = 10, 8 and 5 for γ = 0.6, 0.8 and 1.2 respectively.
Next, let us compare the values of δ 2 under different γ and j, as shown in figure 4.This implies that we can conform to the suggestion in the preceding paragraph when j is other  Figure 5 shows that there is no significant improvement in precision when we add the third and fourth term in the asymptotic expansion.Therefore, we can compute Φ γ by the asymptotic expansion Φ (0) γ,1 when |x| > x 0 and otherwise use the other way to compute ∂ α 2 G (0,0) and G α .Here, x 0 may be chosen as 10 for γ = 0.6 and 0.8, and 6 for 1.2. Figure 6 is the image of δ   under different γ and j, which justifies the choice of x 0 .This will bring a relative error of under 10 −3 .
In order to obtain the boundary data needed in the index function, we need to solve the forward problem.We discretize the temporal variable using the L2 − 1 σ method proposed in [1] when 0 < γ ⩽ 1 and a variant of the L2 method in [44] when 1 < γ < 2, with the step size set to 10 −3 .In the spatial domain, we employ the finite element method of order 2 with a mesh size of 1 400 diam(Ω).This is implemented in C++ code with the MFEM library [30].Finally, we use the fourth-order finite difference formulas to approximate the derivatives.
This example aims to test the validity of the sampling method under Case I in section 2, so we choose I 2,3 for our index function.The index function is calculated on a rectangular grid with a spacing of 0.03.In figure 9, we show the image of the normalized index function I 2,3 (x, 0), i.e. the index function divided by its maximum absolute value on sampling points, for different choices of γ.The location of sources is marked with black crosses in the figure.Note that g 1 (0)/g ′ 1 (0) and g 3 (0)/g ′ 3 (0) are as small as 1/50 and 1/6 respectively.It shows that the sampling method is capable of locating the sources with minor deviations even when the condition in Case I is satisfied weakly.Now we replace the underlying equation from equation ( 2) to (1), and present the graphs of normalized index function I 2,3 (x, 0) in figure 10.It turns out that the sampling method is still effective in this situation.
We apply the sampling method to the TFDE on Ω with a sampling interval of 0.02.The third-order derivative ∂ 3  2 f(x, t) of boundary data is shown in figure 11 and the normalized index function I 2,3 (x, 0) is shown in figure 12 with the location of sources marked with black crosses.This example illustrates that the sampling method can successfully find the location of sources, but fails to recover the shape and direction of them.
Example 3. Now we would like to test the sampling method under Case II in section 2. We retain the parameters in example 1 except that g 2 (t) = 3t, then apply the sampling method to   the TFDE on Ω with an interval of 0.03.The fourth-order derivative ∂ 4 2 f(x, t) of boundary data is shown in figure 13 and the normalized index function I 2,4 (x, 0) is shown in figure 14 with the location of sources marked with black crosses.The sampling method again gives satisfactory outcomes except for γ = 0.6, when the index function vaguely hints that a source may exist in the middle or upper-left side of the region.
Example 4. Finally we test an example where the source g is discontinuous in time.Simultaneously, we aim to determine whether the proposed sampling method can distinguish between closely located sources.
We set Ω to be the circular region as in examples 1 and 3, and set g(x, t) = 2 j =1 g j (t)δ xj (t), where x 1 = (−r, 0) and x 2 = (r, 0) with r determined below.The temporal components g j of the source conform to the normal distribution with mean 1 and variance 0.05, i.e. g 1 = g 2 ∼ N (1, 0.05).We provide its graph in figure 15.In this example, we fix the fractional order γ = 0.8, and consider decreasing source distance by setting r = 1, 0.75, 0.6 and 0.5.We apply the sampling method to the TFDE on Ω with an interval of 0.03.The normalized index function I 2,3 (x, 0) is shown in figure 16 with the location of sources marked with black crosses.This suggests that the proposed sampling method is applicable to discontinuous functions, and it successfully distinguishes between closely located sources.

Conclusion and remarks
In this article, we present a direct sampling method to determine the location of the source term in a TFDE.We introduce a family of index functions {I αβ } α,β∈N based on the property of the source term.We then focus on the computation of the index function when the dimension is not larger than 3, and derive asymptotic expansions to reduce round-off errors and improve the speed of computation.Numerical examples confirm the effectiveness of the sampling method.
Although the location of sources obtained from the sampling method may exhibit some deviation from the true location, it can still be used as an initial guess in an iterative method.Now we would like to provide some remarks and highlight some limitations in the choice of the index function.
• Consider the point source g(x, t) = m j =1 g j (t)δ xj (x).Define g j the zero extension of g j on R, then the boundary value f can be written as ˆR G (x, t; x j , s) g j (s) ds.
Roughly speaking, the index function I α,α+k amplifies the discontinuities of g j (k) over (−∞, T), so we can generalize Cases I and II in section 2 to piecewise regular functions g j .For example, if there is a jump discontinuity in g ′ j (t) at t = t 0 , then there will be a (positive or negative) peak near (x j , t 0 ) in I α,α+2 (x, t).
• In Case II we require g ′ j (0) > 0. This is because we would like to ensure that I α,α+2 has a maximum point (instead of a minimum point) near the point (x j , 0). • Theoretically speaking, as long as g j (0) = 0, we can find out x j using the index function I α,α+1 with a sufficiently large α in Case I.But due to numerical instability, we have to demand that g j (0) is not too small and use a smaller α.α = 2 is recommended.• It is difficult to clarify the assumptions in cases I and II.Examples 1 and 3 in section 4.3 show that both I 2,3 and I 2,4 can recover the sources with the same intensity g 1 and g 3 .On the other hand, these two cases are far from covering all cases, especially for smooth g j , such as g j (t) = e −1/t .
We end this section with a list of possible future research problems.
• Estimate the deviation of the maxima of the index function from the location of the source, at least for a certain collection of the source term.• Estimate the influence of the following factors on the accuracy of the index function: (a) the distance between sources, (b) the size of point-like sources or (c) gathering boundary data from equation (1) instead of equation ( 2).• Optimize the algorithm of the computation of the index function, so that the sampling method can be applied to real-life situations, which are large-scale and three-dimensional.• Find relationships between Φ n,γ for different n, or find the asymptotic expansion of Φ n,γ , where n > 3 is involved.This is an interesting topic, as for n > 3, E γ (−| • | 2 ) is no longer in L 2 .Therefore, Φ n,γ is a distribution.

Figure 8 .
Figure 8.The third derivative of the boundary data ∂ 3 2 f(x, t) for different γ.

Figure 11 .
Figure 11.The third derivative of the boundary data ∂ 3 2 f(x, t) for different γ.

Figure 13 .
Figure 13.The fourth derivative of the boundary data ∂ 4 2 f(x, t) for different γ.

Figure 15 .
Figure 15.Temporal component of the source g j in example 4.