1 Introduction

The time-fractional diffusion equation (tFDE) was derived via the framework of a continuous time random walk under the assumption that the mean waiting time has a power-law decaying tail proportional to t^α, α ∈ (0, 1). The solution of the tFDE accurately describes the power-law decaying behavior in a large number of anomalous diffusion processes [39]. To improve the modeling accuracy of single-term tFDE, diffusion equations with multiple time-derivatives are proposed, which permit to describe also processes whose scaling law changes with time [31, 39]. Such an example is the derived in [43] two time-scale mobile-immobile tFDE model for the subdiffusive transport of solutes in heterogeneous porous media.

The multi-term generalizations of the tFDE, involving time-derivatives of orders distributed in the interval (0, 1], are considered in two forms: the so-called “natural” (or Caputo) form and “modified” (or Riemann-Liouville) form, for the precise definitions see [31, 39] and equations (4.2) and (4.3) below. Initial-boundary value problems for the equation in “natural” form are studied analytically and numerically in a large number of works. Based on eigenfunction expansion, explicit solution is given in [30], the fundamental well-posedness theory is established in [25], maximum principles are derived in [30, 27], smoothing properties are studied in [19], complete monotonicity of the corresponding relaxation functions is discussed in [2, 13]. Regarding numerical treatments, we refer to [19, 44]. For further studies on multi-term equations in Caputo form see e.g. [38], where different types of spatial operators are considered, and [45], where the multi-term time-space Caputo-Riesz fractional diffusion equation on an infinite domain is examined. Concerning the tFDE in “modified” form, explicit expressions of the solutions of two- and three-term equations are derived in terms of series of Fox H-functions in [15, 22, 40]. The Rayleigh-Stokes problem for a generalized second-grade fluid, which is a special case of a two-term “modified”-form equation, is studied analytically and numerically in [7]. Abstract framework for multi-term evolution equations is developed in [23], see also [3].

The multinomial Mittag-Leffler function

where z_j ∈ ℂ, μ_j > 0, β ∈ ℝ, j = 1, …, m, is proposed in [16] and used for solving multi-term fractional differential equations with constant coefficients by operational method in [16, 28, 29]. Originally, it is named “multivariate” Mittag-Leffler function. The current name is attributed to this function in [18]. The multinomial Mittag-Leffler function plays a crucial role in the study of multi-term time-fractional diffusion equations. This is due to the fact that the time-dependent components in the eigenfunction expansion of the solution to initial-boundary value problems for multi-term equations are expressed in terms of multinomial Mittag-Leffler functions, see e.g. [19, 25, 27, 30]. To prove existence and regularity results, estimates for this function are essential.

The following useful estimate is established in [25]: for 0 < β <2, 1 > μ₁ > μ_j > 0, and −K ≤ z_j < 0, where K > 0, j = 2, …, m, and μ₁π/2 < ν < μ₁π, there exists a constant C = C(μ₁, …, μ_m, β, K, ν), such that

Estimate (1.1) is a natural extension of a property of the classical Mittag-Leffler function ([34], Theorem 1.6) and has found numerous applications in the study of initial-boundary value problems for the tFDE in “natural” form, such as in deriving regularity estimates and smoothing properties for the solution [19, 25] and in the study of different types of inverse problems [1, 24, 26, 46]. However, estimate (1.1) is not suitable when the tFDE in “modified” form is considered.

In a series of papers [4, 5, 6] initial-boundary-value problems for diffusion equations with multiple time derivatives and nonlocal boundary conditions are considered. The nonlocal character of one of the boundary conditions leads to a non-selfadjoint problem and multidimensional eigenspaces. This, in turn, implies that the time-dependent components in the generalized eigenfunction expansions of the solutions are expressed in terms of multinomial Mittag-Leffler functions and convolutions of them. It is known that convolution of two classical Mittag-Leffler functions is a Prabhakar function [12]. Therefore, the need of Prabhakar type generalization of the multinomial Mittag-Leffler function naturally emerge in the context of nonlocal boundary value problems for the multi-term tFDE. Such a generalization is defined and used in [4].

For other types of multi-index and multi-variable generalizations of the classical Mittag-Leffler function we refer to the recent surveys [20, 33] and the last edition of the monograph [12].

In this paper the study of the multinomial Prabhakar type function and the multi-term tFDEs in their two forms are intertwined. First, basic properties of the multinomial Prabhakar type function are established with the main emphasis on complete monotonicity. The obtained results are applied to prove well-posedness of the considered equations. As particular examples, the relaxation functions for equations with multiple time-derivatives in the “natural” and “modified” forms are studied in detail and useful estimates are derived. This, in turn, provides some new relations for the multinomial Mittag-Leffler functions. The main tools used in the present work are Laplace transform and Bernstein functions’ technique.

