Recovering the Potential and Order in One-Dimensional Time-Fractional Diffusion with Unknown Initial Condition and Source

This paper is concerned with an inverse problem of recovering a potential term and fractional order in a one-dimensional subdiffusion problem, which involves a Djrbashian-Caputo fractional derivative of order $\alpha\in(0,1)$ in time, from the lateral Cauchy data. In the model, we do not assume a full knowledge of the initial data and the source term, since they might be unavailable in some practical applications. We prove the unique recovery of the spatially-dependent potential coefficient and the order $\alpha$ of the derivation simultaneously from the measured trace data at one end point, when the model is equipped with a boundary excitation with a compact support away from $t=0$. One of the initial data and the source can also be uniquely determined, provided that the other is known. The analysis employs a representation of the solution and the time analyticity of the associated function. Further, we discuss a two-stage procedure, directly inspired by the analysis, for the numerical identification of the order and potential coefficient, and illustrate the feasibility of the recovery with several numerical experiments.

The model (1.1) has been studied extensively in the engineering, physical and mathematical literature due to its extraordinary capability for describing anomalous diffusion phenomena. It can be viewed as the macroscopic counterpart of continuous time random walk in which the waiting time between consecutive particle jumps follows a heavy-tailed distribution with a divergent mean, and the probability density function of the particle appearing at location x at time t > 0 satisfies a model of the form (1.1), in analogy with the classical diffusion equation for Brownian motion [30]. The model (1.1) inherits certain analytic properties of the latter, but also differs considerably due to the presence of the nonlocal fractional derivative term ∂ α t u: it has limited smoothing property in space and slow asymptotic decay in time [38,16]. The list of successful applications is long and still fast growing, including thermal diffusion in fractal domains [33], dispersion in a heterogeneous aquifer [1] and transport in column experiments [11] etc. See the comprehensive reviews [30,29] for the derivation of relevant mathematical models and many applications in physics and biology.
In this work, the inverse problem of interest is to recover the potential q ∈ Q := {q ∈ C(Ω) : q ≥ 0 in Ω} in the elliptic operator A and the order α of derivation from the boundary observational data at the left end point h(t) = u(0, t) for t ∈ [0, T ]. Note that in the model (1.1), besides the potential q and the order α, the space-dependent source term f and the initial data u 0 are both unknown. The situation that the initial data u 0 is inaccessible arises naturally, e.g., in heat conduction in high-temperature furnace [43]. To make the matter worse, only one single boundary observation data (at the left end) is available, which makes the inverse problem much more challenging both mathematically and numerically. The ability of choosing the boundary excitation g(t) is indispensable for the unique recovery, without which the identifiability generally does not hold, as indicated by example 3.1 below. In theorems 3.1 and 3.2, we present a uniqueness result for recovering the potential q in the operator A and the order α of derivation. The proof employs suitable solution representation in proposition 2.1, analyticity in time in proposition 2.2 and Gel'fand-Levitan theory. Further, we discuss the numerical reconstruction by a two-stage procedure inspired directly by the analysis: at stage i, we numerically continuate the boundary observation data h(t) by rational functions, and at stage ii, we perform a standard least-squares procedure to recover the potential q using the conjugate gradient method (with proper early stopping). The simulation study with exact data indicates that the recovery is feasible. The uniqueness result, the two-stage recovery procedure and the numerical verification of the recovery represent the main contributions of this work.
Next we situate the work in existing literature. The recovery of the space-dependent potential q in the classical diffusion equation from lateral Cauchy data has been extensively discussed, and several uniqueness results have been obtained [34,31,41]. The study on related inverse problems for time-fractional models is of more recent origin, starting from [7] (see [19] for an early tutorial) and there are a few works on recovering a spatially dependent potential from lateral Cauchy data [36,37,44,22]. Rundell and Yamamoto [36] showed that the lateral Cauchy data can uniquely determine the spectral data when u 0 ≡ f ≡ 0, and proved the uniqueness of the potential q by the classical Gel'fand-Levitan theory. They also proposed a recovery procedure based on Newton's method and empirically studied the singular value spectrum of the linearized forward map, showing the severe ill-posed nature of the inverse problem. Later, they [37] relaxed the regularity condition on the boundary excitation g(t) (in a suitable Sobolev space in time). Recently, Jing and Yamamoto [22] proved the identifiability of multiple parameters (including order, spatially dependent potential, initial value and Robin coefficients in the boundary condition) simultaneously in the one-dimensional subdiffusion / diffusion-wave (i.e., α ∈ (0, 2)) equation with a zero boundary condition and source, excited by a nontrivial initial condition from the lateral Cauchy data at both end points (cf. remark 3.2 for details). See also the work [44] for relevant results in the diffusion wave case; and [21] for the case of a Robin boundary conditions. In all these existing works, the initial condition / source is assumed to be fully known, so that the forward map is well defined, which differs from the current work. There are two closely related inverse problems to the concerned one. (i) is to recover the spatially dependent potential q from the terminal data u(T ) [47,23,20], which enjoys much better stability estimates (e.g., local Lipschitz stability) and effective iterative algorithms for numerical recovery, e.g., fixed point iterations. (ii) is to recover a time-dependent potential q(t) from time-dependent observations [15] (or a time-dependent source from observation at one point [28]), which behaves similarly to (i) due to the directional alignment of the unknown and observations. The rest of the paper is organized as follows. In section 2, we derive a crucial representation of the solution to the direct problem (1.1). Then in section 3, we prove the unique recovery of the order α and the potential q. In section 4, we describe a two-stage numerical algorithm for recovering the potential q. Last, we present several numerical experiments to show the feasibility of the simultaneous recovery in section 5. Throughout, the notation c denotes a generic constant which may differ at each occurrence, and (·, ·) denotes the standard L 2 (Ω) inner product (or duality pairing).
2 Well-posedness of the direct problem In this section, we collect several preliminary results on the direct problem (1.1), especially the solution representation, which will play an important role in the study.

