1 Introduction

In recent paper by Lv et al. [1] authors studied the question of stochastic representation of anomalous diffusion with space and time dependent drift and diffusion coefficient. This problem was also studied in [2–6], where the author investigated solution of the anomalous diffusion equation with space–time dependent forces. In this paper we show that the result of [1] is true in even more general setting. Namely, the factorization of the space–time dependent drift and diffusion coefficients into a space dependent part and a time dependent part is not necessary. We also point out some drawbacks in the derivation of the main result in [1] and present the correct proof.

In [1] authors analyzed the following fractional Fokker–Planck equation (FFPE)

$$\begin{aligned} \frac{\partial }{\partial t}p(x,t)=\Bigg \{ -\frac{\partial }{\partial x}F(x,t)+\frac{1}{2}\frac{\partial ^2}{\partial x^2} D(x,t)\Bigg \}{_0}{ D_t^{1-\alpha }}p(x,t),\quad 0<\alpha < 1, \end{aligned}$$
(1)

with initial condition \(p(x,0)=\delta (x)\). The operator \({_0}{ D_t^{1-\alpha }}\) is the fractional Riemann-Liouville derivative which is defined as

$$\begin{aligned} {_0}{ D_t^{1-\alpha }} q(t)=\frac{1}{\Gamma (\alpha )}\frac{\partial }{\partial t}\int \limits _0^t (t-s)^{\alpha -1}q(s) ds. \end{aligned}$$

Moreover, the operator in braces is the classical Fokker–Planck operator with \(F(x,t)\) and \(D(x,t)\) being the drift and diffusion parameters, respectively.

Equation (1) describes the evolution in time of the the probability density function (PDF) \(p(x,t)\) of some subdiffusive process with sublinear in time mean square displacement. The operator \({_0}{ D_t^{1-\alpha }}\) introduces a convolution integral with a slowly decaying power-law kernel, which is typical for memory effects in complex systems. Appearance of \({_0}{ D_t^{1-\alpha }}\) in Eq. (1) corresponds to the trapping events in the underlying motion of the test particle. These trapping events are characterictic for subdiffusive dynamics and result in sublinear in time mean square displacement [7].

Equation (1) was first derived in [8] in the framework of continuous-time random walk (CTRW) with heavy-tailed waiting times. In [8] the authors concentrated on the case of space-dependent drift \(F(x,t)=F(x)\) and constant diffusion \(D(x,t)=const\). Since then FFPE became the standard physical equation describing subdiffusive dynamics [7].

Equation (1) for the case of time-dependent force \(F(x,t)=F(t)\) was first derived in [2]. The derivation was based on the generalized master equation with two balance conditions: the probability conservation in a given state and under transition between different states. The main difference from the previous space-dependent case was lying in the fact that the fractional operator \({_0}{ D_t^{1-\alpha }}\) appeared to the right of the Fokker–Planck operator, so that it did not act on the time-dependent force. As a result the force was not modified, which fulfilled the physical requirement that the external time-dependent force cannot be influenced by the environment.

The case of space–time dependent coefficients was introduced and derived in [3, 6] using Langevin and CTRW approaches, which finally clarified issues addressed in [9] for space–time-dependent driving.

In [1], the authors assume that the drift \(F\) and diffusion \(D\) have the following, very special, factorized form \(F(x,t)=F(x)f(t)\) and \(D(x,t)=D(x)\tilde{d}(t)\). The main purpose of [1] was to prove that Eq. (1) describes PDF of the following stochastic process

$$\begin{aligned} X(0)&= 0\nonumber \\ dX(t)&= F\big (X(t)\big )f(t)dS_{\alpha }(t)+\sqrt{D\big (X(t)\big )\tilde{d}(t)}dB(S_{\alpha }(t)), \end{aligned}$$
(2)

where \(F(x), D(x)\in C^{\infty }((-\infty , \infty ))\), \(f(t), \tilde{d}(t)\in C([0, \infty ))\) and \(S_\alpha (t)\) is the inverse \(\alpha \)-stable subordinator [4]

$$\begin{aligned} S_{\alpha }(t)=\inf \{\tau >0: U_{\alpha }(\tau )>t \}. \end{aligned}$$

