On the Time-Fractional Hilfer Telegraph Equation with Variable Coefficients

1. Introduction

Fractional differential equations play a central role in the modelling of processes with memory and hereditary effects, and arise in areas such as anomalous diffusion, viscoelasticity, and wave propagation in complex media. Among these models, fractional telegraph equations are of particular interest, as they interpolate between diffusion and wave phenomena and provide a flexible framework for describing damped propagation and transport processes [9,11,12,14,24].

A significant challenge arises when the governing equations involve variable coefficients. From a modelling perspective, time-dependent coefficients allow the description of non-homogeneous media and processes with evolving properties. From an analytical standpoint, however, the presence of variable coefficients destroys the simple algebraic structure available in constant-coefficient problems, leading to non-trivial operator compositions and making explicit solution formulas considerably more difficult to obtain. This has motivated a substantial literature on fractional differential equations with variable coefficients. One of the early studies in this direction is due to Dzhrbashyan and Nersessyan, who investigated the Cauchy problem for fractional-order differential equations and helped lay the foundation for later developments in fractional calculus with non-constant coefficients [2]. Subsequent work by Atanackovic and Stankovic further clarified the structure of linear fractional differential equations with variable coefficients [1]. These contributions show that the variable-coefficient setting is not only natural from the modelling point of view, but also mathematically rich and technically delicate.

In recent years, substantial progress has been made in the analysis of fractional differential equations with variable coefficients. Several works have focused on obtaining explicit representations. For example, Kim and O studied explicit representations of Green’s functions for linear fractional differential operators with variable coefficients [16]. Pak et al. derived analytical solutions for linear inhomogeneous fractional differential equations with continuous variable coefficients [17]. More recently, Restrepo and his collaborators obtained explicit formulas for linear differential equations with variable coefficients in terms of infinite series of fractional integro-differential operators and also studied Prabhakar-type linear differential equations with variable coefficients [4,5,6,7,19]. In a closely related direction, Restrepo, Ruzhansky, and Suragan developed explicit solutions for linear variable-coefficient fractional differential equations with respect to functions [18].

Nevertheless, for telegraph-type equations involving Hilfer fractional derivatives, explicit solution formulas in the presence of variable coefficients remain relatively scarce. Hilfer fractional derivatives provide a natural and flexible framework for such problems. By interpolating between the Riemann–Liouville and Caputo derivatives, they allow the treatment of fractional initial conditions and provide a unified setting for a wide class of fractional evolution equations. When combined with operator-based techniques, they offer an effective way to analyse linear problems with variable coefficients.

In this work, we investigate a time-fractional telegraph equation with time-dependent coefficients within the Hilfer framework. Our main objective is to derive explicit representations for the solution of the associated Cauchy problem and to characterize the precise influence of these variable coefficients on the solution structure. To achieve this, we reformulate the problem in the Fourier domain and apply the framework established by Restrepo and Surugan [19] to express the solution as a convergent Neumann-type series, where the effect of the time-dependent coefficients is encoded in iterated operator compositions.

A key feature of the method is that it accommodates both commutative and non-commutative operator structures. In the general case, the lack of commutativity leads naturally to expansions involving ordered compositions. In contrast, when the coefficients are constant, the operators commute and the solution series can be summed explicitly, yielding closed-form expressions in the Fourier and space-time domains. We also treat the case of power-law coefficients, which highlights the full complexity of the non-commutative setting.

The paper is organized as follows. In Section 2, we recall the necessary preliminaries from fractional calculus, integral transforms, and Mittag-Leffler functions, alongside establishing key auxiliary results used throughout the text. In Section 3, we formulate the time-fractional telegraph equation with time-dependent coefficients and derive the explicit solution to the corresponding Cauchy problem in the Fourier domain. We then analyze the specific cases of constant and power-law coefficients, obtaining their respective closed-form representations. The general space-time solution is subsequently established by taking the inverse Fourier transform, interpreted in a distributional sense. In the constant-coefficient regime, our Fourier-domain and space-time representation formulas recover established results from the literature, validating the consistency of the approach. Finally, in Section 4, we present concluding remarks and discuss a natural extension of this model to the $\psi$-Hilfer framework.

2. Preliminaries

2.1. Fractional Derivatives and Integral Transforms

In this section, we recall several definitions and properties of fractional operators and special functions required for our main results.

Definition 1.

Let $a,b\in \mathbb{R}$ with $a<b$, and let $\alpha >0$. The left Riemann-Liouville (RL) fractional integral $I_{t,a^{+}}^{\alpha }$ of order $\alpha >0$ of a function $f\in L^1\left([a,b]\right)$ is defined by (cf. [15])

$$\begin{array}{cc}\hfill & I_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_a^t{(t-w)}^{\alpha -1}f\left(w\right)dw,t>a.\hfill \end{array}$$

The following lemma establishes the asymptotic behavior of the Riemann-Liouville fractional integral near the lower integration limit for continuous functions, which will be needed later.

Lemma 1.

Let $a,b\in \mathbb{R}$ with $a<b$, let $\alpha >0$, and let $f\in C\left(\right[a,b\left]\right)$. Then, as $t\to a^{+}$,

$$\begin{array}{c}\hfill I_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{f\left(a\right)}{\Gamma \left(\alpha +1\right)}}{(t-a)}^{\alpha }+o\left({(t-a)}^{\alpha }\right).\end{array}$$

Proof. 

By the change of variable $u=t-w$,

$$\begin{array}{c}\hfill I_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^{t-a}{u}^{\alpha -1}f(t-u)du.\end{array}$$

Since f is continuous at a, $f(t-u)=f\left(a\right)+o\left(1\right)$, uniformly for $u\in [0,t-a]$ as $t\to a^{+}$. Therefore,

$$\begin{array}{c}\hfill I_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^{t-a}{u}^{\alpha -1}\left(f\left(a\right)+o\left(1\right)\right)du.\end{array}$$

As the term $o\left(1\right)$ is uniform, it can be taken outside the integral, yielding

$$\begin{array}{c}\hfill I_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{f\left(a\right)+o\left(1\right)}{\Gamma \left(\alpha \right)}}{\int }_0^{t-a}{u}^{\alpha -1}du={\displaystyle \frac{f\left(a\right)+o\left(1\right)}{\Gamma \left(1+\alpha \right)}}{(t-a)}^{\alpha },\end{array}$$

which is equivalent to

$$\begin{array}{c}\hfill I_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{f\left(a\right)}{\Gamma \left(1+\alpha \right)}}{(t-a)}^{\alpha }+o\left({(t-a)}^{\alpha }\right).\end{array}$$

□

Remark 1.

Under the assumptions of Lemma 1, we have

$$\begin{array}{c}\hfill \underset{t\to a^{+}}{lim}{I}_{t,a^{+}}^{\alpha }f\left(t\right)=0.\end{array}$$

The same conclusion holds for every $f\in L^{\infty }\left([a,b]\right)$, since

$$\begin{array}{c}\hfill \left|I_{t,a^{+}}^{\alpha }f\left(t\right)\right|\le {\displaystyle \frac{{\parallel f\parallel }_{\infty }}{\Gamma \left(1+\alpha \right)}}{(t-a)}^{\alpha }.\end{array}$$

Definition 2.

The left Hilfer (or composite) fractional derivative ${{}^H\partial }_{a^{+}}^{\alpha ,\mu }$ of order $\alpha >0$ and type $0\le \mu \le 1$ of a function f is given by (cf. [11,12,13,25])

$$\begin{array}{cc}\hfill & {{}^H\partial }_{t,a^{+}}^{\alpha ,\mu }f\left(t\right)=I_{t,a^{+}}^{\mu (m-\alpha )}\left({\displaystyle \frac{d^m}{d{t}^m}}\left(I_{t,a^{+}}^{(1-\mu )(m-\alpha )}f\left(t\right)\right)\right),\hfill \end{array}$$

(1)

where $m=\lfloor \alpha \rfloor +1$ and $\lfloor ·\rfloor$ denotes the floor function. Here, $f\in L^1\left([a,b]\right)$ and $I_{t,a^{+}}^{(1-\mu )(m-\alpha )}f\in A{C}^m\left([a,b]\right)$.

Note that this formulation interpolates between classical fractional derivatives:

• For $\mu =0$, it reduces to the left Riemann-Liouville fractional derivative of order $\alpha$, given by

  $$\begin{array}{cc}\hfill & {{}^{RL}\partial }_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\Gamma \left(m-\alpha \right)}}{\displaystyle \frac{d^m}{d{t}^m}}{\int }_a^t{\displaystyle \frac{f\left(w\right)}{{(t-w)}^{1-m+\alpha }}}dw,\hfill \end{array}$$

  where $f\in L^1\left([a,b]\right)$.
• For $\mu =1$ we recover the left Caputo fractional derivative of order $\alpha$, given by

  $$\begin{array}{cc}\hfill & {{}^C\partial }_{t,a^{+}}^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\Gamma \left(m-\alpha \right)}}{\int }_a^t{\displaystyle \frac{f^{\left(m\right)}\left(w\right)}{{(t-w)}^{1-m+\alpha }}}dw,\hfill \end{array}$$

  where $f\in A{C}^m\left([a,b]\right)$.

We shall make use of the following composition rules:

$$\begin{array}{c}\hfill {{}^H\partial }_{t,a^{+}}^{\gamma ,\mu }{I}_{t,a^{+}}^{\gamma }f\left(t\right)=f\left(t\right)\mathrm{and}{{}^H\partial }_{t,a^{+}}^{{\gamma }_1,\mu }{I}_{t,a^{+}}^{{\gamma }_2}f\left(t\right)=I_{t,a^{+}}^{{\gamma }_2-{\gamma }_1}f\left(t\right).\end{array}$$

The proof of the second composition rule follows the same argument used in the proof of equation (3.3) in [19]. The previous definitions can be naturally extended to functions of several variables by employing partial fractional integrals and derivatives (see [20]).

For the power function, we have the following integro-differentiation rules:

$$\begin{array}{cc}\hfill & I_{t,a^{+}}^{\gamma }\left[{(t-a)}^{\lambda -1}\right]={\displaystyle \frac{\Gamma \left(\lambda \right)}{\Gamma \left(\lambda +\gamma \right)}}{(t-a)}^{\lambda -1+\gamma },\hfill \\ \hfill & {{}^H\partial }_{t,a^{+}}^{\gamma ,\mu }\left[{(t-a)}^{\lambda -1}\right]\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & =\left\{\begin{array}{cc}{\displaystyle \frac{\Gamma \left(\lambda \right)}{\Gamma \left(\lambda -\gamma \right)}{(t-a)}^{\lambda -1-\gamma },}\hfill & \lambda -1+(1-\mu )(m-\gamma )>\lfloor \gamma \rfloor \hfill \\ {\displaystyle 0,}\hfill & \lambda -1+(1-\mu )(m-\gamma )=p,p=0,\dots ,m-1\hfill \end{array}\right..\hfill \end{array}$$

(3)

In this work, we make use of the n-dimensional Fourier transform and the Laplace transform. For a real-valued Lebesgue integrable function f on ${\mathbb{R}}^n$, the n-dimensional Fourier transform is defined by (see [15])

$$\begin{array}{cc}\hfill & \mathcal{F}\left\{f\left(x\right)\right\}\left(\xi \right)=\widehat{f}\left(\xi \right)={\int }_{{\mathbb{R}}^n}{e}^{i\xi ·x}f\left(x\right)dx,x,\xi \in {\mathbb{R}}^n,\hfill \end{array}$$

where $\xi ·x$ represents the dot product of vectors $\xi$ and x, while the corresponding inverse Fourier transform is given by

$$\begin{array}{cc}\hfill & f\left(x\right)={\mathcal{F}}^{-1}\left\{\widehat{f}\left(\xi \right)\right\}\left(x\right)={\displaystyle \frac{1}{{\left(2\pi \right)}^n}}{\int }_{{\mathbb{R}}^n}{e}^{-i\xi ·x}\widehat{f}\left(\xi \right)d\xi .\hfill \end{array}$$

We will also use the following well-known Convolution Theorem:

$$\begin{array}{c}\hfill \mathcal{F}\left\{\left(f{*}_{x}g\right)\left(x\right)\right\}\left(\xi \right)=\mathcal{F}\left\{f\right\}\left(\xi \right)\mathcal{F}\left\{g\right\}\left(\xi \right),\end{array}$$

(4)

where the convolution ${*}_x$ is given by

$$\begin{array}{c}\hfill \left(f{*}_{x}g\right)\left(x\right)={\int }_{{\mathbb{R}}^n}f(x-z)g\left(z\right)dz.\end{array}$$

(5)

The n-dimensional Laplace operator ${\Delta }_x={\sum }_{i=1}^n{\displaystyle \frac{{\partial }^2}{\partial x_i^2}}$ has Fourier symbol $-{\left|\xi \right|}^2$, i.e.:

$$\begin{array}{c}\hfill \mathcal{F}\left\{{\Delta }_{x}f\left(x\right)\right\}\left(\xi \right)=-{\left|\xi \right|}^2\mathcal{F}\left\{f\left(x\right)\right\}\left(\xi \right),\end{array}$$

(6)

valid for sufficiently regular functions f (for instance, $f\in \mathcal{S}\left({\mathbb{R}}^n\right)$ or $f\in L^2\left({\mathbb{R}}^n\right)$ with ${\Delta }_{x}f\in L^2\left({\mathbb{R}}^n\right)$). Since convolutions in the time variable will also appear in the representation of the solution, we recall that

$$\begin{array}{c}\hfill \left(f{*}_{t}g\right)\left(t\right)={\int }_0^{t}f(t-w)g\left(w\right)dw.\end{array}$$

(7)

2.2. Special Functions

The multivariate Mittag-Leffler (ML) function $E_{(a_1,\dots ,a_n),\lambda }\left(z_1,\dots ,z_n\right)$ of n complex variables $z_1,\dots ,z_n\in \mathbb{C}$ with complex parameters $a_1,\dots ,a_n,\lambda \in \mathbb{C}$ (with positive real parts) is defined by (see [21]):

$$\begin{array}{c}\hfill {\displaystyle E_{(a_1,\dots ,a_n),\lambda }\left(z_1,\dots ,z_n\right)={\displaystyle \sum _{k=0}^{\infty }}\sum _{\begin{array}{c}{l}_1+\dots +l_n=k\\ l_1,\dots ,l_n\ge 0\end{array}}\left(\genfrac{}{}{0pt}{}{k}{l_1,\dots ,l_n}\right)\frac{{\prod }_{i=1}^n{z}_i^{l_i}}{\Gamma \left(\lambda +{\sum }_{i=1}^n{a}_i{l}_i\right)},}\end{array}$$

where the multinomial coefficients are given by

$$\begin{array}{c}\hfill \left({\displaystyle \genfrac{}{}{0pt}{}{k}{l_1,\dots ,l_n}}\right):={\displaystyle \frac{k!}{l_1!\times \dots \times l_n!}}.\end{array}$$

For $n=2$ and $n=3$, we obtain the bivariate and trivariate ML functions, which can be expressed, respectively, as

$$\begin{array}{cc}\hfill & E_{(a_1,a_2),\lambda }\left(z_1,z_2\right)={\displaystyle \sum _{l_1=0}^{\infty }}{\displaystyle \sum _{l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_2\right)}},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & E_{(a_1,a_2,a_3),\lambda }\left(z_1,z_2,z_3\right)={\displaystyle \sum _{l_1=0}^{\infty }}{\displaystyle \sum _{l_2=0}^{\infty }}{\displaystyle \sum _{l_3=0}^{\infty }}{\displaystyle \frac{(l_1+l_2+l_3)!}{l_1!l_2!l_3!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}{z}_3^{l_3}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_2+a_3{l}_3\right)}}.\hfill \end{array}$$