Preliminaries
First we describe several preliminary results that will be used in deriving the solution representation. We use extensively the two-parameter Mittag-Leffler function E α,β (z) defined by (see, e.g., [10] and [16, Section 3.1]) It is an entire function of order 1 α and type one. The following properties hold [16, Section 3.1]. Lemma 2.1. For any α ∈ (0, 2) and β ∈ R, the following statements hold.
(i) For any ϕ ∈ ( α 2 π, min(π, απ)), the following asymptotics hold (ii) For any λ > 0, the following Laplace transform relation holds (iii) The following differentiation formula holds for any λ ∈ C Next, we introduce Bochner-Sobolev spaces W α,p (0, T ; X), for a UMD space X (see [13,Chapter 4] for the definition of UMD spaces, which include Sobolev spaces W s,p (Ω) with s ≥ 0 and 1 < q < ∞). For any s ≥ 0 and 1 ≤ p < ∞, we denote by W s,p (0, T ; X) the space of functions v : (0, T ) → X, with the norm defined by complex interpolation. Equivalently, the space is equipped with the quotient norm where the infimum is taken over all possible v that extend v from (0, T ) to R, and F denotes the Fourier transform. In case that X = R, we denote W s,p (0, T ; R) by W s,p (0, T ) for convenience. The next lemma provides a norm equivalence result [17,Lemma 2.3].