The rest of the paper is organized as follows. Section 2 is concerned with the definition and basic relations for the multinomial Prabhakar type function. Complete monotonicity and asymptotic behavior is studied in Section 3. Multi-term evolution equations are considered in Section 4. Section 5 is devoted to a detailed study of the properties of the related relaxation functions. In Section 6 the moments of the Green function are expressed in terms of multinomial Prabhakar type functions. Definitions and basic properties of fractional calculus operators and Bernstein functions are listed in an Appendix.

2 Multinomial Prabhakar type function: basic relations

The classical Prabhakar function is defined as [11, 35]

where (δ)_k denotes the Pochhammer symbol

and Γ(·) is the Euler gamma function.

A multinomial generalization of the Prabhakar function (2.1) is defined next. For the sake of brevity we use the vector notation μ→=(μ1,μ2,…,μm).

The multinomial Prabhakar function is defined as follows, [4]:

where zj∈ℂ, μj,β,δ∈ℝ, μ_j > 0, j = 1, …, m.

In general, the parameters μ_j, β, δ, are allowed to assume complex values with ℜμ_j > 0. In this work, however, we restrict our attention to real parameters, which are of particular interest for the considered applications.

The classical Prabhakar function (2.1) is recovered from (2.2) for m = 1. The binomial variant (m = 2) of function (2.2) was recently introduced and studied in [10]. In the special case δ = 1 the Pochhammer symbol yields (1)_k = k! and the function (2.2) is the multinomial Mittag-Leffler function

If δ is a negative integer, then the Prabhakar function (2.2) is defined by a finite sum, and Eμ→, β0(⋅)=1/Γ(β).

Let us note that the double summation in (2.2) can be formally replaced by the multiple summation ∑k1,...,km=0∞:=∑k1=0∞…∑km=0∞. This yields a multiple power series, which converges absolutely and locally uniformly, and thus defines an entire function in each z_j, j = 1, …, m. Therefore, both representations are equivalent.

Applying successive term by term differentiation in (2.2) and using the identity (δ)_k+1 = δ(δ + 1)_k, we deduce the relation

which generalizes a well-known identity for m = 1, see e.g. [34], Eq.(2.1).

In the rest of this work we are concerned only with the following multinomial Prabhakar type function of a single variable t > 0, which is of particular importance for the study of multi-term time-fractional equations

where μ_j > 0, β > 0, δ ∈ ℝ, a_j > 0, j = 1, …, m. For the sake of brevity the short notation ℰμ→,βδ(t;a→) is used for the function (2.4). Definition (2.2) yields the series representation

The first terms in the power series (2.5) give the following asymptotic expansion for t → 0:

We study the multinomial Prabhakar type function (2.4) applying Laplace transform (LT) technique. For this reason we are concerned only with locally integrable functions ℰμ→,βδ(t;a→). Taking into account (2.6), this is guaranteed by the assumptions on the parameters of function (2.4).

The LT pair (2.7) shows that, in general, the representation as a multinomial Prabhakar type function is not unique. For example, the identity ℰμ,β2δ(t;a)=ℰ(μ,2μ),βδ(t;2a,a2) can be proven by the use of (2.7). Moreover, the order of parameters μ_j (resp. a_j) in (2.4) can be changed simultaneously. For clarity, in what follows we choose the representation with minimal m and when a special arrangement of the parameters μ_j (resp. a_j) is assumed, this is explicitly stated.

A reduction of parameters result is established next.

In the particular case δ = 1, m = 2, representation (2.9) appears in the context of two-term time-fractional equations, where instead of the usual in this case binomial Mittag-Leffler function, an infinite series of Prabhakar functions is used, see e.g. [5, 6]. Let us note that further reduction identities can be found in [34], Chapter 5.

The integration, differentiation and convolution properties for the multinomial Prabhakar type functions, given next, extend those for the classical Prabhakar function (see e.g. [11]).

The above identities can be verified directly from the series definition (2.5), or, proving by the use of (2.7) that the Laplace transforms of both sides coincide. Technically, the second method is easier and the proofs are straightforward. For this reason they are omitted here. In the binomial case m = 2 the above identities are proved in detail in [10].

3 Complete monotonicity

This section is devoted to the study of complete monotonicity property of the multinomial Prabhakar type function (2.4) for t > 0. Concerning the classical Prabhakar type function the current most general result states that the function tβ−1Eμ,βδ(−tμ), t > 0, is completely monotone if the parameters satisfy the conditions [11]

A detailed proof can be found in [8]. This result is extended next to the multinomial case. The proof uses the Bernstein functions’ technique, for details see Appendix. It is based on the following result.