(9)

The following addition formula for the bivariate ML function was stated in [9].

Lemma 2

(Addition formula). Let $a_1,a_2,\lambda ,z_1,z_2\in \mathbb{C}$, with $\mathrm{Re}\left(a_1\right),\mathrm{Re}\left(a_2\right)>0$. Then

$$\begin{array}{c}\hfill E_{(a_1,a_2),\lambda }\left(z_1,z_2\right)={\displaystyle \frac{1}{\Gamma \left(\lambda \right)}}+z_1{E}_{(a_1,a_2),\lambda +a_1}\left(z_1,z_2\right)+z_2{E}_{(a_1,a_2),\lambda +a_2}\left(z_1,z_2\right).\end{array}$$

Since this lemma was formulated in [9] without a formal proof, we provide one here.

Proof. 

The proof is based on the series representation of the ML function. We separate the term corresponding to $(l_1,l_2)=(0,0)$ from the remaining terms and use the identity

$$\begin{array}{c}\hfill {\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}={\displaystyle \frac{(l_1+l_2-1)!}{(l_1-1)!l_2!}}+{\displaystyle \frac{(l_1+l_2-1)!}{l_1!(l_2-1)!}},\mathrm{for}{l}_1,l_2\ge 1.\end{array}$$

(10)

Hence,

$$\begin{array}{cc}\hfill E_{(a_1,a_2),\lambda }\left(z_1,z_2\right)& ={\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_2\right)}}\hfill \\ \hfill & ={\displaystyle \frac{1}{\Gamma \left(\lambda \right)}}+\sum _{\begin{array}{c}{l}_1,l_2\ge 0\\ (l_1,l_2)\ne (0,0)\end{array}}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_2\right)}}\hfill \\ \hfill & ={\displaystyle \frac{1}{\Gamma \left(\lambda \right)}}+\sum _{l_1\ge 1,l_2\ge 0}{\displaystyle \frac{(l_1+l_2-1)!}{(l_1-1)!l_2!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_2\right)}}\hfill \\ \hfill & +\sum _{l_1\ge 0,l_2\ge 1}{\displaystyle \frac{(l_1+l_2-1)!}{l_1!(l_2-1)!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_2\right)}}.\hfill \end{array}$$

Now, performing the index shifts $l_1-1\mapsto l_1$ in the first series and $l_2-1\mapsto l_2$ in the second series, we obtain

$$\begin{array}{cc}\hfill E_{(a_1,a_2),\lambda }\left(z_1,z_2\right)& ={\displaystyle \frac{1}{\Gamma \left(\lambda \right)}}+\sum _{l_1,l_2\ge 0}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{z_1^{l_1+1}{z}_2^{l_2}}{\Gamma \left(\lambda +a_1+a_1{l}_1+a_2{l}_2\right)}}\hfill \\ \hfill & +\sum _{l_1,l_2\ge 0}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2+1}}{\Gamma \left(\lambda +a_2+a_1{l}_1+a_2{l}_2\right)}}\hfill \\ \hfill & ={\displaystyle \frac{1}{\Gamma \left(\lambda \right)}}+z_1{E}_{(a_1,a_2),\lambda +a_1}\left(z_1,z_2\right)+z_2{E}_{(a_1,a_2),\lambda +a_2}\left(z_1,z_2\right).\hfill \end{array}$$

□

The following lemma establishes a reduction property relating the trivariate and bivariate ML functions.

Lemma 3

(Reduction property). Let $a_1,a_2,\lambda ,z_1,z_2,z_3\in \mathbb{C}$, with $\mathrm{Re}\left(a_1\right),\mathrm{Re}\left(a_2\right)>0$. Then

$$\begin{array}{c}\hfill E_{(a_1,a_2,a_2),\lambda }\left(z_1,z_2,z_3\right)=E_{(a_1,a_2),\lambda }\left(z_1,z_2+z_3\right).\end{array}$$

(11)

Proof. 

Setting $a_3=a_2$ into (), we have

$$\begin{array}{c}\hfill E_{(a_1,a_2,a_2),\lambda }\left(z_1,z_2,z_3\right)={\displaystyle \sum _{l_1=0}^{\infty }}{\displaystyle \sum _{l_2=0}^{\infty }}{\displaystyle \sum _{l_3=0}^{\infty }}{\displaystyle \frac{(l_1+l_2+l_3)!}{l_1!l_2!l_3!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}{z}_3^{l_3}}{\Gamma \left(\lambda +a_1{l}_1+a_2(l_2+l_3)\right)}}.\end{array}$$

Making the change of variables $l_4=l_2+l_3$ which means that $l_3=l_4-l_2$, we can write

$$\begin{array}{c}\hfill E_{(a_1,a_2,a_2),\lambda }\left(z_1,z_2,z_3\right)={\displaystyle \sum _{l_1=0}^{\infty }}{\displaystyle \sum _{l_4=0}^{\infty }}{\displaystyle \sum _{l_2=0}^{l_4}}{\displaystyle \frac{(l_1+l_4)!}{l_1!l_2!(l_4-l_2)!}}{\displaystyle \frac{z_1^{l_1}{z}_2^{l_2}{z}_3^{l_4-l_2}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_4\right)}}.\end{array}$$

By the binomial formula, we have

$$\begin{array}{c}\hfill {\displaystyle \sum _{l_2=0}^{l_4}}{\displaystyle \frac{z_2^{l_2}{z}_3^{l_4-l_2}}{l_2!(l_4-l_2)!}}={\displaystyle \frac{{(z_2+z_3)}^{l_4}}{l_4!}}.\end{array}$$

Therefore, we get

$$\begin{array}{cc}\hfill E_{(a_1,a_2,a_2),\lambda }\left(z_1,z_2,z_3\right)& ={\displaystyle \sum _{l_1=0}^{\infty }}{\displaystyle \sum _{l_4=0}^{\infty }}{\displaystyle \frac{(l_1+l_4)!}{l_1!l_4!}}{\displaystyle \frac{z_1^{l_1}{(z_2+z_3)}^{l_4}}{\Gamma \left(\lambda +a_1{l}_1+a_2{l}_4\right)}}\hfill \\ \hfill & =E_{(a_1,a_2),\lambda }\left(z_1,z_2+z_3\right).\hfill \end{array}$$

□

3. Time-Fractional Telegraph Equation with Hilfer Derivatives and Variable Coefficients

We consider the following time-fractional telegraph equation with time-dependent coefficients and Hilfer fractional derivatives:

$$\begin{array}{cc}\hfill & {{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}u(x,t)+a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}u(x,t)+a_2\left(t\right){\Delta }_{x}u(x,t)+a_3\left(t\right)u(x,t)=q(x,t),\hfill \end{array}$$

(12)

subject to the initial conditions:

$$\begin{array}{cc}\hfill & \left(I_{t,0^{+}}^{(1-{\mu }_2)(2-{\alpha }_2)}u\right)(x,0^{+})=g_1\left(x\right)\mathrm{and}{\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{(1-{\mu }_2)(2-{\alpha }_2)}u\right)(x,0^{+})=g_2\left(x\right),\hfill \end{array}$$

(13)

where the conditions in (13) are interpreted in the limit as $t\to 0^{+}$. Here, $(x,t)\in {\mathbb{R}}^n\times [0,T]$ and ${\Delta }_x$ denotes the Laplacian with respect to the spatial variable in ${\mathbb{R}}^n$. Moreover, the fractional orders of differentiation satisfy ${\alpha }_1\in (0,1]$ and ${\alpha }_2\in (1,2]$, while the types of the derivatives are ${\mu }_1,{\mu }_2\in [0,1]$. We further assume that the source term q, the initial data $g_1$ and $g_2$, and the variable coefficients $a_i$ satisfy the following regularity conditions:

$$\begin{array}{c}\hfill q\in C([0,T];L^1\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right)),g_1,g_2\in L^1\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right),a_1,a_2,a_3\in C\left([0,T]\right).\end{array}$$

(14)

Here, the time-dependent map $t\mapsto q(·,t)$ belongs to the space of continuous functions on $[0,T]$ taking values in the intersection space $L^1\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right)$. These assumptions allow us to apply the Fourier transform with respect to the spatial variable and establish the operational framework developed below. Although the analysis can be carried for general continuous coefficients, the physical telegraph-type regime corresponds to the sign conditions

$$\begin{array}{c}\hfill a_1\left(t\right)>0,a_2\left(t\right)<0,a_3\left(t\right)\ge 0,\forall t\in [0,T].\end{array}$$

We seek classical solutions $u:{\mathbb{R}}^n\times (0,T]\to \mathbb{R}$ such that, for each $t\in (0,T]$, the spatial profile satisfies $u(·,t)\in C^2\left({\mathbb{R}}^n\right)\cap L^1\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right)$, subject to the far-field decay condition

$$\begin{array}{c}\hfill \underset{\left|x\right|\to \infty }{lim}u(x,t)=0.\end{array}$$

(15)

Regarding its temporal dynamics, the function satisfies $u(x,·)\in C\left(\right(0,T\left]\right)$ and possesses continuous Hilfer fractional derivatives ${{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}u$ and ${{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}u$ on $(0,T]$. At the temporal boundary $t\to 0^{+}$, the solution naturally accommodates the singular behavior dictated by the fractional initial conditions in (13).

3.1. Solution in the Fourier-Time Domain

We now apply the n-dimensional Fourier transform with respect to the spatial variable $x\in {\mathbb{R}}^n$. Using (6), equation (12) becomes

$$\begin{array}{cc}\hfill & {{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}\widehat{u}(\xi ,t)+a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}\widehat{u}(\xi ,t)+\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right)\widehat{u}(\xi ,t)=\widehat{q}(\xi ,t)\hfill \end{array}$$

(16)

subject to the transformed conditions

$$\begin{array}{cc}\hfill & \left(I_{t,0^{+}}^{(1-{\mu }_2)(2-{\alpha }_2)}\widehat{u}\right)(\xi ,0^{+})={\widehat{g}}_1\left(\xi \right)\mathrm{and}{\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{(1-{\mu }_2)(2-{\alpha }_2)}\widehat{u}\right)(\xi ,0^{+})={\widehat{g}}_2\left(\xi \right).\hfill \end{array}$$

(17)

From (14), the Riemann-Lebesgue lemma and Plancherel’s theorem imply that

$$\begin{array}{c}\hfill \widehat{q}\in C([0,T];C_c\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right))\mathrm{and}{\widehat{g}}_1,{\widehat{g}}_2\in C_c\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right),\end{array}$$

(18)

from which we conclude that $\widehat{q}(\xi ,·)\in C\left([0,T]\right),$ for each fixed $\xi \in {\mathbb{R}}^n.$ We solve the Cauchy problem (16)-(17) by using the representation formulas in [19] for Hilfer-type fractional differential equations with variable coefficients. In particular, we combine the series representation of the solution of the associated homogeneous equation with the particular solution obtained in [19]. This leads to the following result.

Theorem 1.

Assume the hypothesis in (14) for the functions q, $g_1,g_2$, and $a_1,a_2,a_3$. For each fixed frequency $\xi \in {\mathbb{R}}^n$ and under the condition