Well-posedness of the direct problem
Now we study the direct problem (1.1), especially the solution representation. One distinct feature of problem (1.1) is that it involves a nonzero Neumann boundary condition, which has not been extensively studied in the literature ( [37,25] for relevant works). Following [37], we exploit the one-dimensional nature of problem (1.1), and derive a series representation of the solution u. The derivation is based on the standard separation of variable technique (see, e.g., [38] and [16,Section 6.2]). Specifically, let A be the realization of the elliptic operator −A in L 2 (Ω), with its domain Let {(λ n , ϕ n )} ∞ n=1 be the eigenpairs of the operator A, i.e., −Aϕ n = λ n ϕ n , in Ω, By the standard Sturm-Liouville theory [27], the spectrum of the operator A consists of a strictly increasing sequence of positive eigenvalues {λ n } ∞ n=1 and the associated eigenfunctions {ϕ n } ∞ n=1 can be chosen to form an orthonormal basis of the space L 2 (Ω). By means of Liouville transformation, we deduce that the eigenvalues λ n grow asymptotically as O(n 2 ) [27] (also known from Weyl's law [45]): The unnormalized eigenfunctions ϕ n satisfy the following asymptotics [27, Section 2 of Chapter 1]: This estimate implies that the L 2 (Ω)-orthonormal eigenfunctions are uniformly bounded. Then we define the fractional power A s , s ≥ 0, by with its domain {v ∈ L 2 (Ω) : A s v ∈ L 2 (Ω)}, and the associated graph norm · D(A s ) given by With these preliminaries, we can now study the direct problem (1.1). In view of the linearity of the problem, we may split the solution u into two parts: u = u i + u b , with u i and u b solving respectively With the eigenexpansion {(λ n , ϕ n )} ∞ n=1 , the solution u i can be represented by (see e.g., [38] and [16, Section 6.2]) When the source f is time independent, by Lemma 2.1(iii), we have Further, we have the following a priori estimate on the solution u i .
Next we turn to the a representation of the solution u b . We need the following identity.
Proof. A variant of the identity can be found in [37,Section 2], and we recap the proof only for completeness. We denote the integral on the left hand side by S. Then changing the order of integration gives where the last line follows from the change of variables η = t − s. Moreover, by the definition of E α,α (−λ n t α ) and applying termwise integration [16, (3.5)], we have 1 Thus, by integration by parts and lemma 2.1(iii), we obtain This completes the proof of the lemma. Now we can derive a representation of the solution u b , corresponding to nonzero boundary conditions. The derivation essentially exploits the one-dimensional nature of the problem. Proposition 2.1. Let α ∈ (0, 1) and p > 4/(3α) . Suppose g 1 ∈ L p (0, T ), g 2 ∈ W α,p (0, T ) and g 2 (0) = 0. Then the solution u to the following initial boundary value problem can be represented by Proof. The derivation proceeds by homogenizing the boundary conditions. First, we assume g 1 , Next we evaluate the boundary conditions: Thus, the function v satisfies ). Since the function v satisfies homogeneous boundary conditions, by the standard separation of variable technique (see e.g., [38] or [16, Section 6.1]), it can be represented by Next we simplify the second term. By the definition off and using lemma 2.4, we deduce Next we evaluate the last two terms in the square bracket. By Green's identity, for any φ ∈ H 2 (Ω), Using this identity, the two terms in the square bracket can be evaluated as Consequently, we arrive at
Proof. The assertion is direct from proposition 2.1 and the identity (2.4). Indeed, we have ). This directly leads to the desired identity.
Remark 2.1. Let the solution operators E(t) and F (t) be defined by Then the solution u can be formally represented with where δ 0 (x) denotes the Dirac delta function concentrated at x = 0. Indeed, we can expand δ 0 (x) in terms of the . Then substituting the identity and collecting the terms lead to the desired identity. Alternatively, the solution can also be represented concisely using the fractional θ functions; see [16, Section 7.1] for relevant discussions.
Next we study the convergence of the series in corollary 2.1. First, note that for f ∈ D(A −s ) with 0 ≤ s < 3 4 , the constant ρ 0 is well defined. Indeed, by (2.2) and the Cauchy-Schwarz inequality, and asymptotics for the eigenvalues {λ n } ∞ n=1 , we have For the remaining series, we give two analyticity results concerning the following two auxiliary functions: They arise naturally in the uniqueness proof, and the analyticity plays an important role in section 3.
(ii) The Laplace transforms of h 1 (t) and h 2 (t) exist and are given respectively by Proof. By lemma 2.1(i), there exists a constant c and θ 0 > 0 such that , ∀x ∈ Ω, z ∈ Σ, and by the Weyl's asymptotics of the eigenvalues λ n and the Cauchy-Schwarz inequality, we have we deduce that the series in analytic in t. The analyticity of h 2 (t) follows similarly as h 1 (t), but with the following estimate from lemma 2.1(i): there exists a constant c and θ 0 > 0 such that Then by Weyl's asymptotics of the eigenvalues λ n , for any γ ∈ ( 1 2 , 1), Since E α,α (−λ n z α ) is analytic in z ∈ Σ, we deduce that the series in analytic in t. These discussions show assertion (i). Next, for any t > 0, both series converge uniformly in [t, ∞), and there holds and the function e −t (z) t −α is integrable in t over (0, ∞) for any fixed z with (z) > 0. By Lebesgue dominated convergence theorem, we can take Laplace transform termwise and by lemma 2.1(ii), we obtain Thus, the Laplace transform of h 1 (t) exists. The Laplace transform of h 2 (t) follows similarly from the estimate Then by Lebesgue dominated convergence theorem and lemma 2.1(ii), we obtain for (z) > 0 This shows assertion (ii), and completes the proof of the proposition.

