It is well known (see, for example, [1] - [4] ) that Brownian motion, discovered in the first half of the 19th century, models motions of molecules in gases, electrons in semiconductors, neutrons in nuclear reactors, and much more. The main difference between Brownian motion and processes obeying Newton's laws is that the diffusion packet spreads according to the law t 1/2 (and not like t). The subdiffusion process is characterized by a fractional exponent ρ ∈ (0, 1), which is included in the diffusion equation as the order of the fractional time derivative.
The theory of differential equations with fractional derivatives has gained considerable popularity and importance in the past few decades, mainly due to its applications in numerous seemingly distant fields of science and technology (see, for example, [1] - [8] ). In turn, the mathematical aspects of fractional differential equations and methods for solving them have been studied by many authors (see, for example, [9] - [33] ).
By inverse problems in the theory of partial differential equations are commonly called problems in which, together with solving a differential equation, it is also necessary to determine a coefficient(s) of the equation or/and the right side (source function). Naturally, in this case additional information should be given to find a new unknown function. Note that the interest in studying inverse problems for equations of mathematical physics is due to the importance of their applications in many branches of modern science, including mechanics, seismology, medical tomography, epidemics, and geophysics, just to mention a few. A significant number of studies are devoted to inverse problems of determining the right-hand side of subdiffusion equations (see, for example, [1] , [14] - [22] and references therein).
The present paper is devoted to the other important type of inverse problems, namely to determining of the order of fractional derivative in a subdiffusion equation, which is considered to govern the anomaly of diffusion. More precisely, this inverse problems is the determination of the unknown order of time-derivative in order to match available data such as u(x0, t), 0 < t < T, at a monitoring point x0 ∈ Ω. One of the practical example is a modeling of COVID-19 outbreak. The data [23] presented by Johns Hopkins University about the outbreak from different countries seem to show fractional order dynamical processes, in which the identification of fractional order rate of change is a key issue [24, 25] .
The problem of identification of fractional order of the model was considered by some researchers. Note that all the publications assumed the fractional derivative of order 0 < ρ < 1 in the sense of Caputo and studied mainly the uniqueness problem. In paper [26] by J. Cheng et al. an inverse problem for determining the order of the Caputo fractional derivative and a coefficient of one-dimensional time-fractional diffusion equation is studied. The authors attached the homogeneous Neumann boundary condition and the initial value given by the Dirac delta function. They proved that the order of derivative and the unknown coefficient are uniquely determined by the known datum u(0, t), 0 < t < T . The uniqueness of a solution of the two parameter inverse problem is considered in paper [27] by Tatar and Ulusoy for the differential equation
The multi-term time-fractional diffusion equation and distributed order fractional diffusion equations considered in papers Li et al. [28] and [29] , correspondingly.
In paper [30] by X. Zheng et al. the authors tried to solve the most difficult problem of determining the variable order of the Caputo fractional differentiation. In this work, as in many other papers, only the question of uniqueness is considered. But in our opinion, Lemma 4.1 of this paper is questionable, since there exist functions (see, for example, [35] ) whose Fourier series converge to zero in a certain region, but not all Fourier coefficients are zero.
The following two papers [31] and [32] deals with the existence problem. J. Janno [31] considered a one-dimensional time-fractional diffusion equation with Caputo derivatives. Giving an extra boundary condition Bu(·, t) = h(t), 0 < t < T the author succeeded to prove the existence theorem for determining the order of the derivative and the kernel of the integral operator in the equation. The complexity of the proof of the existence can be seen from the statement of corresponding theorem (Theorem 7. 2 is formulated on more than one journal page). In the paper of Hatano et al. [32] the equation ∂ ρ t u = ∆u, where ∆ is the Laplace operator, is considered with the Dirichlet boundary condition and the initial function ϕ(x). They proved that if ϕ ∈ C ∞ 0 (Ω) and ∆ϕ(x0) = 0, then
In the recent survey paper [33] by Z. Li et al. in the section of Open Problems they noted: "The studies on inverse problems of the recovery of the fractional orders are far from satisfactory since all the publications either assumed the homogeneous boundary condition or studied this inverse problem by the measurement in t ∈ (0, ∞). It would be interesting to investigate inverse problem by the value of the solution at a fixed time as the observation data".
In the present work we address this problem. Namely, as follows from our main result, in the case of the initial-boundary (Neumann) value problem for the equation with the Riemann-Liouville fractional derivative
where t0 ≥ 1 is an observation time, recovers the order ρ ∈ (0, 1), and if we have two pairs of solutions {u1(x, t), ρ1} and {u2(x, t), ρ2}, then u1(x, t) ≡ u2(x, t) and ρ1 = ρ2.
The paper is organized as follows. In the next section we formulate the main result. In Section 3 we prove the existence of a unique solution to the forward problem. This result will be used to prove the main result in Section 4. Throughout the paper we assume that the fractional order ρ of the main equation is constant. The solution of the forward problem is obtained under this assumption. The Cauchy problem in the case of piece-wise constant ρ was studied in [34] when x ∈ R N . Modifying the result of this paper to the case of bounded domain Ω, one can extend the main result presented in this paper to the case of piece-wise constant ρ, as well.
Let Ω be an arbitrary N -dimensional domain with twice differentiable boundary ∂Ω. Namely, the functions, defining the boundary equation in local coordinates, are continuously twice differentiable. Let the second order differential operator
for all x ∈ Ω and ξi, where a = const > 0 and Dj u = ∂u ∂x j , j = 1, . . . , N. Consider the spectral problem
where n is an external normal to the surface ∂Ω. It is known (see, for example, [36] , p. 100, [37] , p. 111), that if c(x) ≥ 0 and
where ⌊a⌋ stands for the integer part of a, then the corresponding inverse operator is compact, i.e. spectral problem (2.1) -(2.2) has a complete in L2(Ω) set of orthonormal eigenfunctions {v k (x)} and a countable set of nonnegative eigenvalues {λ k }.
The fractional part of our equation will be defined through the Riemann-Liouville fractional derivative ∂ ρ t of order 0 < ρ < 1. To define the Riemann-Liouville fractional derivative, one can define the fractional integration of order ρ < 0 of a function f defined on [0, ∞) by the formula
provided the right-hand side exists. Here Γ(z) is Euler's gamma function. Using this definition one can define the Riemann-Liouville fractional derivative of order ρ, k − 1 < ρ ≤ k, k ∈ N, as (see, for example, [9, 10] )
Note that if ρ = k, then fractional derivative coincides with the ordinary classical derivative of order k :
. Let 0 < ρ < 1 be an unknown number to be determined. Consider the initial-boundary value problem
Under some conditions on initial function ϕ the solution of this problem exists and is unique. This solution obviously depends on ρ. The purpose of this paper is not only to find the solution u(x, t), but also to determine the order ρ ∈ (0, 1) of the time derivative. To do this one needs an extra condition. As was mentioned above, different types of such conditions were considered by a number of authors. We formulate our inverse problem in the following way. Let ω ∈ L2(Ω) be a weight function with a property ||ω|| L 2 (Ω) = 1. Is it possible to determine the order of time fractional derivative 0 < ρ < 1 with the additional information Ω u(x, t0)ω(x)dx = d0, at a fixed time instant t0 ≥ 1? This integral can be considered as the average distribution of the quantity u(t, x) over the region Ω at the time instant t = t0 with the weight function ω.
Below, using the classical Fourier method, we give a positive answer to this question in the case when a weight function is equal to the first eigenfunction of spectral problem (2.1), (2.2): ω(x) = v1(x), i.e. when the latter integral has the form
The quantity f (ρ; t0) is, in fact, the projection of the solution u(x, t0) onto the first eigenfunction, and as is shown below, using the Fourier method, under certain conditions, f (ρ; t0) has a simple form f (ρ; t0) = dt ρ−1 0 , with a constant d depending on ρ. We call the initial-boundary value problem (2.4) -(2.6) (for an arbitrary, but fixed 0 < ρ < 1) the forward problem. The initial-boundary value problem (2.4) -(2.6) for an unknown ρ, together with extra condition (2.7), is called an inverse problem. A pair {u(x, t) , ρ} of the function u(x, t) and the parameter ρ with the properties (1) ρ ∈ (0, 1),
and satisfying all the conditions of problem (2.4) -(2.7) in the classical sense is called a classical solution of inverse problem (2.4) -(2.7).
We will also call the classical solution simply the solution to the inverse problem. We draw attention to the fact, that in the definition of the classical solution the requirement of continuity in the closed domain of all derivatives included in equation (2.4) is not significant. However, on the one hand, the uniqueness of such a solution is proved quite simply, and on the other hand, the solution found by the Fourier method satisfies the above conditions.
Further, let us denote by ϕj the Fourier coefficients of the function ϕ(x) with respect to the system of eigenfunctions {v k (x)}, defined as a scalar product on L2(Ω), i.e. ϕj = (ϕ, vj). Let D α stands for D α 1 1 · · · D α N N . For the seek of simplicity, we describe the proposed method, which is based on the classical Fourier method, for finding the order of fractional differentiation in the case of λ1 = 0 and ϕ1 = 0. Otherwise, the method becomes technically combersome. Now we formulate the main result of this paper. 
(1) Conditions (2.8)-(2.10) are standard for the existence of a solution to the forward problem (see for example, [21] ). Condition λ1 = 0 of the theorem is satisfied, for example, in the case of the Neumann condition on the boundary for the Laplace operator.
(2) Theorem defines the unique ρ from (2.7). Hence, if we define the integral (2.7) at another time instant t1 and get a new ρ1, i.e. f (ρ1; t1) = d1, then from the equality f (ρ1; t0) = d0, by virtue of the theorem, we obtain ρ1 = ρ. In this section we prove existence and uniqueness of the solution of the forward problem by the Fourier method. In accordance with the Fourier method, we will look for a solution to problem (2.4) -(2.6) in the form of a series:
where functions Tj (t) are solutions to the Cauchy type problem
The unique solution of problem (3.2) has the form (see, for example, [10] , p. 16)
where Eρ,µ is the Mittag-Leffler function
. 
which absolutely and uniformly converges on x ∈Ω for each t ∈ (0, T ].
Proof. The uniqueness of the solution can be proved by the standard technique based on completeness of the set of eigenfunctions {v k (x)} in L2(Ω) (see, for example, [21] ) To prove the existence we need to introduce for any real number τ an operatorÂ τ , acting in L2(Ω) in the following wayÂ
Obviously, the operatorÂ τ with the domain of definition
is selfadjoint. If we denote by A the operator in L2(Ω), acting as
and with the domain of definition
then it is not hard to show, that the operatorÂ ≡Â 1 is the selfadjoint extension of the operator A in L2(Ω). In the same way one can define the operator (Â + I) τ , where I is the identity operator in L2(Ω). Further, we use the following lemma (see [39] , p. 453):
Then for any multi-index α satisfying |α| ≤ 2 the operator D α (Â + I) −σ (completely) continuously maps the space L2(Ω) into C(Ω), and moreover, the following estimate holds
Let |α| ≤ 2. First we prove that one can validly apply the operators D α and ∂ ρ t to the series in (3.4) term-by-term. Suppose that the function ϕ(x) satisfies the following condition for some τ > N 4 :
Consider the sum
vj (x)ϕjt ρ−1 Eρ,ρ(−λjt ρ ).
Since (Â + I) −τ −1 vj(x) = (λj + 1) −τ −1 vj (x), we can rewrite the latter in the form
Therefore, by virtue of Lemma 3.3, one has
.
Using the orthonormality of the system {vj }, we have
For the Mittag-Leffler function with a negative argument we have an estimate (see, for example, [10] , p.13)
Applying this inequality we have
Here we used the inequality Eρ,ρ(−λjt ρ ) < C in the case λj < t −ρ , and the inequality Eρ,ρ(−λjt ρ ) < C 1 λ j t ρ in the case λj > t −ρ . Taking into account (3.10), one can rewrite the estimate (3.9) as
Cϕ. This implies uniform convergence on x ∈Ω of the differentiated sum (3.7) with respect to variables xj, j = 1, . . . , N, for each t ∈ (0, T ]. On the other hand, the sum (3.8) converges for any permutation of its members, as well, since these terms are mutually orthogonal. This implies the absolute convergence of the differentiated sum (3.7) on the same interval t ∈ (0, T ].
Further, it is not hard to see that
Absolute and uniform convergence of the latter series can be proved as above.
Obviously, the function in (3.4) satisfies boundary conditions (2.5) . Considering the initial condition as (see, for example, [10] p. 104)
it is not hard to verify, that this condition is also satisfied.
Hence, if the function ϕ(x) satisfies condition (3.6) , then all the statements of Theorem 3.2 hold. As is shown in work [36] by V.A. Il'in (see also [37] p. 111) the fulfillment of conditions (2.3)-(2.9) guarantee the convergence of the series in (3.6). Thus, Theorem 3.2 is completely proved.
In this section we prove the main result. First we prove the following auxiliary lemma: Proof. Since the system of eigenfunctions {vj (x)} are orthonormal and λ1 = 0, then from (3.4) one has
Let Ψ(ρ) be the logarithmic derivative of the gamma function Γ(ρ) (see, for example, [38] ). Then Γ ′ (ρ) = Γ(ρ)Ψ(ρ), and for ρ ∈ (0, 1) we have Γ(ρ) > 0 and Ψ(ρ) < 0. Therefore,
for t0 > 1 (the case t0 = 1 is obvious). Thus function f (ρ; t0) increases or decreases depending on sign of ϕ1. It is easy to verify equalities (4.1).
Proof of Theorem 2.2. First we show existence of the order of the fractional derivative ρ, which satisfies condition (2.7). We have
It follows from Lemma 4.1 and representation (4.2) immediately, that if 0 < d0 ϕ1 < 1, then there exists a unique ρ, which satisfies condition (2.7).
To prove the uniqueness of a solution of inverse problem (2.4)-(2.7) we suppose that there exist two pairs of solutions {u1, ρ1} and {u2, ρ2} such, that 0 < ρ k < 1, k = 1, 2, and
Consider the following functions
Then, for each j = 1, 2, . . . , we have
Therefore (see (3. 3)), for each j = 1, 2, . . . , w j k (t) = ϕj t ρ k −1 Eρ k ,ρ k (−λjt ρ k ), k = 1, 2, and condition (2.7) implies w 1 1 (t0) = w 1 2 (t0). Since λ1 = 0 we obtain
As we have seen above (see Lemma 4.1), this equation yields ρ1 = ρ2. But in this case w j 1 (t) = w j 2 (t) for all t and j.
Hence
for all j. Since the set of eigenfunctions {vj } is complete in L2(Ω), then we finally have u1(x, t) = u2(x, t). Thus, "if part" of the theorem is proved.
To prove "only if" part of the theorem assume that condition (2.11) is not verified. In this case, as it follows evidently from representation (4.2), equation f (ρ; t0) = d0 has no solution on the interval (0, 1). Hence, in this case the inverse problem does not have a solution. The proof of Theorem 2.2 is complete.
As an example of application of Theorem 2.2 consider the following initial-boundary value problem for one-dimensional diffusion equation (4.6) ∂ ρ t u(x, t) − uxx(x, t) = 0, x ∈ (0, π), t > 0; with the initial condition and the boundary condition (4.8) ux(0, t) = 0, ux(π, t) = 0, t ≥ 0.
where 0 < ρ < 1. In this case the corresponding spectral problem has the set of eigenfunctions {cos kx} complete in L2(0, π), and eigenvalues k 2 , k = 0, 1, · · · . Note that the first eigenvalue in this case is λ1 = 0 and the corresponding eigenfunction is v1(x) = 1. Therefore, condition (2.7) takes the form (4.9) π 0 u(x, t0)dx = d0, t0 ≥ 1. The proof of this theorem immediately follows from Theorem 2.2.