$$\begin{array}{c}\hfill \left(\parallel a_1{\parallel }_{max}{I}_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+\parallel -a_2{\left|\xi \right|}^2+a_3{\parallel }_{max}{I}_{t,0^{+}}^{{\alpha }_2}\right)e^{\nu t}\le C{e}^{\nu t},\end{array}$$

(19)

for all $t\in [0,T]$ with ${\parallel ·\parallel }_{max}$ denoting the maximum norm in $C\left(\right[0,T\left]\right)$, and for some fixed $\nu \in {\mathbb{R}}^{+}$ and a constant $0<C<1$ independent of t, the Cauchy problem (16)-(17) admits a unique solution $\widehat{u}(\xi ,·)\in C^{1,{\alpha }_2}\left([0,T]\right)\cup \overline{C^{1,{\alpha }_2}\left((0,T]\right)}$, where

$$\begin{array}{c}\hfill C^{1,{\alpha }_2}\left([0,T]\right)=\left\{y\left(t\right)\in C\left([0,T]\right):{{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}y\left(t\right)\in C\left([0,T]\right)\right\},\end{array}$$

and $\overline{C^{1,{\alpha }_2}\left((0,T]\right)}$ denotes the closure of $C^{1,{\alpha }_2}\left((0,T]\right)$.

The solution $\widehat{u}(\xi ,t)$ can be decomposed as

$$\widehat{u}(\xi ,t)={\widehat{u}}_h(\xi ,t)+{\widehat{u}}_p(\xi ,t),$$

where ${\widehat{u}}_h$ is the solution of the associated homogeneous equation and ${\widehat{u}}_p$ is a particular solution of (16). Moreover,

$$\begin{array}{cc}\hfill & {\widehat{u}}_h(\xi ,t)={\widehat{g}}_1\left(\xi \right){\widehat{y}}_0(\xi ,t)+{\widehat{g}}_2\left(\xi \right){\widehat{y}}_1(\xi ,t),\hfill \end{array}$$

(20)

where, for each $j=0,1$,

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& ={\Psi }_j\left(t\right)-I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right),\hfill \end{array}$$

(21)

with

$$\begin{array}{c}\hfill {\Psi }_j\left(t\right)={\displaystyle \frac{t^{{\theta }_j}}{\Gamma \left(1+{\theta }_j\right)}}\mathrm{and}{\theta }_j=j-(1-{\mu }_2)(2-{\alpha }_2).\end{array}$$

For $j=0$, the representation (21) is valid provided that $(1-{\mu }_1)(1-{\alpha }_1)\ge (1-{\mu }_2)(2-{\alpha }_2)$. The particular solution is given by

$$\begin{array}{cc}\hfill & {\widehat{u}}_p(\xi ,t)=I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\delta \left(t\right){*}_t\widehat{q}(\xi ,t),\hfill \end{array}$$

(22)

where ${*}_t$ denotes the convolution in time (see (7)). The operators $\mathcal{A}$ and $\mathcal{B}$ are given by

$$\begin{array}{c}\hfill \mathcal{A}(\xi ,t):=a_1\left(t\right)I_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right)I_{t,0^{+}}^{{\alpha }_2}\end{array}$$

and

$$\begin{array}{c}\hfill \mathcal{B}(\xi ,t):=a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}+\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right).\end{array}$$

Remark 2.

The appearance of the convolution operator in (22) follows from the convolution theorem for the Laplace transform. Indeed, the inverse of the symbol of the operator in the Fourier-Laplace domain defines a convolution kernel in the time variable (see [26]). Moreover, the distribution $\delta \left(t\right)$ acts as the identity element with respect to Laplace convolution and therefore gives rise to the corresponding fundamental solution associated with the source term q.

The result follows directly from Theorems 3.2 and 3.4 of [19] applied to the transformed equation (16). For completeness, we verify that (20) and (22) define, respectively, the solution of the associated homogeneous equation and a particular solution of (16), both satisfying the initial conditions (17).

Proof. 

We first show that each function ${\widehat{y}}_j$ solves the homogeneous equation associated with (16). Using the identity ${{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{I}_{t,0^{+}}^{{\alpha }_2}f=f$ we obtain, for each $j=0,1$,

$$\begin{array}{cc}\hfill {{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\widehat{y}}_j(\xi ,t)& ={{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\Psi }_j\left(t\right)-{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right).\hfill \end{array}$$

(23)

Next, using the relation ${{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{I}_{t,0^{+}}^{{\alpha }_2}f=I_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}f$, we obtain

$$\begin{array}{c}\hfill {{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\widehat{y}}_j(\xi ,t)={{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\Psi }_j\left(t\right)-I_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right).\end{array}$$

Consequently,

$$\begin{array}{cc}\hfill & a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\widehat{y}}_j(\xi ,t)+\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right){\widehat{y}}_j(\xi ,t)\hfill \\ \hfill & =\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)+{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^{p+1}\mathcal{B}(\xi ,t){\Psi }_j\left(t\right).\hfill \end{array}$$

(24)

Adding (23) and (24), isolating the first term in the first series and shifting the index in the second series, we obtain

$$\begin{array}{cc}\hfill & {{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\widehat{y}}_j(\xi ,t)+a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\widehat{y}}_j(\xi ,t)+\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right){\widehat{y}}_j(\xi ,t)\hfill \\ \hfill & ={{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\Psi }_j\left(t\right)-\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)-{\displaystyle \sum _{p=1}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)\hfill \\ \hfill & +\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)+{\displaystyle \sum _{p=1}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)\hfill \\ \hfill & ={{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\Psi }_j\left(t\right).\hfill \end{array}$$

It remains to show that ${{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\Psi }_j\left(t\right)=0.$ Using the explicit form of ${\Psi }_j\left(t\right)$, the definition of the left Hilfer derivative given in (1), and the differentiation rule for power functions in (), we immediately obtain,

$$\begin{array}{c}\hfill {{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\Psi }_j\left(t\right)=I_{t,0^{+}}^{{\mu }_2(1-{\alpha }_2)}{\displaystyle \frac{d^2}{d{t}^2}}\left({\displaystyle \frac{t^j}{\Gamma \left(1+{\theta }_j\right)}}\right)=0,\end{array}$$

since the second derivative of $t^j$ vanishes for $j=0,1$. Hence, each ${\widehat{y}}_j$ solves the homogeneous equation associated with (16). Since ${\widehat{u}}_h$ is a linear combination of ${\widehat{y}}_0$ and ${\widehat{y}}_1$, it also satisfies the associated homogeneous equation.

We now verify that ${\widehat{u}}_p$ is a particular solution of (16). Using again the composition rule ${{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{I}_{t,0^{+}}^{{\alpha }_2}f=f$, we obtain

$$\begin{array}{cc}\hfill & {{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\widehat{u}}_p(\xi ,t)={\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\delta \left(t\right){*}_t\widehat{q}(\xi ,t).\hfill \end{array}$$

(25)

On the other hand, by the definition of $\mathcal{A}$, we have

$$\begin{array}{cc}\hfill & a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\widehat{u}}_p(\xi ,t)+\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right){\widehat{u}}_p(\xi ,t)\hfill \\ \hfill & =a_1\left(t\right)I_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\delta \left(t\right){*}_t\widehat{q}(\xi ,t)\hfill \\ \hfill & +\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right)I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\delta \left(t\right){*}_t\widehat{q}(\xi ,t)\hfill \\ \hfill & =-{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^{p+1}\delta \left(t\right){*}_t\widehat{q}(\xi ,t).\hfill \end{array}$$

(26)

Adding (25) and (26), separating the first term in the first series, and setting $r=p+1$ in the second series, yields

$$\begin{array}{cc}\hfill & {{}^H\partial }_{t,0^{+}}^{{\alpha }_2,{\mu }_2}{\widehat{u}}_p(\xi ,t)+a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\widehat{u}}_p(\xi ,t)+\left(-a_2\left(t\right){\left|\xi \right|}^2+a_3\left(t\right)\right){\widehat{u}}_p(\xi ,t)\hfill \\ \hfill & =\delta \left(t\right){*}_t\widehat{q}(\xi ,t)+{\displaystyle \sum _{p=1}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\delta \left(t\right){*}_t\widehat{q}(\xi ,t)\hfill \\ \hfill & -{\displaystyle \sum _{r=1}^{\infty }}{(-\mathcal{A}(\xi ,t))}^r\delta \left(t\right){*}_t\widehat{q}(\xi ,t)\hfill \\ \hfill & =\widehat{q}(\xi ,t).\hfill \end{array}$$

Thus, ${\widehat{u}}_p$ is a particular solution of (16).

To complete the proof, we need to verify that the formal solution given by the linear combination

$$\begin{array}{c}\hfill \widehat{u}(\xi ,t)={\widehat{g}}_1\left(\xi \right){\widehat{y}}_0(\xi ,t)+{\widehat{g}}_2\left(\xi \right){\widehat{y}}_1(\xi ,t)+{\widehat{u}}_p(\xi ,t)\end{array}$$

satisfies the initial conditions (17), for each $\xi \in {\mathbb{R}}^n$. For convenience, let $\gamma :=(1-{\mu }_2)(2-{\alpha }_2)\in [0,1[$. We analyze first the behavior of the homogeneous components ${\widehat{y}}_j(\xi ,t)$ as $t\to 0^{+}$. We decompose ${\widehat{y}}_j$ as

$$\begin{array}{c}\hfill {\widehat{y}}_j(\xi ,t)={\Psi }_j\left(t\right)-I_{t,0^{+}}^{{\alpha }_2}{\Phi }_j(\xi ,t),\end{array}$$

(27)

where ${\Psi }_j\left(t\right)={\displaystyle \frac{t^{{\theta }_j}}{\Gamma \left(1+{\theta }_j\right)}}$ with ${\theta }_j=j-\gamma$, and ${\Phi }_j(\xi ,t):={\sum }_{p=0}^{\infty }{(-\mathcal{A}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)$. Applying the fractional integral operator $I_{t,0^{+}}^{\gamma }$ to (27), and using the semigroup property, we obtain

$$\begin{array}{c}\hfill I_{t,0^{+}}^{\gamma }{\widehat{y}}_j(\xi ,t)=I_{t,0^{+}}^{\gamma }{\Psi }_j\left(t\right)-I_{t,0^{+}}^{\gamma +{\alpha }_2}{\Phi }_j(\xi ,t).\end{array}$$

(28)

Using (2), we have

$$\begin{array}{c}\hfill I_{t,0^{+}}^{\gamma }{\Psi }_j\left(t\right)=I_{t,0^{+}}^{\gamma }\left({\displaystyle \frac{t^{j-\gamma }}{\Gamma \left(1+j-\gamma \right)}}\right)={\displaystyle \frac{t^j}{\Gamma \left(1+j\right)}}.\end{array}$$

Thus, for $j=0$, we have $I_{t,0^{+}}^{\gamma }{\Psi }_0\left(t\right)=1$, and for $j=1$, we have $I_{t,0^{+}}^{\gamma }{\Psi }_1\left(t\right)=t$. Then, we obtain

$$\begin{array}{c}\hfill \underset{t\to 0^{+}}{lim}{I}_{t,0^{+}}^{\gamma }{\Psi }_j\left(t\right)={\delta }_{j0}\mathrm{and}\underset{t\to 0^{+}}{lim}{\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{\gamma }{\Psi }_j\left(t\right)\right)={\delta }_{j1},\end{array}$$

(29)

where ${\delta }_{ij}$ denotes the Kronecker delta. Second, we consider the remainder term in (28). Under the conditions of Theorem 1, the series defining ${\Phi }_j(\xi ,·)$ converges in $C\left(\right[0,T\left]\right)$. Hence, ${\Phi }_j(\xi ,·)$ is continuous on $[0,T]$. Since $\gamma +{\alpha }_2>0$, then by Lemma 1, we conclude that

$$\begin{array}{c}\hfill \underset{t\to 0^{+}}{lim}{I}_{t,0^{+}}^{\gamma +{\alpha }_2}{\Phi }_j(\xi ,t)=0.\end{array}$$

For the second initial condition, we obtain:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{\gamma }{\Phi }_j\right)(\xi ,t)=I_{t,0^{+}}^{\gamma +{\alpha }_2-1}{\Phi }_j(\xi ,t).\end{array}$$

Since $\gamma +{\alpha }_2-1>0$, by Lemma 1 we also conclude that

$$\begin{array}{c}\hfill \underset{t\to 0^{+}}{lim}{\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{\gamma }{\Phi }_j(\xi ,t)\right)=0.\end{array}$$

Therefore, evaluating (28) and its time derivative at $t\to 0^{+}$ yields

$$\begin{array}{cc}\hfill \underset{t\to 0^{+}}{lim}{I}_{t,0^{+}}^{\gamma }{\widehat{y}}_0(\xi ,t)={\delta }_{j0}\mathrm{and}& \underset{t\to 0^{+}}{lim}{\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{\gamma }{\widehat{y}}_0(\xi ,t)\right)={\delta }_{j1}.\hfill \end{array}$$

(30)

For the particular solution ${\widehat{u}}_p$, since $\widehat{q}(·,t)\in C\left([0,T]\right)$ and using similar arguments, we can easily conclude that

$$\begin{array}{c}\hfill \underset{t\to 0^{+}}{lim}{I}_{t,0^{+}}^{\gamma }{\widehat{u}}_p(\xi ,t)=0\mathrm{and}\underset{t\to 0^{+}}{lim}{\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{\gamma }{\widehat{u}}_p(\xi ,t)\right)=0.\end{array}$$

(31)

By substituting the limits (30) and (31) into the initial conditions (17), we readily verify:

$$\begin{array}{c}\hfill \underset{t\to 0^{+}}{lim}{I}_{t,0^{+}}^{\gamma }\widehat{u}(\xi ,t)={\widehat{g}}_1\left(\xi \right)\times 1+{\widehat{g}}_2\left(\xi \right)\times 0={\widehat{g}}_1\left(\xi \right),\\ \hfill \underset{t\to 0^{+}}{lim}{\displaystyle \frac{\partial }{\partial t}}\left(I_{t,0^{+}}^{\gamma }\widehat{u}(\xi ,t)\right)={\widehat{g}}_1\left(\xi \right)\times 0+{\widehat{g}}_2\left(\xi \right)\times 1={\widehat{g}}_2\left(\xi \right).\end{array}$$

This completes the proof. □

Remark 3.

To justify the convergence of the series representations in (21) and (22), we introduce the weighted Banach space

$$\begin{array}{c}\hfill C_{\nu }^{*}\left([0,T]\right)=\left\{y\in {C\left([0,T]\right):\parallel y\parallel }_{\nu }:=\underset{t\in [0,T]}{sup}{e}^{-\nu t}\left|y\left(t\right)\right|<\infty \right\},\end{array}$$

where $\nu >0$. In this setting, the operator $\mathcal{A}$ acts continuously on $C_{\nu }^{*}\left([0,T]\right)$. Moreover, we assume that $\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)\in C_{\nu }^{*}\left([0,T]\right)$. The series appearing in (21) and (22) are Neumann series associated with the resolvent operator ${(I+\mathcal{A})}^{-1}$. Consequently, the series converge absolutely in $C_{\nu }^{*}$ whenever $\mathcal{A}$ is a contraction, that is, ${\parallel \mathcal{A}\parallel }_{C_{\nu }^{*}}<1$. This contraction property was established in [19], and condition (19) provides a sufficient criterion for $\mathcal{A}$ to be a contraction on $C_{\nu }^{*}$.

3.2. Particular Cases

3.2.1. Constant Coefficients

We now consider the case of constant coefficients

$$\begin{array}{c}\hfill a_1\left(t\right)={\gamma }_1,a_2\left(t\right)=-{\gamma }_2,a_3\left(t\right)={\gamma }_3,\end{array}$$

(32)

with ${\gamma }_1,{\gamma }_2>0$, and ${\gamma }_3\ge 0$. We begin with the solution of the associated homogeneous equation (20):

$$\begin{array}{c}\hfill {\widehat{u}}_h(\xi ,t)={\widehat{g}}_1\left(\xi \right){\widehat{y}}_0(\xi ,t)+{\widehat{g}}_2\left(\xi \right){\widehat{y}}_1(\xi ,t),\end{array}$$

where, for each $j=0,1$,

$$\begin{array}{c}\hfill {\widehat{y}}_j(\xi ,t)={\Psi }_j\left(t\right)-I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right)\end{array}$$

(33)

and the operators $\mathcal{A}$ and $\mathcal{B}$ are given by

$$\begin{array}{c}\hfill \mathcal{A}(\xi ,t):={\gamma }_1{I}_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)I_{t,0^{+}}^{{\alpha }_2}\mathrm{and}\mathcal{B}(\xi ,t):={\gamma }_1{{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}+\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right).\end{array}$$

Since the coefficients are constant, the fractional integral operators appearing in $\mathcal{A}$ commute with $I_{t,0^{+}}^{{\alpha }_2}$. Therefore, applying the multinomial theorem to expand the powers of $\mathcal{A}$ and using the semigroup property of fractional integrals, we obtain

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& ={\Psi }_j\left(t\right)\hfill \\ \hfill & -{\displaystyle \sum _{p=0}^{\infty }}\sum _{l_1+l_2=p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{l_1,l_2}}\right){(-{\gamma }_1)}^{l_1}{\left(-\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)\right)}^{l_2}{I}_{t,0^{+}}^{{\alpha }_2+l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2}\mathcal{B}(\xi ,t){\Psi }_j\left(t\right),\hfill \end{array}$$

where $(l_1,l_2)\in {\mathbb{N}}_0^2$. Using the differentiation and integration formulas () and (2) applied to ${\Psi }_j$, we find

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& ={\displaystyle \frac{t^{{\theta }_j}}{\Gamma \left(1+{\theta }_j\right)}}-{\gamma }_1{t}^{{\theta }_j+{\alpha }_2-{\alpha }_1}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{{\mathbf{z}}_1^{l_1}{\mathbf{z}}_2^{l_2}}{\Gamma \left(1+{\theta }_j+{\alpha }_2-{\alpha }_1+l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2\right)}}\hfill \\ \hfill & -\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)t^{{\theta }_j+{\alpha }_2}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{{\mathbf{z}}_1^{l_1}{\mathbf{z}}_2^{l_2}}{\Gamma \left(1+{\theta }_j+{\alpha }_2+l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2\right)}},\hfill \end{array}$$

where ${\mathbf{z}}_1=-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1}$ and ${\mathbf{z}}_2=-\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)t^{{\alpha }_2}.$ Recognizing the series representation of the bivariate Mittag-Leffler function (see (8)), we obtain

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& =t^{{\theta }_j}\left({\displaystyle \frac{1}{\Gamma \left(1+{\theta }_j\right)}}+{\mathbf{z}}_1{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_j+{\alpha }_2-{\alpha }_1}\left({\mathbf{z}}_1,{\mathbf{z}}_2\right)+{\mathbf{z}}_2{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_j+{\alpha }_2}\left({\mathbf{z}}_1,{\mathbf{z}}_2\right)\right).\hfill \end{array}$$