Uniqueness
In this section, we study the uniqueness of the inverse problem: given the observation h(t) = u(0, t) at the left end point x = 0, can we uniquely determine the potential q and the order α? Note that without any restriction on the boundary data g, the desired uniqueness result does not hold. This is illustrated by the following example with a zero excitation g ≡ 0.
Now we proceed to the uniqueness. The proof is split into two steps, and both steps rely on the time analyticity result in proposition 2.2. The first step is concerned with the unique recovery of the order α and partial information of the unknown initial data u 0 / source f . The notation K denotes the set {k ∈ N : ρ k = 0}, i.e., the support of the sequence (ρ 0 , ρ 1 , . . .), with the constants ρ k defined in corollary 2.1, and the setK is defined similarly. It is worth noting that the set K is not a priori known, since the elliptic operator A is not fully known (due to the unknown potential q). The condition K = ∅ holds as long as h(t) ≡ constant on [0, T 2 ], and thus it is very mild.

Remark 3.2.
There have been several works on identifying multiple parameters from one single observation [46,24,22]. The recent work [22] is closest to the current one in some sense, which is concerned with the following model on Ω = (0, 1), with α ∈ (0, 2), The inverse problem is to recover u 0 , q, α, h and H from two boundary observations, i.e., u(0, t) and u(1, t). They proved the uniqueness of the recovery under the following condition (for α ∈ (0, 1)): (u 0 , ϕ n ) = 0, for all n ∈ N [22, Theorem 1]. This condition assumes that all the eigenmodes of the initial value u 0 should be nonzero, which is generally restrictive, and can be relaxed using multiple initial conditions [22,Theorem 1']. In contrast, theorem 3.2 relies on the nonzero boundary excitation g(t) for the potential recovery, and thus avoids this assumption. Remark 3.3. There are several potential extensions of the stated uniqueness results. (1) The results hold also for the multi-term time-fractional model, which involves multiple time fractional derivatives, i.e., the term ∂ α t u in the model (1.1) is replaced by N i=1 r i ∂ αi t u, with r i > 0 and 0 < α 1 < . . . < α N < 1. Then the weights r i and α i are uniquely determined, provided that h(0) = 0. (2) One can uniquely determine the diffusion coefficient a when the potential q is known, by a different version of Gel'fand-Levitan theory [7].
(3) The boundary conditions can be of more general Sturm-Liouville form. Then the Robin coefficients in the boundary conditions can also be determined uniquely from lateral Cauchy data [22], cf. remark 3.2.
The preceding analysis indicates that the both steps rely essentially on unique continuation, which is well known to be severe ill-conditioned. A natural question is how the fractional paradigm actually affects the degree of ill-conditioning, measured in terms of the asymptotic decay rate of the singular value spectrum of the associated (linearized) forward map. This issue has been numerically studied for several inverse problems in [19]; see also [36] for the inverse potential problem. However, a theoretical analysis in the context of potential recovery from lateral Cauchy data is still unavailable.

Reconstruction algorithm
Now we describe an algorithm for simultaneously recovering the potential q, the order α, and also u 0 , under the assumption f ≡ 0 (or also f , if u 0 ≡ 0). The procedure is directly inspired by the uniqueness proof, and consists of two steps.

Step 1: order determination and numerical continuation
In the first step, we determine the fractional order α and numerically continuate the trace data h(t) from [0, T 1 ] to the whole interval [0, T ] (to assist the recovery of the potential q). We discuss the two issues separately. The recovery of the order α cannot be carried out in the usual manner by means of least-squares fitting, since the problem data in the direct problem (1.1) over the interval [0, T 1 ] is not fully known. The next result suggests one possible recovery formula for the order α from the small time asymptotics of the observation h(t), under suitable smoothness condition u 0 and f , if the function Au 0 + f does not vanish at x = 0.
Proof. By the definition of the Mittag-Leffler function E α,1 (z), we have This and the solution representation from corollary 2.1 lead to We denote the last sum by I. Since the eigenfunctions {ϕ n } ∞ n=1 forms an orthonormal basis in L 2 (Ω), by integration by parts, we have and ∞ n=1 (f, ϕ n )ϕ n (x) = f (x). By lemma 2.1(i), we bound the sum I by where the last inequality follows from the conditions s ∈ ( 1 4 , 5 4 ] and ∈ (0, s − 1 4 ). Combining the preceding estimates with the Sobolev embedding theorem directly shows the assertion.

By proposition 4.1, under mild conditions, the trace data h(t) satisfies
This motivates a simple procedure: minimize over α (and c 0 and c 1 ) the following objective for some t 0 > 0 sufficiently close to zero. The minimization can be carried out by any stand-alone algorithms, e.g., gradient descent, and Newton method. Note that it is important to take t 0 sufficiently close to zero so that the term o(t α ) is indeed negligible. Next, we numerically continuate the given data h(t) from the interval [0, T 1 ] to [T 1 , T ], in order to extract the combined information on u 0 and f . Mathematically, this amounts to recovering {(ρ k , λ k )} k∈K (with the index set K defined in section 3) from h: By theorem 3.1, {(ρ k , λ k )} k∈K can indeed be uniquely determined by h(t), t ∈ [0, T 1 ]. This problem is also known as an infinite-dimensional spectral estimation problem for α = 1, for which the issue is to recover (ρ k , λ k ) of an exponential family [4] and there are several efficient methods for recovery, e.g., matrix pencil method [12] and MUSIC (MUltiple SIgnal Classification) [39]. However, for α = 1, to the best of our knowledge, there is no known analogue of these methods. This is essentially due to the inequality ∞). Instead, we resort to the classical rational approximation for numerical continuation, i.e., where r ∈ N is the polynomial order. The approximation h r (t) can be constructed efficiently when h(t) is accurate using the AAA algorithm [32], despite the well-known ill-posed nature of analytic continuation. This choice is in part motivated by the fact that the function E α,1 (−λt α ) admits excellent rational approximations [16,Theorem 3.6]. Our numerical experiments indicate that the procedure is indeed viable for exact data.