Here, \(U_{\alpha }\) is the \(\alpha \)-stable subordinator [10, 11] with Laplace transform \({E}\left[ e^{-kU_{\alpha }(\tau )} \right] =e^{-\tau k^\alpha }\), \(0<\alpha <1\). Moreover, \(C^{\infty }((-\infty , \infty ))\) denotes the space of functions with well defined derivatives of any order on \((-\infty , \infty )\),whereas \(C([0, \infty ))\) is the space of continuous functions on \([0, \infty )\).

The process \(X(t)\) in (2) can be equivalently represented in the following subordinated form

$$\begin{aligned} X(t)=Y(S_\alpha (t)), \end{aligned}$$

where \(Y(t)\) is given by

$$\begin{aligned}&Y(0)=0\\&dY(t)=F(Y(t))f(U_\alpha (t))dt+\sqrt{D(Y(t))\tilde{d}(U_\alpha (t))}dB(t),\quad t\ge 0, \end{aligned}$$

Let us now explain the structure and physical meaning of the process \(X(t)=Y(S_\alpha (t))\). Note first that subordination is a change of time. The subordinator \(S_\alpha (t)\) becomes a new operational time of the system and it introduces the trapping events to the motion of the particle (see [3]). Moreover, changing the time does not affect the spatial coefficients \(F(x)\) and \(D(x)\), it can only affect the temporal ones. Note that the parent process \(Y(t)\) is actually a diffusion process, except from the fact that \(U_\alpha (t)\) appears in the temporal parts of its drift and diffusion coefficients. Because of this, after subordination the temporal parts of the drift and diffusion coefficients are given by \(f(U_\alpha (S_\alpha (t)))\) and \(\tilde{d}(U_\alpha (S_\alpha (t)))\). Now, since in every moment of jump we have \(U_\alpha (S_\alpha (t))=t\) (see [3]), the time-dependent biases acting on the particle are equal to \(f(U_\alpha (S_\alpha (t)))=f(t)\) and \(\tilde{d}(U_\alpha (S_\alpha (t)))=\tilde{d}(t)\). This means that whenever the particle does not rest in the trap, the coefficients \(f(t)\) and \(\tilde{d}(t)\) are not modified by the subordinator \(S_\alpha (t)\). This fulfills the physical requirement that the external time-dependent force cannot be influenced by the environment and corresponds to the fact that the fractional operator \({_0}{ D_t^{1-\alpha }}\) in (1) appears to the right of the Fokker–Planck operator. Summing up, the process \(X(t)=Y(S_\alpha (t))\) describes correctly the subdiffusive motion of the particle driven by space–time dependent coefficients \(F(x,t)\) and \(D(x,t)\).

Although the main result presented in [1] is correct, its derivation is wrong. In this paper we point out problems occurring in the proof of this result. Next, we present the correct proof in even more general case, when both drift and diffusion can have any, not necessarily factorized, form.

The paper is organized as follows. In Sect. 2 we point out critical errors in the proof presented in [1]. In Sect. 3 we present a correct proof of the main result of [1] in general setting. Applying the method from [6], we derive the stochastic representation of FFPE with arbitrary drift and diffusion parameters.

2 Problem in the Proof

Let us look at the second part of the proof of Theorem 2 in [1], where the authors consider the case \(f(t)\ne \tilde{d}(t)\). First of all, the statement that \(U_\alpha (\tau )\in \mathcal {F}_{S_\alpha (t)}\) is incorrect. Note that \(\mathcal {F}_{S_\alpha (t)}\) is the \(\sigma \)-algebra generated by \(S_\alpha (t)\) (the family of sets containing all the information about \(S_\alpha (t)\)). Since both sets \(\{U_\alpha (\tau )\ge t \}\) and \(\{ S_\alpha (t)\le \tau \}\) are equal [4, 12], so the \(\sigma \)-algebra generated by \(U_\alpha (\tau )\) must contain all the sets of the form \(\{ S_\alpha (t_0)\le \tau \}\), where \(t_0\ge 0\). Therefore, \(U_\alpha (\tau )\) is not \(\mathcal {F}_{S_\alpha (t)}\)-measurable. It is rather measurable with respect to the \(\sigma \)-algebra generated by the whole process \(S_\alpha \). Indeed, to know the value of \(U_\alpha (\tau )\) we need to know the whole trajectory of the process \(S_\alpha \) and not only its one particular value \(S_\alpha (t)\) at the single time point \(t\).