Using the addition formula for the bivariate Mittag-Leffler function (see Lemma 2), we finally obtain

$$\begin{array}{c}\hfill {\widehat{y}}_j(\xi ,t)=t^{{\theta }_j}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_j}\left({\mathbf{z}}_1,{\mathbf{z}}_2\right).\end{array}$$

(34)

Hence, we get

$$\begin{array}{c}\hfill {\widehat{u}}_h(\xi ,t)={\widehat{g}}_1\left(\xi \right)t^{{\theta }_0}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_0}\left({\mathbf{z}}_1,{\mathbf{z}}_2\right)+{\widehat{g}}_2\left(\xi \right)t^{{\theta }_1}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_1}\left({\mathbf{z}}_1,{\mathbf{z}}_2\right),\end{array}$$

(35)

with ${\theta }_0=-(1-{\mu }_2)(2-{\alpha }_2)$ and ${\theta }_1=1-(1-{\mu }_2)(2-{\alpha }_2)$. This coincides with the Fourier-domain representation of the homogeneous solution obtained in [26] (see formula (41), second and third terms), with $c_2=1$, $c_1={\gamma }_1$, $c_0^2={\gamma }_2$, $d^2={\gamma }_3$, and $\psi \left(t\right)=t$.

The solution (34) can also be obtained if we simplify (33). Indeed, for $j=0$, the condition

$$\begin{array}{c}\hfill (1-{\mu }_1)(1-{\alpha }_1)\ge (1-{\mu }_2)(2-{\alpha }_2)\end{array}$$

required in Theorem 1 ensures that

$$\begin{array}{c}\hfill {{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\Psi }_j\left(t\right)={\displaystyle \frac{t^{{\theta }_j-{\alpha }_1}}{\Gamma \left(1+{\theta }_j-{\alpha }_1\right)}}=I_{t,0^{+}}^{-{\alpha }_1}{\Psi }_j\left(t\right)\end{array}$$

(36)

holds for each $j=0,1$. Therefore, since $I_{t,0^{+}}^{{\alpha }_2}$ commutes with all operators appearing in $\mathcal{A}$ and using (36), we can relate the operators $\mathcal{A}$ and $\mathcal{B}$ by

$$\begin{array}{c}\hfill I_{t,0^{+}}^{{\alpha }_2}\mathcal{B}(\xi ,t)={\gamma }_1{I}_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)I_{t,0^{+}}^{{\alpha }_2}=\mathcal{A}(\xi ,t).\end{array}$$

Hence, (33) simplifies to

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& ={\Psi }_j\left(t\right)-{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^{p+1}{\Psi }_j\left(t\right).\hfill \end{array}$$

Reindexing the series, we obtain

$$\begin{array}{c}\hfill {\widehat{y}}_j(\xi ,t)={\Psi }_j\left(t\right)-{\displaystyle \sum _{p=1}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p{\Psi }_j\left(t\right)={\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p{\Psi }_j\left(t\right).\end{array}$$

Hence, the homogeneous solution simplifies to

$$\begin{array}{c}\hfill {\widehat{u}}_h(\xi ,t)={\displaystyle \sum _{j=0}^1}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p{\Psi }_j\left(t\right){\widehat{g}}_{j+1}\left(\xi \right).\end{array}$$

(37)

Now, applying the multinomial theorem to expand the powers of $\mathcal{A}$ and using the semigroup property of fractional integrals and (2), we obtain

$$\begin{array}{cc}\hfill {\widehat{u}}_h(\xi ,t)& ={\displaystyle \sum _{j=0}^1}{\displaystyle \sum _{p=0}^{\infty }}\sum _{l_1+l_2=p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{l_1,l_2}}\right){(-{\gamma }_1)}^{l_1}{\left(-\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)\right)}^{l_2}{I}_{t,0^{+}}^{l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2}{\Psi }_j\left(t\right){\widehat{g}}_{j+1}\left(\xi \right)\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^1}{t}^{{\theta }_j}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \frac{{\mathbf{z}}_1^{l_1}{\mathbf{z}}_2^{l_2}}{\Gamma \left(1+{\theta }_j+l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2\right)}}{\widehat{g}}_{j+1}\left(\xi \right)\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^1}{t}^{{\theta }_j}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_j}\left({\mathbf{z}}_1,{\mathbf{z}}_2\right){\widehat{g}}_{j+1}\left(\xi \right),\hfill \end{array}$$

(38)

which coincides with (35). This shows that, in the constant-coefficient case, the homogeneous solution ${\widehat{u}}_h$ admits the simpler representation (37).

For the particular solution ${\widehat{u}}_p$, since the fractional integral operator $I_{t,0^{+}}^{{\alpha }_2}$ commutes with all operators appearing in $\mathcal{A}$, we have

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)=I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p\delta \left(t\right){*}_t\widehat{q}(\xi ,t)={\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p{I}_{t,0^{+}}^{{\alpha }_2}\delta \left(t\right){*}_t\widehat{q}(\xi ,t).\end{array}$$

Since

$$\begin{array}{c}\hfill I_{t,0^{+}}^{{\alpha }_2}\delta \left(t\right)={\displaystyle \frac{t^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}}\end{array}$$

(39)

it follows that

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)={\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{A}(\xi ,t))}^p{\displaystyle \frac{t^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}}{*}_t\widehat{q}(\xi ,t).\end{array}$$

Now, applying the multinomial expansion together with the semigroup property of fractional integrals and (2), we obtain

$$\begin{array}{cc}\hfill & {\widehat{u}}_p(\xi ,t)\hfill \\ \hfill & ={\displaystyle \sum _{p=0}^{\infty }}\sum _{\begin{array}{c}{l}_1+l_2=p\\ l_1,l_2\in {\mathbb{N}}_0\end{array}}{(-1)}^{l_1+l_2}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{l_1,l_2}}\right){\gamma }_1^{l_1}{\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)}^{l_2}{I}_{t,0^{+}}^{l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2}{\displaystyle \frac{t^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}}{*}_t\widehat{q}(\xi ,t)\hfill \\ \hfill & =t^{{\alpha }_2-1}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{(-{\gamma }_1)}^{l_1}{\left(-\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)\right)}^{l_2}{\displaystyle \frac{t^{l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2}}{\Gamma \left({\alpha }_2+l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2\right)}}{*}_t\widehat{q}(\xi ,t).\hfill \end{array}$$

Using the series representation (8), we rewrite

$$\begin{array}{cc}\hfill {\widehat{u}}_p(\xi ,t)& =\left(t^{{\alpha }_2-1}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),{\alpha }_2}\left(-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1},-\left({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3\right)t^{{\alpha }_2}\right)\right){*}_t\widehat{q}(\xi ,t).\hfill \end{array}$$

(40)

This coincides with the Fourier-domain representation of the particular solution derived in [26] (see formula (41), last term), with $c_2=1$, $c_1={\gamma }_1$, $c_0^2={\gamma }_2$, $d^2={\gamma }_3$, and $\psi \left(t\right)=t$.

3.2.2. Power-Law Coefficients

We now consider the case of power-law coefficients

$$\begin{array}{c}\hfill a_1\left(t\right)={\gamma }_1{t}^{{\beta }_1},a_2\left(t\right)=-{\gamma }_2{t}^{{\beta }_2},a_3\left(t\right)={\gamma }_3{t}^{{\beta }_3},\end{array}$$

(41)

where ${\beta }_i\ge 0$ for $i=1,2,3$, ${\gamma }_1,{\gamma }_2>0$, and ${\gamma }_3\ge 0$. We next verify that condition (19) is satisfied. Recall that

$$\begin{array}{c}\hfill I_{t,0^{+}}^{\alpha }{e}^{\nu t}={\displaystyle \frac{\gamma (\alpha ,\nu t)}{\Gamma \left(\alpha \right){\nu }^{\alpha }}}{e}^{\nu t},\end{array}$$