4.2
Step 2: recovering q (and u 0 ) by iterative regularization With the analytic continuation in Step 1, we can proceed to the reconstruction of the potential q, as in the proof of theorem 3.2. Specifically, leth which represents the reduced data for the boundary excitation g only (supported on the interval [T 1 , T ], by construction). This naturally motivates approximately minimizing with F (q) = u(q)(0, t), where u(q) denotes the solution to the direct problem (1.1) corresponding to the elliptic operator A, with u 0 ≡ f ≡ 0 and given g. The map F is nonlinear, and one may apply standard iterative regularization methods [8,14], e.g., (nonlinear) conjugate gradient method. In the numerical experiments, we employ the conjugate gradient method [2], which generally enjoys fast convergence.
Once the potential q is determined, one can also attempt recovering the initial data u 0 from the observation h(t) over the interval [0, T 1 ], if f ≡ 0. This can be achieved by approximately minimizing with F (u 0 ) = u(0, t), where u denotes the solution to the direct problem (1.1) corresponding to the elliptic operator A with the recovered q, and f ≡ 0 (over the interval (0, T 1 )). The optimization can be carried out efficiently by standard gradient type methods. In practice, the gradients of the functionals J(q) and J(u 0 ) can be computed efficiently by the adjoint technique. We provide relevant details in the appendix.