Let us note that the condition β ≤ 1 is also necessary for complete monotonicity property (3.2). Indeed, ℰμ→,βδ(t;a→)∈Cℳℱ implies that the asymptotic expansions of this function for t → 0 as well as for t → +∞ should be completely monotone functions. We see from (2.6) that at t → 0 the function ℰμ→,βδ(t;a→) behaves as t^β−1/Γ(β), which is completely monotone only when β ≤ 1.

Next we derive the asymptotic expansion for t → ∞. To this end we need the expansion of ℰ^μ→,βδ(s;a→) for s → 0. Let μ₁ > μ₂ > … > μ_m > 0. Then for s → 0,

and, therefore

From the asymptotic behavior of the Prabhakar function (see e.g. [11], eq. (3.13)) the leading term as t → +∞ is obtained as follows

We observe that the leading terms in (3.4) are completely monotone functions under the assumptions of Theorem 3.2.

Let us point out that (3.4) can be guaranteed only when a_j > 0 for each j = 1, …, m. In the classical case m = 1 this is known [11]. A relevant counterexample concerning the two-term case is provided in [25], Remark 4.1.

We also note that, according to (P3) and (P9), the complete monotonicity property (3.2) implies that ℰ^μ→,βδ(s;a→) can be analytically extended to the whole complex plane cut along the negative real axis. Therefore, the function sμ1+amsμ1−μm+…+a2sμ1−μ2+a1 should not have any zeros there. This is guaranteed by the assumptions μ_j < μ₁ ≤ 1 and a_j > 0. The question whether these conditions are also necessary for complete monotonicity property (3.2) in the multinomial case needs further investigation. Let us mention here the study in [41] of the zeros of such functions with μ₁ > 1, which could be helpful.

Further, let us note that identity (2.7) implies

Therefore, according to property (P5) in the Appendix, ℰ^μ→,βδ(s;a→)∈Sℱ if and only if ℰ^μ→,1−β−δ(s;a→)∈Sℱ. If β ∈ (0, 1) then both Laplace transforms vanish as s → ∞ and according to (P3) ℰμ→,βδ(t;a→)∈Cℳℱ if and only if ℰμ→,1−β−δ(t;a→)∈Cℳℱ. In other words, identity (3.5) implies that ℰμ→,βδ(t;a→) and ℰμ→,1−β−δ(t;a→) are Sonine kernels, that is

and the complete monotonicity of the one implies the complete monotonicity of the other. In this way we obtained the following corollary.

4 Multi-term time-fractional evolution equations

Let DCtα and Dtα denote the fractional time-derivatives in the Caputo and Riemann-Liouville sense, respectively, and let A be a generator of a bounded analytic semigroup (see e.g. [9]). In this section we are concerned with the two types of multi-term generalizations of the single-term fractional evolution equation (t > 0, x ∈ ℝⁿ)

Let 1 ≥ α > α₁ > ... > α_m > 0, b_j > 0, j = 1, ..., m. We consider the multi-term tFDE in the so-called “natural” (or Caputo) form

and in the so-called “modified” (or Riemann-Liouville) form

Let us point out that in our considerations of equations (4.2) and (4.3) the case α = 1 is included in order to cover important models, such as the two time-scale mobile-immobile model for the subdiffusive transport of solutes in heterogeneous porous media [43], and the Rayleigh-Stokes problem for a generalized second grade fluid [7]. Therefore, it is not possible to use for the study of equations (4.2) and (4.3) the framework of general fractional derivative proposed in [21]. Indeed, if for example, the multi-term derivative operator in (4.2) with α = 1 is represented as a general fractional derivative, the corresponding kernel of this derivative would contain a Dirac delta function, see also [17] for a related discussion.

For a unified approach to the two types of multi-term tFDEs, (4.2) and (4.3), we rewrite them for f ≡ 0 as a Volterra integral equation

where the kernel κ(t) = κ₁(t) in the case of equation (4.2) and κ(t) = κ₂(t) in the case of equation (4.3). The Laplace transforms of the kernels obey κi^(s)=1/gi(s), i = 1, 2, where

Therefore, taking into account (2.7), we deduce

The kernels κi(t)∈C(ℝ+)∩Lloc1(ℝ+) are completely monotone functions, see Theorem 3.2. Then [36], Corollary 2.4, implies that the problem (4.4) is well-posed and admits a bounded analytic solution operator S(t).

We observe that the functions g₁(s) and g₂(s) in (4.5) are complete Bernstein functions (the second one as inverse of a Stieltjes function, see (P4)). Moreover, according to Proposition 3.1, a stronger property is satisfied: gi(s)1/α∈Cℬℱ, i = 1, 2. This together with (P6) also implies