where $\gamma (\alpha ,·)$ denotes the lower incomplete Gamma function. Since $\gamma (\alpha ,x)\le \Gamma \left(\alpha \right),$ for all $x\ge 0,$ it follows that, for all $t\ge 0$,

$$\begin{array}{c}\hfill I_{t,0^{+}}^{\alpha }{e}^{\nu t}\le {\displaystyle \frac{e^{\nu t}}{{\nu }^{\alpha }}}.\end{array}$$

Therefore, condition (19) reduces to

$$\begin{array}{c}\hfill \left({\displaystyle \frac{\parallel a_1{\parallel }_{max}}{{\nu }^{{\alpha }_2-{\alpha }_1}}}+{\displaystyle \frac{\parallel -a_2{\left|\xi \right|}^2+a_3{\parallel }_{max}}{{\nu }^{{\alpha }_2}}}\right)e^{\nu t}\le C{e}^{\nu t}.\end{array}$$

(42)

Hence, for each fixed $\xi \in {\mathbb{R}}^n$, there exists $\nu >0$ such that

$$\begin{array}{c}\hfill {\displaystyle \frac{\parallel a_1{\parallel }_{max}}{{\nu }^{{\alpha }_2-{\alpha }_1}}}+{\displaystyle \frac{\parallel -a_2{\left|\xi \right|}^2+a_3{\parallel }_{max}}{{\nu }^{{\alpha }_2}}}<1\end{array}$$

and consequently condition (42) holds for all $t\in [0,T]$, with a constant $0<C<1$ independent of t.

We consider first the particular solution (22):

$$\begin{array}{cc}\hfill & {\widehat{u}}_p(\xi ,t)=I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p{\left({\gamma }_1{t}^{{\beta }_1}{I}_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+{\gamma }_2{\left|\xi \right|}^2{t}^{{\beta }_2}{I}_{t,0^{+}}^{{\alpha }_2}+{\gamma }_3{t}^{{\beta }_3}{I}_{t,0^{+}}^{{\alpha }_2}\right)}^p\delta \left(t\right){*}_t\widehat{q}(\xi ,t).\hfill \end{array}$$

(43)

For convenience, we introduce the operators

$$\begin{array}{c}\hfill Y_i={\gamma }_i{t}^{{\beta }_i}{I}_{t,0^{+}}^{{\delta }_i},i=1,2,3,\end{array}$$

where ${\delta }_1={\alpha }_2-{\alpha }_1,{\delta }_2={\delta }_3={\alpha }_2,$ and define the operator

$$\begin{array}{cc}\hfill & \mathcal{L}(\xi ,t):=Y_1+{\left|\xi \right|}^2{Y}_2+Y_3.\hfill \end{array}$$

Since the operators $Y_i$ do not commute, the expansion of ${\mathcal{L}}^p(\xi ,t)$ requires the use of the non-commutative multinomial formula. More precisely, we have

$$\begin{array}{c}\hfill {\left(Y_1+\dots +Y_r\right)}^p=\sum _{\left|l\right|=p}\sum _{\pi \in \mathfrak{S}\left(l\right)}{Y}_{i_{\pi \left(1\right)}}\dots Y_{i_{\pi \left(p\right)}},\end{array}$$

where $l=(l_1,l_2,\dots ,l_r)\in {\mathbb{N}}_0^r$ is a multi-index and $\mathfrak{S}\left(l\right)$ denotes the set of distinct permutations of the multiset $\{i_1,\dots ,i_p\}$. Applying this expansion, we obtain

$$\begin{array}{c}\hfill {\mathcal{L}}^p(\xi ,t)=\sum _{\left|l\right|=p}{\left|\xi \right|}^{2l_2}\sum _{\pi \in \mathfrak{S}\left(l\right)}{Y}_{i_{\pi \left(1\right)}}\dots Y_{i_{\pi \left(p\right)}},\end{array}$$

(44)

where $l=(l_1,l_2,l_3)\in {\mathbb{N}}_0^3$. Substituting (44) into (43), we find

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)={\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\left|l\right|=p}{\left|\xi \right|}^{2l_2}{I}_{t,0^{+}}^{{\alpha }_2}\sum _{\pi \in \mathfrak{S}\left(l\right)}{Y}_{i_{\pi \left(1\right)}}\dots Y_{i_{\pi \left(p\right)}}\delta \left(t\right){*}_t\widehat{q}(\xi ,t).\end{array}$$

By using (39), we obtain

$$\begin{array}{c}\hfill Y_{i_{\pi \left(p\right)}}\delta \left(t\right)={\gamma }_{i_{\pi \left(p\right)}}{\displaystyle \frac{t^{{\beta }_{i_{\pi \left(p\right)}}+{\delta }_{i_{\pi \left(p\right)}}-1}}{\Gamma \left({\delta }_{i_{\pi \left(p\right)}}\right)}}.\end{array}$$

Since each operator $Y_{i_{\pi \left(p\right)}}$ consists of a multiplication operator followed by a fractional integral operator, a recursive application of (2) yields the explicit formula

$$\begin{array}{c}\hfill I_{t,0^{+}}^{{\alpha }_2}{Y}_{i_{\pi \left(1\right)}}\dots Y_{i_{\pi \left(p\right)}}\delta \left(t\right)={\gamma }_{\pi }{\displaystyle \prod _{m=0}^{p-1}}{\displaystyle \frac{\Gamma \left({\lambda }_m\right)}{\Gamma \left({\lambda }_m+{\delta }_{i_{\pi \left(m\right)}}\right)}}{\displaystyle \frac{t^{{\lambda }_0+{\alpha }_2-1}}{\Gamma \left({\delta }_{i_{\pi \left(p\right)}}\right)}},\end{array}$$

(45)

where ${\gamma }_{\pi }={\gamma }_{i_{\pi \left(1\right)}}\dots {\gamma }_{i_{\pi \left(p\right)}}={\gamma }_1^{l_1}{\gamma }_2^{l_2}{\gamma }_3^{l_3}:={\gamma }_l$ and ${\lambda }_m={\sum }_{r=m+1}^p({\beta }_{i_{\pi \left(r\right)}}+{\delta }_{i_{\pi \left(r\right)}})$, for $m=0,\dots ,p-1$, with ${\delta }_{i_{\pi \left(0\right)}}={\alpha }_2$. Hence, we obtain

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)=t^{{\alpha }_2-1}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\left|l\right|=p}{\left|\xi \right|}^{2l_2}{\gamma }_l\sum _{\pi \in \mathfrak{S}\left(l\right)}{\displaystyle \prod _{m=0}^{p-1}}{\displaystyle \frac{\Gamma \left({\lambda }_m\right)}{\Gamma \left({\lambda }_m+{\delta }_{i_{\pi \left(m\right)}}\right)}}{\displaystyle \frac{t^{{\lambda }_0}}{\Gamma \left({\delta }_{i_{\pi \left(p\right)}}\right)}}{*}_t\widehat{q}(\xi ,t),\end{array}$$

(46)

We now verify that (46) reduces to (40) when ${\beta }_1={\beta }_2={\beta }_3=0$. In this case, ${\lambda }_m={\sum }_{r=m+1}^p{\delta }_{i_{\pi \left(r\right)}}$, which implies the recursive relation ${\lambda }_m={\lambda }_{m+1}+{\delta }_{i_{\pi (m+1)}}$. Consequently, in the product appearing in (45), the numerator of each factor cancels with the denominator of the subsequent one. Therefore, (46) simplifies to

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)=t^{{\alpha }_2-1}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\left|l\right|=p}{\left|\xi \right|}^{2l_2}{\gamma }_l\left|\mathfrak{S}\left(l\right)\right|{\displaystyle \frac{t^{{\lambda }_0}}{\Gamma \left({\lambda }_0+{\alpha }_2\right)}}{*}_t\widehat{q}(\xi ,t).\end{array}$$

Since

$$\begin{array}{c}\hfill \left|\mathfrak{S}\left(l\right)\right|={\displaystyle \frac{(l_1+l_2+l_3)!}{l_1!l_2!l_3!}}\mathrm{and}{\lambda }_0=l_1{\delta }_1+l_2{\delta }_2+l_3{\delta }_3=l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2+l_3{\alpha }_2\end{array}$$

we obtain

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)=t^{{\alpha }_2-1}{\displaystyle \sum _{l_1,l_2,l_3=0}^{\infty }}{\displaystyle \frac{(l_1+l_2+l_3)!}{l_1!l_2!l_3!}}{\displaystyle \frac{{(-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1})}^{l_1}(-{\gamma }_2{\left|\xi \right|}^2{t}^{{\alpha }_2}{)}^{l_2}{(-{\gamma }_3{t}^{{\alpha }_2})}^{l_3}}{\Gamma \left({\alpha }_2+l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2+l_3{\alpha }_2\right)}}{*}_t\widehat{q}(\xi ,t).\end{array}$$

(47)

Recognizing the series definition of the trivariate ML function (see ()), we obtain

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)=\left(t^{{\alpha }_2-1}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2,{\alpha }_2),{\alpha }_2}\left(z_1,z_2,z_3\right)\right){*}_t\widehat{q}(\xi ,t).\end{array}$$

with $z_1=-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1}$, $z_2=-{\gamma }_2{\left|\xi \right|}^2{t}^{{\alpha }_2}$, and $z_3=-{\gamma }_3{t}^{{\alpha }_2}$. Now, applying the reduction formula (11), we get

$$\begin{array}{c}\hfill {\widehat{u}}_p(\xi ,t)=\left(t^{{\alpha }_2-1}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),{\alpha }_2}\left(z_1,z_2+z_3\right)\right){*}_t\widehat{q}(\xi ,t),\end{array}$$

which coincides exactly with (40).

We now turn to the solution of the associated homogeneous equation

$$\begin{array}{c}\hfill {\widehat{u}}_h(\xi ,t)={\widehat{g}}_1\left(\xi \right){\widehat{y}}_0(\xi ,t)+{\widehat{g}}_2\left(\xi \right){\widehat{y}}_1(\xi ,t).\end{array}$$

Substituting (41) into (21), we have, for each $j=0,1$,

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& ={\Psi }_j\left(t\right)-I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-\mathcal{L}(\xi ,t))}^p\mathcal{B}(\xi ,t){\Psi }_j\left(t\right).\hfill \end{array}$$

The application of $\mathcal{B}$ to ${\Psi }_j$ yields

$$\begin{array}{c}\hfill \mathcal{B}(\xi ,t){\Psi }_j\left(t\right)={\displaystyle \frac{{\gamma }_1{t}^{{\theta }_j-{\alpha }_1+{\beta }_1}}{\Gamma \left(1+{\theta }_j-{\alpha }_1\right)}}+{\left|\xi \right|}^2{\displaystyle \frac{{\gamma }_2{t}^{{\theta }_j+{\beta }_2}}{\Gamma \left(1+{\theta }_j\right)}}+{\displaystyle \frac{{\gamma }_3{t}^{{\theta }_j+{\beta }_3}}{\Gamma \left(1+{\theta }_j\right)}}.\end{array}$$

Now, using arguments similar to those leading to (45), a recursive application of (2) yields

$$\begin{array}{cc}\hfill & {\widehat{y}}_j(\xi ,t)={\Psi }_j\left(t\right)\hfill \\ \hfill & -{\gamma }_1{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\begin{array}{c}\left|l\right|=p\end{array}}{\left|\xi \right|}^{2l_2}{\gamma }_l\sum _{\pi \in \mathfrak{S}\left(l\right)}{\displaystyle \prod _{m=0}^p}{\displaystyle \frac{\Gamma \left({\lambda }_m+{\theta }_j-{\alpha }_1+{\beta }_1+1\right)}{\Gamma \left({\lambda }_m+{\delta }_{i_{\pi \left(m\right)}}+{\theta }_j-{\alpha }_1+{\beta }_1+1\right)}}{\displaystyle \frac{t^{{\alpha }_2+{\lambda }_0+{\theta }_j-{\alpha }_1+{\beta }_1}}{\Gamma \left(1+{\theta }_j-{\alpha }_1\right)}}\hfill \\ \hfill & -{\gamma }_2{\left|\xi \right|}^2{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\begin{array}{c}\left|l\right|=p\end{array}}{\left|\xi \right|}^{2l_2}{\gamma }_l\sum _{\pi \in \mathfrak{S}\left(l\right)}{\displaystyle \prod _{m=0}^p}{\displaystyle \frac{\Gamma \left({\lambda }_m+{\theta }_j+{\beta }_2+1\right)}{\Gamma \left({\lambda }_m+{\delta }_{i_{\pi \left(m\right)}}+{\theta }_j+{\beta }_2+1\right)}}{\displaystyle \frac{t^{{\alpha }_2+{\lambda }_0+{\theta }_j+{\beta }_2}}{\Gamma \left(1+{\theta }_j\right)}}\hfill \\ \hfill & -{\gamma }_3{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\begin{array}{c}\left|l\right|=p\end{array}}{\left|\xi \right|}^{2l_2}{\gamma }_l\sum _{\pi \in \mathfrak{S}\left(l\right)}{\displaystyle \prod _{m=0}^p}{\displaystyle \frac{\Gamma \left({\lambda }_m+{\theta }_j+{\beta }_3+1\right)}{\Gamma \left({\lambda }_m+{\delta }_{i_{\pi \left(m\right)}}+{\theta }_j+{\beta }_3+1\right)}}{\displaystyle \frac{t^{{\alpha }_2+{\lambda }_0+{\theta }_j+{\beta }_3}}{\Gamma \left(1+{\theta }_j\right)}},\hfill \end{array}$$