Numerical results and discussions
Now we present some numerical results to illustrate the feasibility of simultaneously recovering the coefficient q and the fractional order α, without fully knowing the direct problem (1.1). The domain Ω is taken to be the unit interval Ω = [0, 1], and the final time T = 1, and T 1 = 0.5. The direct and adjoint problems are all discretized by the standard continuous piecewise linear Galerkin method in space, and backward Euler convolution quadrature in time [18]. The domain Ω is divided into M subintervals each of width 1/M . For the inversion step, we take M = 200 and N = 2000. The exact data h † on the lateral boundary (0, T ) is obtained by solving the direct problem (1.1) on a finer mesh. It is known that due to the severe ill-conditioning of the inverse problem, the numerical recovery in the presence of data noise is very challenging. Indeed, for the inverse potential problem, it was observed numerically in [36] that there are only a few significant singular values, and this is also partly confirmed by [40]. This is further complicated by the unknown problem data in the present context. Thus, our experiments below focus on exact data. We illustrate on the following two settings, with u 0 and f being unknown: x) and f ≡ 0, and g = χ [T1,T ] ; (ii) a ≡ 1, q † = min(x, 1 − x), u 0 = cos( 3π 2 x) and f ≡ 0, and g = χ [T1,T ] , where χ S denotes the characteristic function of the set S. The initial condition u 0 is taken to be in D(A) so that the asymptotic expansion in proposition 4.1 is indeed valid. Case (i) involves a smooth potential, and case (ii) a nonsmooth potential.
First we study the recovery of the fractional order α using a least-squares fitting as described in section 4. This procedure relies on the validity of the asymptotic expansion in proposition 4.1. The recovered orders are presented in Table 1, where the minimization is carried out by the L-BFGS-B [5], with the box constraint α ∈ [0, 1], using the public implementation https://ww2.mathworks.cn/matlabcentral/ fileexchange/35104-lbfgsb-l-bfgs-b-mex-wrapper (last accessed on May 20, 2021). Note that the least-squares functional is fraught with many local minimum, and a good initial guess for α is needed in order to recover the correct order. It is observed that the accuracy of the recovery tends to improve as the interval (0, t 0 ) used in the least-squares formulation shrinks, since the model function in proposition 4.1 represents an increasingly better approximation as t → 0 + . Further, as α increases, the size of the interval (0, t 0 ) can be increased without sacrificing the accuracy of the recovery since the asymptotic expansion is then valid in a larger neighborhood. Thus one may conclude that with the interval (0, t 0 ) chosen properly (and of course only for very accurate data), the order α can indeed be recovered reliably by the least-squares fitting. These observations hold for both cases (i) and (ii), and thus the smoothness of the potential q does not seem to influence much the recovery of the order α.  One step of the recovery procedure is analytic continuation, extending the observation data h by a rational model h r from the interval [0, T 1 ] to [T 1 , T ]. This step extracts relevant information from unknown initial condition u 0 (and source f ), and plays a central role in formulating the optimization problem for recovering the potential q. This is illustrated in Fig. 1 for the two cases at α = 0.5, where the rational approximation h r is constructed by the AAA algorithm [32] using the MATLAB implementation given therein with a tolerance 1e-9, with the resulting h r of degree r = 11. The pointwise error is evaluated against the ground-truth h * (i.e., h * = u(0, t), t ∈ [0, T ], u being the solution of the direct problem (1.1) with g ≡ 0) over the interval [T 1 , T ]. (Numerically, larger tolerances, e.g., 1e-6, can still give an accurate approximation.) Clearly, h r does give a fairly accurate approximation to h * , and the accuracy degrades as one moves away from the interpolating interval [0, T 1 ]. It is noted that the continuation results for other cases exhibit very similar behavior. Thus, the rational approximation is a very effective approach for analytic continuation when exact data is available.
The reconstructions of the potential q by the conjugate gradient (CG) method, based on the reduced datah(t), are shown in Fig. 2 (with exact order) and Table 2. The maximum number of CG iterations is fixed at 200, and it is stopped so that the error is smallest possible. Throughout, for a reconstructionq, we measure the residual r(q) and the L 2 error e(q), defined respectively by where q † denotes the exact potential. The accuracy of the reconstructions actually does not depend on very much on the order α, and all the reconstructions represent a reasonable but not perfect approximation to the true potential q † . This observation is consistent with prior numerical results for similar problems [40,36], and might be attributed to severe ill-conditioning of the inverse problem. The CG method can steadily decreases the value of the objective (i.e., the residual r), with the first few steps converging fairly rapidly and then slowing down considerably. Nonetheless, the error e trajectory exhibits an unusual oscillating pattern during the iteration: the error e first decreases, and then increases and then further decreases again, and there is also a flat region for which the error e stays nearly constant. This behavior differs drastically from the typical steady error convergence observed for other inverse problems, e.g., inverse source problems [17]. The precise mechanism of the behavior remains elusive. It is worth noting that all these changes occur after the residual r reaches a relatively small magnitude (and flat region), indicating a potential numerical "identifiability" issue, despite the uniqueness in theorem 3.2. This also indicates that in the presence of data noise, the magnitude of the noise has to be very small so that not to wash away these tiny transitions in order to have a fair recovery.
In the current context, the order α is numerically recovered, which incurs inevitable errors. This error can potentially impact the subsequent inversion of the potential q. To examine the influence, we perturb the order α in the optimization problem (4.1) by δα, and repeat the numerical experiments with α + δα. The results are summarized in Table 2, where k * denotes the iteration index at which the error is smallest, and e * and r * denote the corresponding error and residual. The presence of perturbation δα does not affect very much the attainable accuracy, although the error e * increases steadily with the perturbation δα; and generally it takes fewer iterations to reach the optimal accuracy. This observation is largely valid for both cases with all fractional orders under consideration.
Last, we examine the recovery of the initial data u 0 , using the recovered potential q by the CG method (terminated after 200 iterations). The related numerical results are summarized in Table 3 and Fig. 3, where we have assumed that the order α has been estimated reliably. Interestingly, despite the inaccuracy of the recovered potential q, the initial data u 0 can still be recovered with a good accuracy, for all fractional orders. Further, the convergence behavior of the CG method agrees well with that for other inverse problems (but contrasts sharply with that for potential recovery): the method decreases the residual e steadily, and the reconstruction error e exhibits a typical semiconvergence phenomenon, i.e., the iterates first converge and then diverge, due to the inherent ill-posed nature of the inverse problem.
In sum, the numerical experiments demonstrate the following empirical observations: (1) The order α can be recovered from the observation data by a least-squares procedure; (2) the rational function approach  represents a simple method for analytically continuate the data; (3) The CG method can produce fair approximations to the potential q, even under small perturbations of the order α, partly confirming the feasibility of the recovery, but the convergence behavior of the algorithm exhibits an unusual oscillating feature that remains to be further examined; (4) the CG method can produce good approximations of the initial data u 0 , based on the recovered potential q. In particular, the experiments show that the simultaneous recovery of the order, potential, and initial data (or source) is indeed feasible, provided that accurate lateral Cauchy data is available, thereby corroborating the uniqueness results in section 3.