Moreover, the authors incorrectly assumed that both functions \(G(\tau )=F(Y(\tau ))e^{ikY(\tau )}\) and \(g(\tau )=D(Y(\tau ))e^{ikY(\tau )}\) are known integral functions on the finite interval \([0,t]\), if we condition on \(\sigma \)-algebra \(\mathcal {F}_{S_\alpha (t)}\). As one can see, process \(Y\) defined by Eq. (20) in [1] has another source of randomness apart of \(U_\alpha (\tau )\). Namely Brownian motion \(B(\tau )\), which is independent from \(U_\alpha (\tau )\). So \(Y\) is still a random process, even if we condition on \(S_{\alpha }(t)\). As a consequence \(G(\tau )\) and \(g(\tau )\) are complex random processes, so Corollary 1 in [1] does not apply. Recall that Corollary 1 applies only to deterministic functions, and cannot be easily extended to random processes. Thus the crucial step from Eq.

$$\begin{aligned} \frac{\partial }{\partial t}E\left[ ik \int \limits _0^tf(t_1)F(Y(S_\alpha (t_1)))e^{ikY(S_\alpha (t_1))}dS_\alpha (t_1)\right] , \end{aligned}$$
(3)

to

$$\begin{aligned} ikf(t){_0}{ D_t^{1-\alpha }}E\left[ F(X(t))e^{ikX(t)}\right] , \end{aligned}$$
(4)

in [1], which is based on Corollary 1, is unjustified. As a consequence, the whole proof of the authors is incorrect.

3 Stochastic Representation of Generalized FFPE

In this section we give a correct proof of the stochastic representation of FFPE (1) presented in [1]. Applying the methods presented in [6], we obtain the result in even more general setting, i.e. when the drift and diffusion are arbitrary functions. In particular, if we choose \(D(x,t)\) and \(F(x,t)\) to have factorized form, we obtain the statement from [1].

Theorem 1

Let \(S_\alpha (t)\) be the inverse \(\alpha \)-stable subordinator, independent of the standard Brownian motion \(B(\tau )\). Let \(Y(t)\) be the solution of the stochastic differential equation

$$\begin{aligned}&Y(0)=0\nonumber \\&dY(t)=F(Y(t),U_\alpha (t))dt+\sqrt{D(Y(t),U_\alpha (t))}dB(t),\quad t\ge 0, \end{aligned}$$
(5)

where the functions \(D(x,t)\ge 0\) and \(F(x,t)\) are twice differentiable with respect to both variables and satisfy the usual Lipschitz condition. Then, the PDF of the process

$$\begin{aligned} X(t)=Y(S_\alpha (t)) \end{aligned}$$

is the solution of the fractional Fokker–Planck equation

$$\begin{aligned} \frac{\partial }{\partial t} p(x,t)=\left\{ -\frac{\partial }{\partial x} F(x,t) + \frac{1}{2} \frac{\partial ^2}{\partial x^2} D(x,t) \right\} {_0}{D_t^{1-\alpha }}p(x,t) ,\quad t\ge 0 \end{aligned}$$
(6)

with the initial condition \(p(x,0)=\delta (x)\).

Proof

We can rewrite Eq. (5) as a system of stochastic equations:

$$\begin{aligned} \left( \begin{array}{l} dY(t)\\ dZ(t)\\ \end{array}\right) =\left( \begin{array}{l} F(Y(t),Z(t))\\ 0\\ \end{array}\right) dt+ \left( \begin{array}{c} \sqrt{D(Y(t),Z(t))} dB_t\\ dU_\alpha (t)\\ \end{array}\right) . \end{aligned}$$
(7)

Equation (7) belong to the general class of stochastic processes driven by Lévy noise, thus (see [13]) the infinitesimal generator for the process (7) has the form