(48)

where ${\gamma }_l={\gamma }_1^{l_1}{\gamma }_2^{l_2}{\gamma }_3^{l_3}$ and ${\lambda }_m={\sum }_{r=m+1}^p({\beta }_{i_{\pi \left(r\right)}}+{\delta }_{i_{\pi \left(r\right)}})$, for $m=0,\dots ,p$, with the convention ${\lambda }_p=0$ and ${\delta }_{i_{\pi \left(0\right)}}={\alpha }_2$.

We now verify that (48) reduces to (34) when ${\beta }_1={\beta }_2={\beta }_3=0$. In this case, ${\lambda }_m={\sum }_{r=m+1}^p{\delta }_{i_{\pi \left(r\right)}}$, which implies the recursive relation ${\lambda }_m={\lambda }_{m+1}+{\delta }_{i_{\pi (m+1)}}$. Consequently, in each product appearing in (48), the numerator of each factor cancels with the denominator of the subsequent one. Therefore, the sums over $\mathfrak{S}\left(l\right)$ collapse to

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& ={\Psi }_j\left(t\right)-{\gamma }_1{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\begin{array}{c}\left|l\right|=p\end{array}}{\left|\xi \right|}^{2l_2}{\gamma }_l\left|\mathfrak{S}\left(l\right)\right|{\displaystyle \frac{t^{{\alpha }_2+{\lambda }_0+{\theta }_j-{\alpha }_1}}{\Gamma \left({\lambda }_0+{\alpha }_2+{\theta }_j-{\alpha }_1+1\right)}}\hfill \\ \hfill & -({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3){\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\sum _{\begin{array}{c}\left|l\right|=p\end{array}}{\left|\xi \right|}^{2l_2}{\gamma }_l\left|\mathfrak{S}\left(l\right)\right|{\displaystyle \frac{t^{{\alpha }_2+{\lambda }_0+{\theta }_j}}{\Gamma \left({\lambda }_0+{\alpha }_2+{\theta }_j+1\right)}}.\hfill \end{array}$$

Since

$$\begin{array}{c}\hfill \left|\mathfrak{S}\left(l\right)\right|={\displaystyle \frac{(l_1+l_2+l_3)!}{l_1!l_2!l_3!}}\mathrm{and}{\lambda }_0=l_1{\delta }_1+l_2{\delta }_2+l_3{\delta }_3=l_1({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2+l_3{\alpha }_2\end{array}$$

and recognizing the series definition of the trivariate ML function (see ()), we obtain

$$\begin{array}{cc}\hfill {\widehat{y}}_j(\xi ,t)& ={\Psi }_j\left(t\right)-{\gamma }_1{t}^{{\theta }_j+{\alpha }_2-{\alpha }_1}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2,{\alpha }_2),1+{\theta }_j+{\alpha }_2-{\alpha }_1}\left(z_1,z_2,z_3\right)\hfill \\ \hfill & -({\gamma }_2{\left|\xi \right|}^2+{\gamma }_3)t^{{\theta }_j+{\alpha }_2}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2,{\alpha }_2),1+{\theta }_j+{\alpha }_2}\left(z_1,z_2,z_3\right)\hfill \end{array}$$

with $z_1=-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1}$, $z_2=-{\gamma }_2{\left|\xi \right|}^2{t}^{{\alpha }_2}$, and $z_3=-{\gamma }_3{t}^{{\alpha }_2}$. Substituting ${\Psi }_j\left(t\right)$ and applying the reduction formula (11), we obtain

$$\begin{array}{c}\hfill {\widehat{y}}_j(\xi ,t)=t^{{\theta }_j}\left({\displaystyle \frac{1}{\Gamma \left(1+{\theta }_j\right)}}+z_1{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_j+{\alpha }_2-{\alpha }_1}\left(z_1,{\mathbf{z}}_2\right)+{\mathbf{z}}_2{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_j+{\alpha }_2}\left(z_1,{\mathbf{z}}_2\right)\right)\end{array}$$

where ${\mathbf{z}}_2=z_2+z_3$. Finally, applying the addition formula for the bivariate ML function (see Lemma 2), the terms inside the brackets combine into a single bivariate Mittag-Leffler function, yielding

$$\begin{array}{c}\hfill {\widehat{y}}_j(\xi ,t)=t^{{\theta }_j}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),1+{\theta }_j}\left(z_1,{\mathbf{z}}_2\right).\end{array}$$

This result coincides with the expression (34), which confirms the consistency of our results.

3.3. Solution in the Space-Time Domain

We now recover the solution to (12)-(13) in the space-time domain by applying the inverse Fourier transform to its representation in the Fourier domain. Since the Fourier transform acts only on the spatial variable, the time variable is treated as a parameter during the inversion. Hence,

$$\begin{array}{c}\hfill u(x,t)={\mathcal{F}}^{-1}\left\{{\widehat{u}}_h(\xi ,t)\right\}\left(x,t\right)+{\mathcal{F}}^{-1}\left\{{\widehat{u}}_p(\xi ,t)\right\}\left(x,t\right),\end{array}$$

(49)

where ${\widehat{u}}_h$ and ${\widehat{u}}_p$ are given by (20) and (22), respectively.

We first consider the associated homogeneous solution. By linearity of the inverse Fourier transform and the convolution theorem (see (4)), we obtain

$$\begin{array}{cc}\hfill u_h(x,t)& ={\displaystyle \sum _{j=0}^1}{\Psi }_j\left(t\right)\delta \left(x\right)-I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{j=0}^1}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p\left[{\mathcal{A}}^{*p}(x,t){*}_x\mathcal{B}(x,t){\Psi }_j\left(t\right)\right]{*}_x{g}_{j+1}\left(x\right)\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^1}\left({\Psi }_j\left(t\right)\delta \left(x\right)-I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p{\mathcal{A}}^{*p}(x,t){*}_x\mathcal{B}(x,t){\Psi }_j\left(t\right)\right){*}_x{g}_{j+1}\left(x\right),\hfill \end{array}$$

(50)

where ${*}_x$ denotes the spatial convolution defined in (5). The operators $\mathcal{A}$ and $\mathcal{B}$ in the space-time domain are given by

$$\begin{array}{c}\hfill \mathcal{A}(x,t)={\mathcal{F}}^{-1}\left\{\mathcal{A}(\xi ,t)\right\}\left(x,t\right)=\left(a_1\left(t\right)I_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+(a_2\left(t\right){\Delta }_x+a_3\left(t\right))I_{t,0^{+}}^{{\alpha }_2}\right)\delta \left(x\right)\end{array}$$

(51)

and

$$\begin{array}{c}\hfill \mathcal{B}(x,t)={\mathcal{F}}^{-1}\left\{\mathcal{B}(\xi ,t)\right\}\left(x,t\right)=\left(a_1\left(t\right){{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}+a_2\left(t\right){\Delta }_x+a_3\left(t\right)\right)\delta \left(x\right).\end{array}$$

(52)

Since we are working in the distributional setting, the Dirac distribution $\delta \left(x\right)$ is attached to the spatial differential operator ${\Delta }_x$ so that the convolution products are well defined.

For $p\in \mathbb{N}$, the symbol ${\mathcal{A}}^{*p}$ denotes the p-fold spatial convolution of $\mathcal{A}$ with itself, namely

$$\begin{array}{c}\hfill {\mathcal{A}}^{*p}(x,t):=\left(\underset{p\mathrm{times}}{\underbrace{\mathcal{A}(·,t){*}_x\dots {*}_x\mathcal{A}(·,t)}}\right)(x,t)\end{array}$$

with the convention ${\mathcal{A}}^{*0}(x,t)=\delta \left(x\right)$.

The particular solution is obtained analogously. Applying the inverse Fourier transform to (22) and using the convolution theorem, we obtain

$$\begin{array}{c}\hfill u_p(x,t)=I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p{\mathcal{A}}^{*p}(x,t)\delta \left(t\right)\delta \left(x\right){*}_{x,t}q(x,t),\end{array}$$

(53)

where ${*}_{x,t}$ denotes the space-time convolution, that is, the temporal convolution ${*}_t$ followed by the spatial convolution ${*}_x$. We summarize the above results in the following theorem.

Theorem 2

(Distributional space-time representation). Let the assumptions of Theorem 1 be satisfied. Then the solution of problem (12)-(13) admits the distributional representation

$$\begin{array}{c}\hfill u(x,t)=u_h(x,t)+u_p(x,t),\end{array}$$

where

$$\begin{array}{c}\hfill u_h(x,t)={\displaystyle \sum _{j=0}^1}\left({\Psi }_j\left(t\right)\delta \left(x\right)-I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p{\mathcal{A}}^{*p}(x,t){*}_x\mathcal{B}(x,t){\Psi }_j\left(t\right)\right){*}_x{g}_{j+1}\left(x\right),\end{array}$$

(54)

and

$$\begin{array}{c}\hfill u_p(x,t)=I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{p=0}^{\infty }}{(-1)}^p{\mathcal{A}}^{*p}(x,t)\delta \left(t\right)\delta \left(x\right){*}_{x,t}q(x,t).\end{array}$$

(55)

Here, the operators $\mathcal{A}(x,t)$ and $\mathcal{B}(x,t)$ are given by (51) and (52), respectively.

3.3.1. Constant Coefficients

For constant coefficients, the general distributional representations (54) and (55) simplify.

Operator representation

We first derive a compact operator representation of the homogeneous solution. Let $a_1\left(t\right)={\gamma }_1$, $a_2\left(t\right)=-{\gamma }_2$, and $a_3\left(t\right)={\gamma }_3$, with ${\gamma }_1,{\gamma }_2>0$ and ${\gamma }_3\ge 0$. We decompose the operator $\mathcal{A}$ as

$$\begin{array}{c}\hfill \mathcal{A}(x,t)={\mathcal{A}}_1\left(t\right)\delta \left(x\right)+{\mathcal{A}}_2\left(t\right)(-{\Delta }_x)\delta \left(x\right),\end{array}$$

where

$$\begin{array}{c}\hfill {\mathcal{A}}_1\left(t\right)={\gamma }_1{I}_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+{\gamma }_3{I}_{t,0^{+}}^{{\alpha }_2}\mathrm{and}{\mathcal{A}}_2\left(t\right)={\gamma }_2{I}_{t,0^{+}}^{{\alpha }_2}.\end{array}$$

Since both ${\mathcal{A}}_1$ and ${\mathcal{A}}_2$ are linear combinations of Riemann-Liouville fractional integral operators, they commute by the semigroup property. Consequently, the convolution powers ${\mathcal{A}}^{*p}$ can be identified with the corresponding operator powers ${\mathcal{A}}^p$ acting on $\delta \left(x\right)$. The operator $\mathcal{B}$ can be written as

$$\begin{array}{c}\hfill \mathcal{B}(x,t)={\mathcal{B}}_1\left(t\right)\delta \left(x\right)+{\gamma }_2(-{\Delta }_x)\delta \left(x\right),\end{array}$$

with

$$\begin{array}{c}\hfill {\mathcal{B}}_1\left(t\right)={\gamma }_1{{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}+{\gamma }_3.\end{array}$$

Using the binomial expansion and recalling that $\delta \left(x\right)$ acts as the identity for convolution, we obtain

$$\begin{array}{c}\hfill {\mathcal{A}}^p(x,t){*}_x\left(\mathcal{B}(x,t){\Psi }_j\left(t\right)\right)=\sum _{l_1+l_2=p}{\displaystyle \frac{p!}{l_1!l_2!}}{\mathcal{A}}_1{\left(t\right)}^{l_1}{\mathcal{A}}_2{\left(t\right)}^{l_2}{(-{\Delta }_x)}^{l_2}\delta \left(x\right){*}_x\mathcal{B}(x,t){\Psi }_j\left(t\right).\end{array}$$

(56)

Taking into account the expression of $\mathcal{B}(x,t)$ and the fact that ${\Psi }_j$ depends only on the variable t, we obtain

$$\begin{array}{c}\hfill {(-{\Delta }_x)}^{l_2}\delta \left(x\right){*}_x\mathcal{B}(x,t){\Psi }_j\left(t\right)=\left({\mathcal{B}}_1\left(t\right){(-{\Delta }_x)}^{l_2}\delta \left(x\right)+{\gamma }_2{(-{\Delta }_x)}^{l_2+1}\delta \left(x\right)\right){\Psi }_j\left(t\right).\end{array}$$

(57)

Substituting (56) and (57) into (50), we obtain

$$\begin{array}{cc}\hfill & u_h(x,t)\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^1}({\Psi }_j\left(t\right)\delta \left(x\right)-I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{(-1)}^{l_1+l_2}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\mathcal{A}}_1{\left(t\right)}^{l_1}{\mathcal{A}}_2{\left(t\right)}^{l_2}{\mathcal{B}}_1\left(t\right){(-{\Delta }_x)}^{l_2}\delta \left(x\right){\Psi }_j\left(t\right)\hfill \\ \hfill & -{\gamma }_2{I}_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{(-1)}^{l_1+l_2}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\mathcal{A}}_1{\left(t\right)}^{l_1}{\mathcal{A}}_2{\left(t\right)}^{l_2}{(-{\Delta }_x)}^{l_2+1}\delta \left(x\right){\Psi }_j\left(t\right)){*}_x{g}_{j+1}\left(x\right).\hfill \end{array}$$