Concluding remarks
In this work, we have studied an inverse problem of simultaneously recovering the fractional order and the space-dependent potential in a one-dimensional subdiffusion model from the observation data at the end point, when the initial data / source is not fully known. We proved that both order and potential can be uniquely determined, if the Neumann boundary condition satisfies a mild condition. Further, one of the space-dependent source or initial condition can be uniquely determined, if the other is known. The analysis lends itself to an effective two-stage reconstruction algorithm. Numerical results also show the feasibility of the recovery. There are many related theoretical and numerical issues awaiting further research. First, it is of interest to extend the results to the case of a time-dependent potential. One obstacle in the extension is that the time-dependence of the potential precludes a direct application of the separation of variable technique, an important tool in the current analysis. Second, it is natural to analyze more complex subdiffusion models, e.g., multi-term and variable orders (e.g., α(t), α(x) or α(x, t)). We believe that the results remain largely valid for the multi-term case. However, for variable-order models, the solution theory is still far from complete, and substantially new analytical tools are needed. Third, the extension to the multi-dimensional case is very challenging, and requires more data for a unique determination, e.g., restricted Neumannto-Dirichlet map [6] or one specially designed excitation [24]. Fourth and last, the design and analysis of relevant reconstruction algorithms can depart enormously from the more traditional (penalized) least-squares approach. The latter might not be directly applicable, due to the presence of unknown problem data (and thus the very forward model in the least-squares formulation is also unknown).
A The computation of the gradients J (q) and J (u 0 ) To apply the conjugate gradient method, one has to compute the gradient. This can be done efficiently using the adjoint technique. Below we give relevant computation details for completeness. We have the following representations of the gradients J (q) and J (u 0 ). Note that the adjoint problem for v and w satisfies a nonlocal terminal condition. The notation t I 1−α T v(t) and t R ∂ α T v are defined by [16] t I 1−α with v ≡ v(q) and w solving respectively Then taking φ = v ∈ X in (A.1) and φ = u δ ∈ X in (A.2), using the following integration by parts formula (see, e.g., [26, p.   Meanwhile, the space-time weak formulation for the adjoint solution w is given by This gives the expression of J (u 0 ).