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] ). As the most important applications of this theory one may consider the recent investigations on modeling of COVID-19 outbreak [2] , [3] . The data [9] presented by Johns Hopkins University on the outbreak from different countries seem to show fractional order dynamical processes.
In tern, the well-deserved popularity of the theory attracted the attention of specialists, causing a large number of investigations on mathematical aspects of fractional differential equations and methods to solve them (see, for example, [1] and references therein, [10] - [13] ).
In this paper we will investigate by the Fourier method the solvability (in the classical sense) of initial-boundary value problems for a time-fractional subdiffusion equation with elliptic differential operators of any order, defined on an arbitrary N -dimensional bounded domain Ω with smooth boundary ∂Ω. The fractional part of our equation will be the Riemann -Liouville fractional derivative ∂ ρ t of order 0 < ρ ≤ 1. Many special cases of this problem have been considered by a number of authors using different methods. It has been mainly considered the case of one spatial variable x ∈ R and subdiffusion equation with "elliptical part" uxx (see, for example, handbook [1] , book of A.A. Kilbas et al. [6] and monograph of A. V. Pskhu [11] , and references in these works). In multidimensional case (x ∈ R N ) instead of the differential expression uxx it has been considered either the Laplace operator ( [6] , [14] ), or the elliptic differential or pseudo-differential operator in the whole space R N with constant coefficients [10] . In both cases the authors investigated the Cauchy type problems applying either the Laplace transform or the Fourier transform. In his recent paper [15] A. V. Pskhu considered an initial-boundary value problem for subdiffusion equation with the Laplace operator and domain Ω -a multidimensional rectangular region. The author succeeded to construct the Green's function. Naturally, based on the physical meaning of the problem, it is of interest to consider similar problems for arbitrary bounded spatial domains Ω.
In an arbitrary domain Ω initial-boundary value problems for subdiffusion equations (the fractional part of the equation is a multi-term and initial conditions are non-local) with the Caputo derivatives has been investigated by M. Ruzhansky et al. [16] . The authors proved the existence and uniqueness of the generalized solution to the problem.
Let A(x, D) = |α|≤m aα(x)D α be an arbitrary positive formally selfadjoint (symmetric) elliptic differential operator of order m = 2l with sufficiently smooth coefficients aα(x) in Ω, where α = (α1, α2, ..., αN ) -multi-index and D = (D1, D2, ..., DN ), Dj = ∂ ∂x j . Recall, an operator A(x, D) is elliptic in Ω, if for all x ∈ Ω and ξ ∈ R N one has |α|=m aα(x)ξ α > 0 ξ = 0.
The fractional integration of order ρ < 0 of a function f defined on [0, ∞) in the Riemann -Liouville sense is defined by the formula
provided the right-hand side exists. Here Γ(ρ) 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, [11] , p. 14)
Note if ρ = k, then fractional derivative coincides with the ordinary classical derivation:
Let ρ ∈ (0, 1] be a constant number. Consider the differential equation conditions, where f (x, t), ϕ(x) and coefficients bα,j(x) are given functions. We draw attention to the fact, that in the definition of the regular solution the requirement of continuity in the closed domain of all derivatives included in equation (1.1) is not caused by the merits. However, on the one hand, the uniqueness of just such a solution is proved quite simply, and on the other, the solution found by the Fourier method satisfies the above conditions. We will also call the regular solution simply the solution to the boundary value problem. Application of the Fourier method to problem (1.1) -(1.3) leads us to consider the following spectral problem
S. Agmon [17] found the necessary conditions for boundary ∂Ω of the domain Ω and for coefficients of operators A and Bj , which guarantee compactness of the inverse operator, i.e. the existence of a complete in L2(Ω) system of orthonormal eigenfunctions {v k (x)} and a countable set of positive eigenvalues λ k of spectral problem (1.4) -(1.5). We will call these conditions as condition (A).
In accordance with the Fourier method, we will look for a solution to problem (1.1) -(1.3) in the form of a series:
where functions Tj (t) are solutions of the Cauchy type problem
Here we denoted by fj (t) and ϕj the Fourier coefficients of functions f (x, t) and ϕ(x) with respect to the system of eigenfunctions {v k (x)} correspondingly, defined as a scalar product on L2(Ω), i.e. for example, ϕj = (ϕ, vj ). The unique solution of problem (1.6) has the form (see, for example, [11] , p. 16)
where Eρ,µ is the Mittag-Leffler function of the form
. Throughout what follows, we assume that the condition (A) is satisfied.
To formulate the existence theorem we need to introduce for any real number τ an operator A τ , acting in L2(Ω) in the following waŷ
Obviously, operatorÂ τ with the domain of definition
is selfadjoint. If we denote by A the operator in L2(Ω), acting as Ag(x) = A(x, D)g(x) and with the domain of definition D(A) = {g ∈ C m (Ω) : Bj g(x) = 0, j = 1, ..., l, x ∈ ∂Ω}, then it is not hard to show, that operatorÂ ≡Â 1 is a selfadjoint extension in L2(Ω) of operator A. 
which absolutely and uniformly converges on x ∈Ω and for each t ∈ (0, T ]. 
and the boundary condition is the same as (1.3).
It should be noted that Sh.A. Alimov in his paper [19] presents sufficient conditions for a given function to belong to the domain of definition of operatorÂ τ in terms of various classes of differentiable functions. At the end of this paper we consider initial-boundary value problem (1.1)-(1.3) in torus T N when operator A(x, D) has constant coefficients. We will see that in this case the domain of definition D(Â τ ) coincides with the corresponding Sobolev spaces.
In this section we prove Theorem 1.2. Suppose that all the conditions of this theorem are satisfied and let initial-boundary value problem (1.1) -(1.3) have two regular solutions u1(x, t) and u2(x, t). Our aim is to prove that u(x, t) = u1(x, t) − u2(x, t) ≡ 0. Since the problem is linear, then we have the following homogenous problem for u(x, t):
Let u(x, t) be a regular solution of problem (2.1)-(2.3) and v k be an arbitrary eigenfunction of the problem (1.4)-(1.5) with the corresponding eigenvalue λ k . Consider the function
By definition of the regular solution we may write
or, integrating by parts,
Since condition (1.2) can be rewritten as (see, for example, [11] p. 104)
then, using in (2.4) the homogenous initial condition (2.2) in this form, we have the following Cauchy problem for w k (t):
This problem has the unique solution; therefore, the function defined by (2.4), is identically zero: w k (t) ≡ 0 (see (1.7)). From completeness in L2(Ω) of the system of eigenfunctions {v k (x)}, we have u(x, t) = 0 for all x ∈Ω and t > 0. Hence Theorem 1.2 is proved.
Here we have borrowed some original ideas from the method developed in the work of M.A. Krasnoselski et al. [18] . The following lemma plays a key role in this method ( [18] , p.453).
Then for any |α| ≤ m operator D αÂ−σ (completely) continuously maps from L2(Ω) into C(Ω) and moreover the following estimate holds true
Using this lemma, we prove that one can validly apply the operators D α with |α| ≤ m and ∂ ρ t to the series (1.8) term-by-term. For the Mittag-Leffler function with a negative argument we have an estimate |Eρ,µ(−t)| ≤ C 1+t (see, for example, [11] , p.13). Therefore, since all eigenvalues λj are positive,
where 0 < ε < ρ. Indeed, let t ρ λj < 1, then
and if t ρ λj > 1, then 1 λjt (t ρ λj) ε/ρ > 1 λjt .
Note the series (1.8) is in fact the sum of two series. Consider the first series:
and suppose that function ϕ satisfies the condition of Theorem 1.3, i.e. for some τ > N 2m
We choose a small ε > 0 in such a way, that
Therefore by virtue of Lemma 3.1 one has
Using the orthonormality of the system {vj }, we will have
Therefore we can rewrite the estimate (3.5) as
Cϕ. This implies uniformly on x ∈Ω convergence of the differentiated sum (3.3) with respect to the variables xj for each t ∈ (0, T ]. On the other hand, the sum (3.4) 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.3) on the same interval t ∈ (0, T ]. Now we consider the second part of the series (1.8):
and suppose that function f (x, t) satisfies all the conditions of Theorem 1.3, i.e. the following series converges uniformly on t ∈ [0, T ] for some τ > N 2m :
We have
Then by virtue of Lemma 3.1 one has
Using the orthonormality of the system {vj }, we will have
Now we use estimate (3.2) and apply the generalized Minkowski inequality. Then
Hence, using the same argument as above, we see that the differentiated sum (3.6) with respect to the variables xj converges absolutely and uniformly on (
It is not hard to see that
Absolutely and uniformly convergence of the latter series can be proved as above.
Obviously, function (1.8) satisfies the boundary conditions (1.3). Considering the initial condition as in (2.5), it is not hard to verify, that this condition is also satisfied.
Thus Theorem 1.3 is completely proved. Observe, taking α = 0 in the above estimations, one may obtain an estimate for ||u(x, t)|| C(Ω) , which gives the stability of the solution to problem (1.1)-(1.3).
In the proof of Theorem 1.3 we only use the fact that elliptic operator A has a complete in L2(Ω) set of orthonormal eigenfunctions and Lemma 3.1. This lemma reduces the study of uniform convergence to the study of convergence in L2, where the Parseval equality gives a solution to the problem immediately. Therefore, similar to Theorem 1.3 statement holds true for any operator with these properties. As an example we may consider differential operator with involution and the Bessel operator (see the paper of M. Ruzhansky et al. [16] ), or the initial-boundary value problem in N -dimensional torus:
Let us consider the last case. So let A(D) = |α|=m aαD α be a homogeneous symmetric positive elliptic differential operator with constant coefficients. Let the differential equation and initial condition have the form (0 < ρ ≤ 1)
Instead of boundary conditions (1.3) we consider the 2π-periodic in each argument xj functions and suppose that ϕ(x) and f (x, t) are also 2π-periodic functions in xj. Let A stands for the operator A(D), defined on 2π-periodic functions from C m (R N ). The closureÂ of this operator in L2(T N ) is selfadjoint. It is not hard to see that operatorÂ has a complete (in L2(T N )) set of eigenfunctions {(2π) −N/2 e inx }, and corresponding eigenvalues A(n). Therefore, by virtue of J. von Niemann theorem, for any τ > 0 operatorÂ τ acts as
A τ (n)gne inx , where gn is Fourier coefficients of g ∈ L2(T N ). The domain of definition of this operator is defined from the conditionÂ τ g(x) ∈ L2(T N ) and has the form
In order to define the domain of definition of operatorÂ τ in terms of the Sobolev spaces, we remind the definition of these spaces (see, for example, [20] ): we say that function g ∈ L2(T N ) belongs to the Sobolev space L a 2 (T N ) with the real number a > 0, if the norm
is bounded. When a is not integer, then this space is also called the Liouville space. It is not hard to show the existence of constants c1 and c2 such that
Therefore, D(Â τ ) = L τ m 2 (T N ). Repeating verbatim the proof of Theorem 1.3, the following statement is proved. Remark 4.2. Note, when τ > N 2 , according to the Sobolev embedding theorem, all functions in L τ 2 (T N ) are 2π-periodic continuous functions. The fulfillment of the inverse inequality τ ≤ N 2 , admits the existence of unbounded functions in L τ 2 (T N ) (see, for example, [20] ). Therefore, condition τ > N 2 of this theorem is not only sufficient for the statement to be hold, but it is also necessary.
A 2π-periodical continuous function ϕ(x), defined on T N , is said to be a Hölder continuous with exponent β ∈ (0, 1], if for some constant C, |ϕ(x) − ϕ(y)| ≤ C|x − y| β , for all x, y ∈ T N .
We define C a (T N ) to be the Hölder space of 2π-periodical functions in C [a] (T N ), [a] is the integer part of a, all derivatives D α ϕ, |α| = [a], of which are Hölder continuous with exponent a − [a]. It follows from the embedding theorem C a (T N ) → L a−ε 2 (T N ), ε > 0, that the conditions of the Theorem 4.1 can be formulated in terms of the Hölder spaces, replacing the classes L τ 2 (T N ) by the classes C τ (T N ). Again, condition τ > N 2 is precise: if τ = N 2 , then there exists a function in C τ (T N ), such that the Fourier series of which diverges at some point (see, for example, [20] ).
The authors convey thanks to Sh. A. Alimov and S. R. Umarov for discussions of these results. We wish also to thank B. Turmetov for making us aware of the paper of M. Ruzhansky et al. [16] .