$$\begin{aligned}&Af(y,z)= F(y,z)\frac{\partial }{\partial y} f(y,z)+ \frac{1}{2}D(y,z)\frac{\partial ^2}{\partial y^2} f(y,z) \nonumber \\&\quad + \int \limits _0^\infty \left[ f(y,z+z')-f(y,z) \right] \frac{\alpha }{\Gamma (1-\alpha )}{z'}^{-1-\alpha }dz'. \end{aligned}$$
(8)

The infinitesimal generator of the stochastic process \((Y(t),Z(t))\) is a linear operator defined as

$$\begin{aligned} Af(y,z)=\lim _{t \rightarrow 0} \frac{E\Big [f\Big (\left( Y(t),Z(t)\right) +(y,z)\Big )-f(y,z)\Big ]}{t}. \end{aligned}$$

The set of functions \(f:\mathbf {R}^2 \rightarrow \mathbf {R}\) for which the above limit exists for all \((y,z)\in \mathbf {R}^2\) is the domain of \(A\). The generator carries a lot of information about the process. It is an important and useful tool in analyzing evolution equations. Let \(q_t(y,z)\) denote the PDF of the joint distribution \((Y(t),Z(t))\). Then the standard Fokker–Planck equation for \(q_t(y,z)\) can be written as (see [13])

$$\begin{aligned} \frac{\partial }{\partial t}q_t(y,z)=A^+q_t(y,z), \end{aligned}$$
(9)

where \(A^+\) is the formal adjoint of A. It is the operator that satisfies \(\langle Ag,h\rangle = \langle g,A^+h\rangle \) where \(\langle \cdot ,\cdot \rangle \) is the inner product in \(\mathcal {L}^2(\mathbf {R}^2)\) and \(g\), \(h\) are sufficiently smooth functions. In particular we have

$$\begin{aligned} \int \limits _{\mathbf {R}^2} Af(y,z)q_t(y,z)dydz= \int \limits _{\mathbf {R}^2} f(y,z)A^+q_t(y,z)dydz , \end{aligned}$$
(10)

for all \(f \in \mathbf {C}_c^\infty (\mathbf {R}^2)\)—infinitely differentiable functions from \(\mathbf {R}^2\) to \(\mathbf {R}\) with compact support, also known as bump functions. Using the above relation we will find the operator \(A^+\). The right hand side of Eq. (10) yields

$$\begin{aligned} \int \limits _{\mathbf {R}^2} Af(y,z) q_t(y,z)dydz&= \int \limits _{\mathbf {R}^2} \Bigg (F(y,z)\frac{\partial }{\partial y} f(y,z) + \frac{1}{2}D(y,z)\frac{\partial ^2}{\partial y^2} f(y,z) \\&+ \int \limits _0^\infty \left[ f(y,z+z')\!-\!f(y,z)\! \right] \! \frac{\alpha }{\Gamma (1\!-\!\alpha )}{z'}^{-1-\alpha } d{z'}\!\Bigg )q_t(y,z)dydz . \end{aligned}$$

We will deal with each summand separately. Integration by parts and the fact that the support of \(f\) is compact yield

$$\begin{aligned} \int \limits _{\mathbf {R}^2} F(y,z)\frac{\partial }{\partial y} f(y,z)q_t(y,z)dydz=-\int \limits _{\mathbf {R}^2} \frac{\partial }{\partial y}\left( F(y,z)q_t(y,z)\right) f(y,z)dydz. \end{aligned}$$
(11)

Similarly we use integration by parts two times to get

$$\begin{aligned}&\int \limits _{\mathbf {R}^2} \frac{1}{2}D(y,z)\frac{\partial ^2}{\partial y^2} f(y,z) q_t(y,z)d(y,z)\nonumber \\&\quad = \int \limits _{\mathbf {R}^2} \frac{\partial ^2}{\partial y^2} \left( \frac{1}{2}D(y,z)q_t(y,z)\right) f(y,z) dydz. \end{aligned}$$
(12)

For the last summand let us observe that

$$\begin{aligned}&\int \limits _0^\infty \left[ f(y,z+z')-f(y,z) \right] \frac{\alpha }{\Gamma (1-\alpha )}{z'}^{-1-\alpha }dz' = \\&\quad -\left[ f(y,z+z')-f(y,z) \right] \frac{{z'}^{-\alpha }}{\Gamma (1-\alpha )}\Big |_{z'=0}^\infty +\int \limits _0^\infty \frac{\partial }{\partial z'} f(y,z+z') \frac{{z'}^{-\alpha }}{\Gamma (1-\alpha )}dz', \end{aligned}$$

where

$$\begin{aligned} \left[ f(y,z+z')-f(y,z) \right]&\frac{{z'}^{-\alpha }}{\Gamma (1-\alpha )}\Big |_{z'=0}^\infty \\&= \lim _{z' \rightarrow 0} \frac{f(y,z+z')-f(y,z)}{z'} \frac{{z'}^{1-\alpha }}{\Gamma (1-\alpha )} =0, \end{aligned}$$

since \(f\) is differentiable. After substitution \(z+z'=w\) we obtain:

$$\begin{aligned} \int \limits _0^\infty \frac{\partial }{\partial z'} f(y,z+z') \frac{{z'}^{-\alpha }}{\Gamma (1-\alpha )}dz'=\int \limits _z^\infty \frac{\partial }{\partial w} f(y,w) \frac{(w-z)^{-\alpha }}{\Gamma (1-\alpha )}dw. \end{aligned}$$

Now we use the above result and Fubini theorem to the following integral:

$$\begin{aligned} \int \limits _{\mathbf {R}} \int \limits _0^\infty&\left[ f(y,z+z') -f(y,z) \right] q_t(y,z)\frac{\alpha }{\Gamma (1-\alpha )}{z'}^{-1-\alpha }dz'dz \nonumber \\&= \int \limits _\mathbf {R}\int \limits _z^\infty \frac{\partial }{\partial w} f(y,w) \frac{(w-z)^{-\alpha }}{\Gamma (1-\alpha )}q_t(y,z)dwdz \nonumber \\&= \int \limits _{-\infty }^{\infty } \frac{\partial }{\partial w} f(y,w) \int \limits _{-\infty }^w \frac{(w-z)^{-\alpha }}{\Gamma (1-\alpha )}q_t(y,z)dzdw\nonumber \\&= -\int \limits _{-\infty }^{\infty } f(y,w) \frac{\partial }{\partial w}\int \limits _{-\infty }^w \frac{(w-z)^{-\alpha }}{\Gamma (1-\alpha )}q_t(y,z)dzdw. \end{aligned}$$
(13)

Since \(q_t(y,z)=0\) when \(z<0\), we get that

$$\begin{aligned} -\int \limits _{-\infty }^{\infty } f(y,w) \frac{\partial }{\partial w}\int \limits _{-\infty }^w \frac{(w-z)^{-\alpha }}{\Gamma (1-\alpha )}q_t(y,z)dzdw=-\int \limits _{\mathbf {R}} f(y,w) {_0}{D_w^{\alpha }}q_t(y,w)dw \end{aligned}$$

and hence

$$\begin{aligned}&\int \limits _{\mathbf {R}^2} \int \limits _0^\infty \left[ f(y,z+z')-f(y,z) \right] \frac{\alpha }{\Gamma (1-\alpha )}{z'}^{-1-\alpha } q_t(y,z)dz'dydz \nonumber \\&\quad =-\int \limits _{\mathbf {R}^2} f(y,z) {_0} {D_z^{\alpha }}q_t(y,z)dydz. \end{aligned}$$
(14)

Thus the adjoint operator has the form

$$\begin{aligned} A^+f(y,z)=- \frac{\partial }{\partial y}\left( F(y,z)f(y,z)\right) +\frac{\partial ^2}{\partial y^2} \left( \frac{1}{2}D(y,z)f(y,z)\right) -{_0}{D_z^{\alpha }}f(y,z). \end{aligned}$$
(15)

We take advantage of Eq. (9) to obtain

$$\begin{aligned} \frac{\partial }{\partial t}q_t(y,z)=- \frac{\partial }{\partial y}\left( F(y,z)q_t(y,z)\right) +\frac{\partial ^2}{\partial y^2} \left( \frac{1}{2}D(y,z)q_t(y,z)\right) -{_{0}}{D_z^{\alpha }}q_t(y,z). \end{aligned}$$
(16)

Now we denote by \(p_t(x)\) the PDF of \(X(t)\). We need to relate densities \(p_t(x)\) and \(q_t(y,z)\). In order to do this, for each \(\omega \in \Omega \) we denote the corresponding trajectories of the processes as \(X(t,\omega )\) and \((Y(t,\omega ), Z(t,\omega ))\). Let \(\delta _I(x)\) denote the usual indicator function for some fixed interval \(I\), that is the function having value \(1\) for \(x\in I\) and \(0\) otherwise. Then it is easy to observe that

$$\begin{aligned} \int \limits _Ip_t(x)dx=E \left[ \delta _I(X(t,\omega ))\right] . \end{aligned}$$
(17)

Following [6] we define an auxiliary function

$$\begin{aligned} H_t(t',\omega ,u)=\left\{ \begin{array}{ll} \delta _I(Y(t',\omega )) &{}\quad \mathrm{if } \; Z(t'-,\omega )\le t \le Z(t'-,\omega )+u \\ 0 &{}\quad \mathrm{otherwise} \end{array} \right. . \end{aligned}$$
(18)

For \(t',\omega \) such that \(Z(t'-,\omega )\le t\) this function takes value \(\delta _I(Y(t',\omega ))\) only when \(u\) is sufficiently large, that is \(u\) must exceed the gap between \(t\) and \(Z(t'-,\omega )\). Otherwise the function vanishes. Let \(\Delta Z(t,\omega )=Z(t,\omega )-Z(t-,\omega )\) denote the jumps of the process \(Z(t)\). Now we can write

$$\begin{aligned} \delta _I(X(t,\omega ))=\sum _{t'>0}H_t(t',\omega ,\Delta Z(t',\omega )), \end{aligned}$$
(19)

since all summands of the above are equal zero except from \(t'=S_\alpha (t,\omega )\), when we have \(Y(t',\omega )=Y(S_\alpha (t,\omega ),\omega )=X(t,\omega )\). In other words this right hand side of (19) expresses \(X(t,\omega )\) as a value of \(Y(t',\omega )\), where \(t'\) is the moment when the process \(Z(t',\omega )\) exceeds the level \(t\) for the first time. One can notice that jumps \(\Delta Z\) form a Poisson point process (Ppp) with values in \((0,\infty )\) and characteristic measure \(\nu \). This measure is defined as

$$\begin{aligned} \nu (A) = \frac{1}{t}E\left[ \sum _{s\le t}{\delta _A\left( \Delta Z(s)\right) }\right] , \quad t>0 \end{aligned}$$

and does not depend on \(t\) since Ppp is homogeneous in time. It is a known fact (see [11]) in theory of Lévy processes that the Lévy measure is the characteristic measure of the jump process. Since \(Z(t)=U_\alpha (t)\), \(\nu \) has the density \(\frac{\alpha }{\Gamma (1-\alpha )}u^{-1-\alpha }\). This enables us to use the compensation formula (see [14]) to write

$$\begin{aligned} E\left[ \sum _{t'>0}H_t(t',\omega ,\Delta Z(t',\omega ))\right] =E\left[ \int \limits _0^\infty \int \limits _{0}^\infty H_t(t',\omega ,u)\frac{\alpha u^{-1-\alpha }}{\Gamma (1-\alpha )} dudt' \right] . \end{aligned}$$
(20)

Left hand side of Eq. (20) yields

$$\begin{aligned}&E\left[ \sum _{t'>0}H_t(t',\omega ,\Delta Z(t',\omega ))\right] \nonumber \\&\quad =E\left[ \int \limits _{0}^\infty \delta _I(Y(t',\omega ))\delta _{[0,t]}\left( Z(t',\omega )\right) \frac{[t-Z(t',\omega )]^{-\alpha }}{\Gamma {(1-\alpha )}}dt' \right] \nonumber \\&\quad =\int \limits _I \int \limits _0^\infty \int \limits _0^t \frac{(t-z)^{-\alpha }}{\Gamma {(1-\alpha )}}q_{t'}(y,z)dzdt'dy \end{aligned}$$
(21)

Combining Eqs. (17), (19) and (21) gives

$$\begin{aligned} \int \limits _I p_t(x)dx&= \int \limits _I \int \limits _0^\infty \int \limits _0^t \frac{(t-z)^{-\alpha }}{\Gamma {(1-\alpha )}}q_{t'}(y,z)dzdt'dy \nonumber \\&= \int \limits _I \int \limits _0^\infty {_0} {I_t^{1-\alpha }}q_{t'}(y,t)dt'dy, \end{aligned}$$
(22)

where \({_0} {I_t^{1-\alpha }}\) denotes the Riemann–Liouville fractional integral of order \(1-\alpha \) [15]. Since the above equality holds for each interval \(I\), we have

$$\begin{aligned} p_t(x)=\int \limits _0^\infty {_0} {I_t^{1-\alpha }}q_{t'}(x,t)dt'. \end{aligned}$$
(23)

It also follows, using Fubini theorem, that

$$\begin{aligned} _{0}D_t^{1-\alpha }p_t(x)&= \frac{\partial }{\partial t} \int \limits _0^t \int \limits _0^\infty {_0}{ I_u^{1-\alpha }}q_{t'}(x,u)dt'\frac{(t-u)^{-1+\alpha }}{\Gamma {(\alpha )}}du \nonumber \\&= \int \limits _0^\infty \frac{\partial }{\partial t}\int \limits _0^t {_0}{I_u^{1-\alpha }}q_{t'}(x,u)\frac{(t-u)^{-1+\alpha }}{\Gamma {(\alpha )}}dudt' \nonumber \\&= \int \limits _0^\infty \frac{\partial }{\partial t} {_0} {I_t^\alpha } {_0}{I_t^{1-\alpha }}q_{t'}(x,t)dt'= \int \limits _0^\infty q_{t'}(x,t)dt'. \end{aligned}$$
(24)

Differentiation of Eq. (23) with respect to \(t\) yields

$$\begin{aligned} \frac{\partial }{\partial t}p_t(x)=\int \limits _0^\infty \frac{\partial }{\partial t} {_0}{I_t^{1-\alpha }} q_{t'}(x,t)dt'= \int \limits _0^\infty {_0} {D_t^{\alpha }} q_{t'}(x,t)dt'. \end{aligned}$$
(25)

Taking into account Eq. (16), the facts that \(\lim _{t'\rightarrow \infty }q_{t'}(x,t)=0\) and also \(q_{0}(x,t)=\delta _{(0,0)}(x,t)\), we obtain for \(t>0\)

$$\begin{aligned}&\frac{\partial }{\partial t}p_t(x) \nonumber \\&\quad =\int \limits _0^\infty \left( \frac{\partial ^2}{\partial x^2} \left( \frac{1}{2}D(x,t)q_{t'}(x,t)\right) - \frac{\partial }{\partial x}\left( F(x,t)q_{t'}(x,t)\right) -\frac{\partial }{\partial t'}q_{t'}(x,t)\right) dt' \nonumber \\&\quad =\frac{\partial ^2}{\partial x^2}\left( \frac{1}{2}D(x,t)\int \limits _0^\infty q_{t'}(x,t)dt'\right) - \frac{\partial }{\partial x}\left( F(x,t)\int \limits _0^\infty q_{t'}(x,t)dt'\right) . \end{aligned}$$
(26)

Finally, from Eq. (24) we get that

$$\begin{aligned} \frac{\partial }{\partial t}p_t(x)= \frac{\partial ^2}{\partial x^2}\left( \frac{1}{2}D(x,t){_0} {D_t^{1-\alpha }}p_t(x)\right) - \frac{\partial }{\partial x}\left( F(x,t) {_0}{D_t^{1-\alpha }}p_t(x)\right) , \end{aligned}$$
(27)

which ends the proof.

4 Conclusions

In this paper we have proved the stochastic representation of FFPE (1). In the proof we have applied the methods presented in [6] to the case of arbitrary drift \(F(x,t)\) and diffusion \(D(x,t)\) coefficient. Our result extends the one presented in [1], where the special factorized form of drift \(F(x,t)\) and diffusion \(D(x,t)\) was assumed. Moreover, we have pointed out some drawbacks in the derivation of the stochastic representation of FFPE studied in [1] and presented the correct proof in the general case.