Due to property (4.8) the following subordination result can be established in the same way as Theorem 5.1 in [3].

Let us note that subordination results in a more general setting are established in [37] by constructing an appropriate transmutation operator.

5 Relaxation functions

Setting A = −λ, λ > 0, in equations (4.2) and (4.3) leads to two forms of multi-term relaxation equations. In this section we study the properties of the relaxation functions, obtained as solutions of these equations.

By the use of Laplace transform we deduce that the solution of the relaxation equation in “natural” form

is given by

and the solution of the relaxation equation in “modified” form

is represented as

where the functions u₁(t; λ), v₁(t; λ), and u₂(t; λ) satisfy the following Laplace transform identities

with g₁(s) and g₂(s) defined in (4.5).

The functions u₁(t; λ) and v₁(t; λ) are the relaxation functions related to problem (5.1) and u₂(t; λ) is the relaxation function related to problem (5.3). Laplace transform inversion in (5.5) by the use of (2.7) yields the following explicit representations in terms of multinomial Mittag-Leffler functions

In the single term case the relaxation functions reduce to the classical Mittag-Leffler functions ui(t;λ)=Eα(−λtα), i = 1, 2, and v1(t;λ)=tα−1Eα,α(−λtα).

Subordination identities for the relaxation functions u_i(t; λ) can be derived from the scalar version of Theorem 4.1, where S(t) = u_i(t; λ), Sβ(t)=Eβ(−λtβ). In particular, for β = 1 it follows

where the functions φ_i(t, τ) are nonnegative and normalized. A subordination result for the third relaxation function v₁(t; λ) is given next.

By fractional integration of (5.10) and taking into account (5.8) and identity (2.10) we deduce the following representation for completely monotone multinomial Mittag-Leffler functions, which is of independent interest.

Some properties of the relaxation functions for the equations in “natural” and “modified” forms, including useful estimates, are collected in the next theorem.

Let us note that all statements in Theorem 5.2 extend known properties of the classical Mittag-Leffler functions Eα(−λtα) and tα−1Eα,α(−λtα).

6 Moments of the Green functions

As an application of the multinomial Prabhakar type functions (2.5), in this section we derive expressions for the moments of the Green functions of the multiterm tFDEs in terms of such functions. Consider the Cauchy problem for the multi-term tFDEs (4.2) and (4.3), where A=(∂∂x)2, x ∈ ℝ. The fundamental solution (Green function) G(x,t) is defined by assuming the initial and boundary conditions

where δ(·) is the Dirac delta function. Applying as usual Laplace transform with respect to the temporal variable and Fourier transform with respect to the spatial variable, we derive for the Green function G(x,t) in Fourier-Laplace domain (see e.g. [31, 39])

Here g(s) = g₁(s) in the case of “natural”-form equation (4.2) and g(s) = g₂(s) in the case of “modified”-form equation (4.3), and the definitions of these functions are given in (4.5). By Fourier inversion in (6.1) the Laplace transform of the solution is obtained as follows

Let γ > 0. Representation (6.2) implies for the Laplace transforms of the moments 〈|x|^γ(t)〉 of order γ

where the formula ∫0∞xb−1e−ax dx=Γ(b)a−b is used. Taking inverse Laplace transform, we obtain by the use of (2.7)

for the “natural”-form diffusion equation (4.2), where C₁ = Γ(γ + 1), and

for the “modified”-form equation (4.3), where C2=Γ(γ+1)bmγ/2.

Let us note that the indices in the brackets of the above multinomial Prabhakar type functions are specially arranged, so that the first index, α − α_m, is the largest. The obtained representations for the moments, together with the properties (2.10), (3.2), and (3.6), imply that the moments of the Green functions of both equations are Bernstein functions (integrals of completely monotone functions) provided αγ ≤ 2.

The corresponding mean squared displacements 〈|x|i2(t)〉 are derived by setting γ = 2. This yields

where the functions l₁(t) and l₂(t) are defined in (5.18) and (5.19). As we see, l₂(t) is a finite sum, and this is the case for all moments of even order for the equation in “modified” form.

The asymptotic behavior of the derived moments can be deduced from the asymptotic expansions (2.6) and (3.4) for the multinomial Prabhakar type functions. In this way we obtain 〈|x|1γ(t)〉~tαγ/2 as t → 0 and 〈|x|1γ(t)〉~tαmγ/2 as t → ∞ for the equation (4.2) in “natural” form, while for the equation (4.3) in “modified” form the opposite behavior is observed: 〈|x|2γ(t)〉~tαmγ/2 as t → 0 and 〈|x|2γ(t)〉~tαγ/2 as t → ∞. The established results are in agreement with those given in [31, 39] for particular cases.