Recalling that

$$\begin{array}{c}\hfill {{}^H\partial }_{t,0^{+}}^{{\alpha }_1,{\mu }_1}{\Psi }_j\left(t\right)=I_{t,0^{+}}^{-{\alpha }_1}{\Psi }_j\left(t\right),\end{array}$$

we have

$$\begin{array}{c}\hfill I_{t,0^{+}}^{{\alpha }_2}{\mathcal{B}}_1\left(t\right)={\gamma }_1{I}_{t,0^{+}}^{{\alpha }_2-{\alpha }_1}+{\gamma }_3{I}_{t,0^{+}}^{{\alpha }_2}={\mathcal{A}}_1\left(t\right).\end{array}$$

Using this identity, we obtain

$$\begin{array}{cc}\hfill u_h(x,t)& ={\displaystyle \sum _{j=0}^1}({\Psi }_j\left(t\right)\delta \left(x\right)-{\mathcal{A}}_1\left(t\right){\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{(-{\mathcal{A}}_1\left(t\right))}^{l_1}{\left({\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{l_2}\delta \left(x\right){\Psi }_j\left(t\right)\hfill \\ \hfill & +{\mathcal{A}}_2\left(t\right){\Delta }_x{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{(-{\mathcal{A}}_1\left(t\right))}^{l_1}{\left({\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{l_2}\delta \left(x\right){\Psi }_j\left(t\right)){*}_x{g}_{j+1}\left(x\right).\hfill \end{array}$$

Using the binomial identity (10) and arguments similar to those employed in Lemma 2, the sum inside the parentheses reduces to

$$\begin{array}{c}\hfill u_h(x,t)={\displaystyle \sum _{j=0}^1}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{(-{\mathcal{A}}_1\left(t\right))}^{l_1}{\left({\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{l_2}\delta \left(x\right){\Psi }_j\left(t\right){*}_x{g}_{j+1}\left(x\right).\end{array}$$

(58)

This expression can also be obtained by inverting the Fourier transform in (37) and expanding the operator powers ${\mathcal{A}}^p$. This shows that all the representations obtained are consistent with one another. Observing that the double series appearing in (58) corresponds to a binomial-type expansion, we recognize that it can be written formally as

$$\begin{array}{c}\hfill {\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{(-{\mathcal{A}}_1\left(t\right))}^{l_1}{\left({\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{l_2}={\left(I+{\mathcal{A}}_1\left(t\right)-{\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{-1}.\end{array}$$

Therefore,

$$\begin{array}{c}\hfill u_h(x,t)={\displaystyle \sum _{j=0}^1}{\left(I+{\mathcal{A}}_1\left(t\right)-{\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{-1}{\Psi }_j\left(t\right){*}_x{g}_{j+1}\left(x\right).\end{array}$$

The previous representation gives a compact operator form of the homogeneous solution, where the inverse operator is understood in terms of its convergent Neumann-type series expansion, whose convergence has already been established. An alternative closed-form operational representation for $u_h$ is given by the Fourier inversion of (38):

$$\begin{array}{c}\hfill u_h(x,t)={\displaystyle \sum _{j=0}^1}{t}^{{\theta }_j}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),{\theta }_j+1}\left(-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1},-\left(-{\gamma }_2{\Delta }_x+{\gamma }_3\right)t^{{\alpha }_2}\right)\delta \left(x\right){*}_x{g}_{j+1}\left(x\right).\end{array}$$

(59)

For the particular solution (53) we get the operator representation:

$$\begin{array}{c}\hfill u_p(x,t)=I_{t,0^{+}}^{{\alpha }_2}{\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{(-{\mathcal{A}}_1\left(t\right))}^{l_1}{\left({\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{l_2}\delta \left(t\right)\delta \left(x\right){*}_{x,t}q(x,t).\end{array}$$

(60)

Alternatively, from the Fourier inversion of (40), we get the closed-form operational representation:

$$\begin{array}{c}\hfill u_p(x,t)=\left(t^{{\alpha }_2-1}{E}_{({\alpha }_2-{\alpha }_1,{\alpha }_2),{\alpha }_2}\left(-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1},-\left(-{\gamma }_2{\Delta }_x+{\gamma }_3\right)t^{{\alpha }_2}\right)\delta \left(x\right)\right){*}_{x,t}q(x,t).\end{array}$$

Space-time fundamental solution kernels.

To derive explicit representations of the kernels generating the homogeneous solution, we define

$$\begin{array}{c}\hfill G_{h,j}(x,t)={\displaystyle \sum _{l_1,l_2=0}^{\infty }}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{(-{\mathcal{A}}_1\left(t\right))}^{l_1}{\left({\mathcal{A}}_2\left(t\right){\Delta }_x\right)}^{l_2}\delta \left(x\right){\Psi }_j\left(t\right),j=0,1,\end{array}$$

(61)

so that

$$\begin{array}{c}\hfill u_h(x,t)={\displaystyle \sum _{j=0}^1}{G}_{h,j}(·,t){*}_x{g}_{j+1}\left(x\right).\end{array}$$

Each kernel $G_{h,j}$ is associated with the role of a fundamental solution associated with the initial condition $g_{j+1}$. We now derive an explicit representation of $G_{h,j}$. Expanding ${\mathcal{A}}_1^{l_1}\left(t\right)$ via the binomial theorem and applying the semigroup property of Riemann-Liouville fractional integrals, we obtain

$$\begin{array}{cc}\hfill {\mathcal{A}}_1^{l_1}\left(t\right){\mathcal{A}}_2^{l_2}\left(t\right)& ={\displaystyle \sum _{l_3=0}^{l_1}}\left({\displaystyle \genfrac{}{}{0pt}{}{l_1}{l_3}}\right){\gamma }_1^{l_3}{\gamma }_3^{l_1-l_3}{I}_{t,0^{+}}^{l_3({\alpha }_2-{\alpha }_1)}{I}_{t,0^{+}}^{(l_1-l_3){\alpha }_2}{\gamma }_2^{l_2}{I}_{t,0^{+}}^{l_2{\alpha }_2}\hfill \\ \hfill & ={\displaystyle \sum _{l_3=0}^{l_1}}\left({\displaystyle \genfrac{}{}{0pt}{}{l_1}{l_3}}\right){\gamma }_1^{l_3}{\gamma }_2^{l_2}{\gamma }_3^{l_1-l_3}{I}_{t,0^{+}}^{(l_1+l_2){\alpha }_2-l_3{\alpha }_1}.\hfill \end{array}$$

(62)

Now, its action on ${\Psi }_j$ gives

$$\begin{array}{c}\hfill {\mathcal{A}}_1^{l_1}\left(t\right){\mathcal{A}}_2^{l_2}\left(t\right){\Psi }_j\left(t\right)={\displaystyle \sum _{l_3=0}^{l_1}}\left({\displaystyle \genfrac{}{}{0pt}{}{l_1}{l_3}}\right){\gamma }_1^{l_3}{\gamma }_2^{l_2}{\gamma }_3^{l_1-l_3}{\displaystyle \frac{t^{{\theta }_j+(l_1+l_2){\alpha }_2-l_3{\alpha }_1}}{\Gamma \left(1+{\theta }_j+(l_1+l_2){\alpha }_2-l_3{\alpha }_1\right)}}.\end{array}$$

Putting this result back into (61), we obtain:

$$\begin{array}{cc}\hfill & G_{h,j}(x,t)\hfill \\ \hfill & ={\displaystyle \sum _{l_1,l_2=0}^{\infty }}{(-1)}^{l_1}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \sum _{l_3=0}^{l_1}}\left({\displaystyle \genfrac{}{}{0pt}{}{l_1}{l_3}}\right){\gamma }_1^{l_3}{\gamma }_2^{l_2}{\gamma }_3^{l_1-l_3}{\displaystyle \frac{t^{{\theta }_j+(l_1+l_2){\alpha }_2-l_3{\alpha }_1}}{\Gamma \left(1+{\theta }_j+(l_1+l_2){\alpha }_2-l_3{\alpha }_1\right)}}{\Delta }_x^{l_2}\delta \left(x\right).\hfill \end{array}$$

(63)

By absolute convergence of the series, we may reindex the series by setting $l_1=l_3+l_4$. Then the indices $l_3$ and $l_4$ run independently from 0 to ∞ and we get

$$\begin{array}{cc}\hfill & G_{h,j}(x,t)\hfill \\ \hfill & ={\displaystyle \sum _{l_2,l_3,l_4=0}^{\infty }}{(-1)}^{l_3+l_4}{\displaystyle \frac{(l_3+l_4+l_2)!}{l_2!l_3!l_4!}}{\gamma }_1^{l_3}{\gamma }_2^{l_2}{\gamma }_3^{l_4}{\displaystyle \frac{t^{{\theta }_j+(l_2+l_3+l_4){\alpha }_2-l_3{\alpha }_1}}{\Gamma \left(1+{\theta }_j+(l_2+l_3+l_4){\alpha }_2-l_3{\alpha }_1\right)}}{\Delta }_x^{l_2}\delta \left(x\right).\hfill \end{array}$$

While equation (59) establishes the fundamental solution from an algebraic operator perspective, an explicit geometric mapping in the real spatial domain ${\mathbb{R}}^n$ requires evaluating the direct action of the powers of the Laplacian ${\left({\Delta }_x\right)}^{l_2}$ on the Dirac delta distribution. In the framework of Riesz distributions (see [10]), one has

$$\begin{array}{c}\hfill {\left({\Delta }_x\right)}^{l_2}\delta \left(x\right)={(-1)}^{l_2}{c}_{n,l_2}{\left|x\right|}^{-n-2l_2},\end{array}$$

(64)

where

$$\begin{array}{c}\hfill c_{n,l_2}={\displaystyle \frac{2^{2l_2}\Gamma \left(\frac{n}{2}+l_2\right)}{{\pi }^{\frac{n}{2}}\Gamma \left(-l_2\right)}}.\end{array}$$

(65)

Here, ${\left|x\right|}^{-n-2l_2}$ is understood as the Riesz distribution obtained by analytic continuation of the family ${\left|x\right|}^{-n-2s}$, $s\in \mathbb{C}$. Substituting (64)-(65) into (63), we obtain the following series representation in the real space domain:

$$\begin{array}{c}\hfill G_{h,j}(x,t)={\displaystyle \frac{t^{{\theta }_j}}{{\pi }^{\frac{n}{2}}{\left|x\right|}^n}}{\displaystyle \sum _{l_2,l_3,l_4=0}^{\infty }}{\displaystyle \frac{(l_2+l_3+l_4)!}{l_2!l_3!l_4!}}{\displaystyle \frac{\Gamma \left(\frac{n}{2}+l_2\right)}{\Gamma \left(-l_2\right)}}{\displaystyle \frac{z_1^{l_3}{z}_2^{l_2}{z}_3^{l_4}}{\Gamma \left(1+{\theta }_j+l_3({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2+l_4{\alpha }_2\right)}},\end{array}$$

(66)

where the complex variables are defined as $z_1=-{\gamma }_1{t}^{{\alpha }_2-{\alpha }_1},$$z_2=-4{\gamma }_2{t}^{{\alpha }_2}{\left|x\right|}^{-2}$, and $z_3=-{\gamma }_3{t}^{{\alpha }_2}$.

At first glance, the triple series representation (66) presents an analytical difficulty, as the term ${\displaystyle \frac{1}{\Gamma \left(-l_2\right)}}$ formally vanishes for all non-negative integer indices $l_2\in {\mathbb{N}}_0$. Consequently, this series cannot be interpreted in a naive, term-by-term manner, since the solution is nontrivial. As we shall establish below, a rigorous regularization of this representation is achieved via its analytic continuation into the complex plane. By reformulating the series as a Mellin-Barnes contour integral, the problematic term ${\displaystyle \frac{1}{\Gamma \left(-s_2\right)}}$ in the denominator is precisely canceled by an identical Gamma factor in the numerator, thereby restoring the non-trivial physical and spatial structure of the fundamental solution.

For the particular solution, using (39) and (62) in (60), we obtain the following series representation of the corresponding fundamental solution:

$$\begin{array}{c}\hfill G_p(x,t)={\displaystyle \sum _{l_1,l_2=0}^{\infty }}{(-1)}^{l_1}{\displaystyle \frac{(l_1+l_2)!}{l_1!l_2!}}{\displaystyle \sum _{l_3=0}^{l_1}}\left({\displaystyle \genfrac{}{}{0pt}{}{l_1}{l_3}}\right){\gamma }_1^{l_3}{\gamma }_2^{l_2}{\gamma }_3^{l_1-l_3}{\displaystyle \frac{t^{{\alpha }_2-1+(l_1+l_2){\alpha }_2-l_3{\alpha }_1}}{\Gamma \left({\alpha }_2+(l_1+l_2){\alpha }_2-l_3{\alpha }_1\right)}}{\Delta }_x^{l_2}\delta \left(x\right).\end{array}$$

Now, using (64) and reindexing the series by setting $l_1=l_3+l_4$, we get

$$\begin{array}{c}\hfill G_p(x,t)={\displaystyle \frac{t^{{\alpha }_2-1}}{{\pi }^{\frac{n}{2}}{\left|x\right|}^n}}{\displaystyle \sum _{l_2,l_3,l_4=0}^{\infty }}{\displaystyle \frac{(l_2+l_3+l_4)!}{l_2!l_3!l_4!}}{\displaystyle \frac{\Gamma \left(\frac{n}{2}+l_2\right)}{\Gamma \left(-l_2\right)}}{\displaystyle \frac{z_1^{l_3}{z}_2^{l_2}{z}_3^{l_4}}{\Gamma \left({\alpha }_2+l_3({\alpha }_2-{\alpha }_1)+l_2{\alpha }_2+l_4{\alpha }_2\right)}}.\end{array}$$

(67)

Mellin-Barnes representations

Using Mellin-Barnes representations together with the residue theorem, we can rewrite the triple series defining $G_{h,j}(x,t)$ as the following Mellin-Barnes contour integral:

$$\begin{array}{cc}\hfill G_{h,j}(x,t)& ={\displaystyle \frac{t^{{\theta }_j}}{{\pi }^{\frac{n}{2}}{\left|x\right|}^n}}{\displaystyle \frac{1}{{\left(2\pi i\right)}^3}}{\int }_{{\mathcal{L}}_{s_3}}{\int }_{{\mathcal{L}}_{s_2}}{\int }_{{\mathcal{L}}_{s_1}}{\displaystyle \frac{\Gamma \left(1+s_1+s_2+s_3\right)\Gamma \left(-s_1\right)\Gamma \left(-s_2\right)\Gamma \left(\frac{n}{2}+s_2\right)\Gamma \left(-s_3\right)}{\Gamma \left(1+{\theta }_j+s_1({\alpha }_2-{\alpha }_1)+s_2{\alpha }_2+s_3{\alpha }_2\right)\Gamma \left(-s_2\right)}}\hfill \\ \hfill & \times z_1^{s_1}{z}_2^{s_2}{z}_3^{s_3}d{s}_{1}d{s}_{2}d{s}_3.\hfill \end{array}$$

Upon cancellation of the common Gamma factors in the integrand, the Mellin-Barnes representation reduces to

$$\begin{array}{cc}\hfill G_{h,j}(x,t)& ={\displaystyle \frac{t^{{\theta }_j}}{{\pi }^{\frac{n}{2}}{\left|x\right|}^n}}{\displaystyle \frac{1}{{\left(2\pi i\right)}^3}}{\int }_{{\mathcal{L}}_{s_3}}{\int }_{{\mathcal{L}}_{s_2}}{\int }_{{\mathcal{L}}_{s_1}}{\displaystyle \frac{\Gamma \left(1+s_1+s_2+s_3\right)\Gamma \left(-s_1\right)\Gamma \left(\frac{n}{2}+s_2\right)\Gamma \left(-s_3\right)}{\Gamma \left(1+{\theta }_j+s_1({\alpha }_2-{\alpha }_1)+s_2{\alpha }_2+s_3{\alpha }_2\right)}}\hfill \\ \hfill & \times z_1^{s_1}{z}_2^{s_2}{z}_3^{s_3}d{s}_{1}d{s}_{2}d{s}_3.\hfill \end{array}$$

The resulting Mellin-Barnes integral provides the analytic continuation of the original triple series. According to the definition of multivariate Fox H-functions (see Appendix C in [23]), the above Mellin-Barnes triple integral can be identified with the trivariate Fox H-function

$$\begin{array}{c}\hfill G_{h,j}(x,t)={\displaystyle \frac{t^{{\theta }_j}}{{\pi }^{n/2}{\left|x\right|}^n}}{H}_{1,1;0,1;1,0;0,1}^{0,1;1,0;0,1;1,0}\left[\begin{array}{cc}\begin{array}{c}{\displaystyle z_1}\\ {\displaystyle z_2}\\ {\displaystyle z_3}\end{array}& \begin{array}{c}{\displaystyle (0;1,1,1);0.5ex0.4cm0.4pt;\left(1-\frac{n}{2},1\right);0.5ex0.4cm0.4pt}\\ {\displaystyle (-{\theta }_j;{\alpha }_2-{\alpha }_1,{\alpha }_2,{\alpha }_2);(0,1);0.5ex0.4cm0.4pt;(0,1)}\end{array}\end{array}\right],\end{array}$$

(68)

under the standard admissibility and convergence conditions for multivariate Fox H-functions.

Proceeding analogously, we obtain the following Mellin-Barnes contour integral representation for $G_p$:

$$\begin{array}{cc}\hfill G_p(x,t)& ={\displaystyle \frac{t^{{\alpha }_2-1}}{{\pi }^{\frac{n}{2}}{\left|x\right|}^n}}{\displaystyle \frac{1}{{\left(2\pi i\right)}^3}}{\int }_{{\mathcal{L}}_{s_3}}{\int }_{{\mathcal{L}}_{s_2}}{\int }_{{\mathcal{L}}_{s_1}}{\displaystyle \frac{\Gamma \left(1+s_1+s_2+s_3\right)\Gamma \left(-s_1\right)\Gamma \left(\frac{n}{2}+s_2\right)\Gamma \left(-s_3\right)}{\Gamma \left({\alpha }_2+s_1({\alpha }_2-{\alpha }_1)+s_2{\alpha }_2+s_3{\alpha }_2\right)}}\hfill \\ \hfill & \times z_1^{s_1}{z}_2^{s_2}{z}_3^{s_3}d{s}_{1}d{s}_{2}d{s}_3\hfill \end{array}$$

(69)

which can be identified with the following trivariate Fox H-function

$$\begin{array}{c}\hfill G_p(x,t)={\displaystyle \frac{t^{{\alpha }_2-1}}{{\pi }^{\frac{n}{2}}{\left|x\right|}^n}}{H}_{1,1;0,1;1,0;0,1}^{0,1;1,0;0,1;1,0}\left[\begin{array}{cc}\begin{array}{c}{\displaystyle z_1}\\ {\displaystyle z_2}\\ {\displaystyle z_3}\end{array}& \begin{array}{c}{\displaystyle (0;1,1,1);0.5ex0.4cm0.4pt;\left(1-\frac{n}{2},1\right);0.5ex0.4cm0.4pt}\\ {\displaystyle (1-{\alpha }_2;{\alpha }_2-{\alpha }_1,{\alpha }_2,{\alpha }_2);(0,1);0.5ex0.4cm0.4pt;(0,1)}\end{array}\end{array}\right].\end{array}$$

Fourier-Laplace verification

To verify the consistency between the operational representation and the explicit space-time representation, we apply the joint Fourier-Laplace transform. We apply the temporal Laplace transform ($t\to s$) combined with the n-dimensional spatial Fourier transform ($x\to \xi$). In the Fourier domain, the spatial transform diagonalizes the differential structure of the Laplacian, converting the complex distributional Riesz kernel from (66) into a simple scalar multiplier via the standard relation $\mathcal{F}\left\{{\left({\Delta }_x\right)}^{l_2}\delta \left(x\right)\right\}={(-1)}^{l_2}{\left|\xi \right|}^{2l_2}$. The Laplace transform converts the fractional time operators into algebraic powers of s. Substituting the transformed expressions into either representation yields the following algebraic series representation in the transformed domain:

$$\begin{array}{c}\hfill {\tilde{\widehat{G}}}_{h,j}(\xi ,s)=s^{-{\theta }_j-1}{\displaystyle \sum _{l_2,l_3,l_4=0}^{\infty }}{\displaystyle \frac{(l_2+l_3+l_4)!}{l_2!l_3!l_4!}}{\left(-{\gamma }_1{s}^{{\alpha }_1-{\alpha }_2}\right)}^{l_3}{\left(-{\gamma }_2{\left|\xi \right|}^2{s}^{-{\alpha }_2}\right)}^{l_2}{\left(-{\gamma }_3{s}^{-{\alpha }_2}\right)}^{l_4}.\end{array}$$

Applying the generalized binomial expansion, we obtain:

$$\begin{array}{cc}\hfill {\tilde{\widehat{G}}}_{h,j}(\xi ,s)& =s^{-{\theta }_j-1}\left({\displaystyle \frac{1}{1-\left(-{\gamma }_1{s}^{{\alpha }_1-{\alpha }_2}-{\gamma }_2{\left|\xi \right|}^2{s}^{-{\alpha }_2}-{\gamma }_3{s}^{-{\alpha }_2}\right)}}\right)\hfill \\ \hfill & ={\displaystyle \frac{s^{-{\theta }_j-1}}{1+{\gamma }_1{s}^{{\alpha }_1-{\alpha }_2}+{\gamma }_2{\left|\xi \right|}^2{s}^{-{\alpha }_2}+{\gamma }_3{s}^{-{\alpha }_2}}}.\hfill \end{array}$$

Finally, multiplying both the numerator and the denominator by $s^{{\alpha }_2}$ normalizes the fraction and eliminates the relative negative fractional exponents, yielding the characteristic symbol:

$$\begin{array}{c}\hfill {\tilde{\widehat{G}}}_{h,j}(\xi ,s)={\displaystyle \frac{s^{{\alpha }_2-{\theta }_j-1}}{s^{{\alpha }_2}+{\gamma }_1{s}^{{\alpha }_1}+{\gamma }_2{\left|\xi \right|}^2+{\gamma }_3}}.\end{array}$$

(70)

The denominator of (70) coincides with the characteristic symbol associated with equation (12). Hence, the operational representation (59) and the explicit space-time representation (66) determine the same Fourier-Laplace symbol and are therefore equivalent. An analogous verification applies to the particular solution. Combining the homogeneous and particular components, we obtain the following result.

Theorem 3.

The solution of the time-fractional telegraph equation (12) with constant coefficients (32) and subject to the initial conditions (15)-(13) is given by $u(x,t)=u_h(x,t)+u_p(x,t)$ where

$$\begin{array}{c}\hfill u_h(x,t)={\displaystyle \sum _{j=0}^1}{\int }_{{\mathbb{R}}^n}{G}_{h,j}(x-z,t)g_{j+1}\left(z\right)dz\end{array}$$

(71)

and

$$\begin{array}{c}\hfill u_p(x,t)={\int }_{{\mathbb{R}}^n}{\int }_0^t{G}_p(x-z,t-w)q(z,w)dwdz,\end{array}$$

(72)

where the fundamental kernels $G_{h,j}$ and $G_p$ are given respectively by the trivariate Fox H-functions (68) and (69).

Remark 4.

If we consider ${\gamma }_3=0$ in (66) and in (67), then the triple series reduce to a double series. Consequently, the trivariate Fox-H functions in (71) and (72) reduce to bivariate Fox-H functions of two variables, allowing us to recover the second, third, and fourth terms in [26] with $c_2=1$, $c_1={\gamma }_1$, $c_0^2={\gamma }_2$, $d=0$, and $\psi \left(t\right)=t$.

4. Conclusions and Final Remark

In this paper, we have derived explicit representation formulas for the solution of a time-fractional telegraph equation with time-dependent coefficients within the framework of the Hilfer fractional derivative. By formulating the problem as a Cauchy problem in the Fourier domain, the general space-time solution was obtained in a distributional sense as convolutions of convergent Neumann-type series involving iterated compositions of non-commuting operators. For the specific cases of constant and power-law coefficients, the precise structure of the solutions was explicitly derived. Furthermore, in the constant-coefficient regime, our formulas successfully recover established results from the literature, validating the consistency and correctness of our approach.

The results obtained in this paper can be naturally extended to the $\psi$-fractional setting, where the temporal operators are defined with respect to a strictly increasing and positive scale function $\psi$. For completeness, we briefly recall the foundational definitions required for this framework.

Definition 3.(cf. [22]) Let $[a,b]$ be a finite or infinite interval on the real line $\mathbb{R}$ and $\alpha >0$. Also, let ψ be an increasing and positive monotone function on $(a,b)$. The left Riemann-Liouville fractional integral of a function f with respect to another function ψ on $[a,b]$ is defined by

$$\begin{array}{c}\hfill \left(I_{a^{+}}^{\alpha ;\psi }f\right)\left(t\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_a^t{\psi }^{\prime }\left(w\right){(\psi \left(t\right)-\psi \left(w\right))}^{\alpha -1}f\left(w\right)dw,t>a.\end{array}$$

Definition 4.(cf. [22]) Let $\alpha >0$, $m=\lfloor \alpha \rfloor +1$, $I=[a,b]$ be a finite or infinite interval on the real line and $f,\psi \in C^m[a,b]$ two functions such that ψ is a positive monotone increasing function and ${\psi }^{\prime }\left(t\right)\ne 0$, for all $t\in I$. The ψ-Hilfer left fractional derivative $H{\mathbb{D}}_{t,0^{+}}^{\alpha ,\mu ;\psi }$ of order α and type $\mu \in [0,1]$ is defined by

$$\begin{array}{c}\hfill \left({{}^{H}D}_{a^{+}}^{\alpha ,\mu ;\psi }f\right)\left(t\right)=I_{a^{+}}^{\mu (m-\alpha );\psi }{\left({\displaystyle \frac{1}{{\psi }^{\prime }\left(t\right)}}{\displaystyle \frac{d}{dt}}\right)}^m{I}_{a^{+}}^{(1-\mu )(m-\alpha );\psi }f\left(t\right).\end{array}$$

By applying the transmutation method of $\psi$-fractional calculus, as established in [3,8], the solution to the corresponding $\psi$-Hilfer problem can be structurally mapped to $u(x,\psi (t\left)\right)$. This formulation implies that the generalized solutions are retrieved by evaluating the classical time-fractional solutions along the deterministic time-change governed by the scale function $\psi \left(t\right)$, which dictates the kernel of the underlying fractional operators.