\section{\label{sec1}Introduction}
Recently, the diffusion equations that generalize the usual one have
received considerable attention due to the broadness of their
physical applications, in particular, to the anomalous diffusion. In
fact, fractional diffusion equations and the nonlinear fractional
diffusion equations have been successfully applied to several
physical situations such as percolation of gases through porous
media [1], thin saturated regions in porous media [2], standard
solid-on-solid model for surface growth [3], thin liquid films
spreading under gravity [4], in the transport of fluid in porous
media and in viscous fingering [5], modeling of non-Markvian
dynamical processes in protein folding [6], relaxation to
equilibrium in system (such as polymer chains and membranes) with
long temporal memory [7], and anomalous transport in disordered
systems [8], diffusion on fractals [9], and the multi-physical
transport in porous media, such as electroosmosis[10-11]. The
properties concerning these equations have also been investigated.
For instance, in [12] boundary values problems for fractional
diffusion equations are studied, in [13] a fractional Fokker-Planck
equation is derived from a generalized master equation, in [14] the
behavior of fractional diffusion at the origin is analyzed and a
connection between the Fox H functions and the fractional diffusion
equations was investigated in [15]. Also a generalization of
Brownian motion to multidimensional anomalous diffusion is
considered by using fractional differential equation in [16], and a
fractional radial diffusion was considered in [17].

In this direction, we dedicate this work to investigate a fractional
diffusion equation which employs space and time fractional
derivatives by taking a time-dependent diffusion coefficient, an
external force and $\mathcal{N}$-dimensional into account. More
precisely, we focus our attention on the following equation
\begin{equation}
\frac{\pa^{\gamma}}{\pa
t^{\gamma}}\rho(x,t)=\int_{0}^{t}dt'\frac{1}{x^{\mathcal{N}-1}}\frac{\pa}{\pa
x}\{x^{\mathcal{N}-1}D(x,t-t')\frac{\pa^{\mu-1}}{\pa
x^{\mu-1}}[\rho(x,t)]^{\nu}\}-\frac{1}{x^{\mathcal{N}-1}}\frac{\pa}{\pa
x}\{x^{\mathcal{N}-1}F(x)\rho(x,t)\},
\end{equation}
where $0 <\gamma \leq1,\mu,\nu\in R$, $\mathcal{N}$ is the number of
dimensions, and the diffusion coefficient is given by
$D(x,t)=D(t)|x|^{-\theta}$, which is a spatial time-dependent
diffusion coefficient, and $F(x)=-\frac{\pa V(x)}{\pa x}$ is an
external force (drift) associated with the potential $V(x)$. Here,
we use the Caputo operator for the fractional derivative with
respect to time $t$ and the Riesz-Weyl operator for the fractional
derivative with respect to spatial $x$ [18]. We work with the
positive spatial variable $x$. Later on, we will extend the results
to the entire real $x$-axis by the use of symmetry (in other words,
we are working with ${\pa}/{\pa|x|}$ and ${\pa^{\mu
-1}}/{\pa|x|^{\mu -1}}$). Also, we employ, in general, the initial
condition $\rho(x,0)=\tilde{\rho}(x)$ ($\tilde{\rho}(x)$ is a given
function), and the boundary condition
$\rho(x\!\to\!\pm\infty,t)\rightarrow 0$. For Eq.(1), one can prove
that $\int_{-\infty}^{+\infty}dx x^{\mathcal{N}-1}\rho(x,t)$ is time
independent (hence, if $\rho (x,t)$ is normalized at $t=0$,it will
remain so forever). Indeed,if we write Eq.(1) as
$\pa_{t}^{\gamma}=-1/x^{\mathcal{N}-1}\pa_{x}(x^{\mathcal{N}-1}J)$,
and, for simplicity, assume the boundary conditions
$J(\pm\infty)=0$, it can be shown that $\int_{-\infty}^{+\infty}dx
x^{\mathcal{N}-1}\rho(x,t)$ is a constant of motion (see [19] and
references therein). Note that Eq.(1) recovers the usual radial
diffusion equation with memory effect for $\gamma=1$. When
$(\mu,\gamma,\nu)=(2,1,1)$ and $\mathcal{N}=1$, Eq.(1) recovers the
standard Fokker-Planck equation in the presence of a drift taking
memory effects and $\mathcal{N}$-dimensional into account. The
particular case $F(x)=0$ (no drift) and $D(x,t)=D\delta(t)$ with
$(\mu,\gamma)=(2,1)$ has been considered by spohn [3]. The case
$D(x,t)=D\delta(t)|x|^{-\theta}$ with $(\mu,\nu)=(2,1)$ and the case
$D(x,t)=D(t)$ with $(\mu,\nu)=(2,1)$ have been investigated in [20]
and [21], respectively. In [22], the solution for Eq.(1) with the
boundary condition $\rho(0,t)=\rho(L,t)=0$ was investigated.

Explicit solutions play an important role in analyzing physical
situations, since they contain, in principle, precise information
about the system. In particular, they can be used as a useful guide
to control the accuracy of numerical solutions. For these reasons,
we dedicated to this work to investigate the solutions to Eq.(1). We
consider different scenarios involving the diffusion coefficient and
external force. Firstly, in the case of the absence of external
force, we consider a spatial time-dependent diffusion coefficient
given by $D(x,t)=D|x|^{-\theta}t^{\alpha-1}/\Gamma(\alpha)$ and
$D(x,t)=D|x|^{-\theta}\delta(t)$. Secondly, the presence of external
force is investigated, which is given by $F(x)\propto
x|x|^{\alpha-1}$. Thirdly, a mixing between the spatial and time
fractional derivatives case is investigated. Lastly, we consider a
particular case as $\gamma=1$, $D(t)=D\delta(t)$, $F(x)=-\mathcal
{K}x$, and $\theta$, $\mu$, $\nu$ arbitraries. In all the above
situations, Eq.(1) satisfies the initial condition
$\rho(x,0)=\tilde{\rho}(x)$ ($\tilde{\rho}(x)$ is a given function),
and the boundary condition $\rho(\pm \infty,t)=0$.The remainder of
this paper goes as follow. In Sec.2, we obtain the exact solutions
for the previous cases. In Sec.3, we present our conclusions.

\section{\label{sec2} Exact solutions for different case}

In this section, we start our discussion by considering Eq.(1) in
the absence of external force with
$D(t)=Dt^{\alpha-1}/\Gamma(\alpha)$ (
$D(t)=D\delta(t)$),$(\mu,\nu)=(2,1)$ and $\gamma$, $\theta$
arbitrary. For this case, Eq.(1) reads
\begin{equation}
\frac{\pa^{\gamma}}{\pa
t^{\gamma}}\rho(x,t)=\int_{0}^{t}dt'D(t-t')\frac{1}{x^{\mathcal{N}-1}}\frac{\pa}{\pa
x}\{x^{\mathcal{N}-1-\theta}\frac{\pa}{\pa x}{\rho(x,t)}\}.
\end{equation}
Here, we use the Caputo operator [18] for the fractional derivative
with respect to time $t$. By employing the Laplace transform in
Eq.(2), we obtain
\begin{equation}
\frac{1}{x^{\mathcal{N}-1}}\tilde{D}(s)\frac{\pa}{\pa
x}\{x^{\mathcal{N}-1-\theta}\frac{\pa}{\pa
x}\tilde{\rho}(x,s)\}-s^{\gamma}\tilde{\rho}(x,s)=-s^{\gamma-1}\rho(x,0),
\end{equation}
where $\tilde{\rho}(x,s)=\mathscr{L} \{\rho(x,t)\}$,
$\tilde{D}(s)=\mathscr{L}\{D(t)\}$, and $
\mathscr{L}\{f(t)\}=\int_0^\infty dte^{-st}f(t)$ denotes the Laplace
transform of the function $f$. This equation can be solved by Green
function method [23]. By substituting
\begin{equation}
\tilde{\rho}(x,s)=\int dx'x'^{\mathcal{N}-1}\tilde{\mathcal
{G}}(x-x',s)\tilde{\rho}(x')
\end{equation}
into Eq.(3), which yields
\begin{equation}
\frac{1}{x^{\mathcal{N}-1}}\tilde{D}(s)\frac{\pa}{\pa
x}\{x^{\mathcal{N}-1-\theta}\frac{\pa}{\pa x}\tilde{\mathcal
{G}}(x,s)\}-s^{\gamma}\tilde{\mathcal
{G}}(x,s)=-s^{\gamma-1}\mathcal {G}(x,0).
\end{equation}
In order to solve Eq.(5), it is convenient to perform the transform
[24]
\begin{equation}\label{6}
y=A(s)x^{v},\quad \mathcal {G}(x,s)=y^{\delta}Z(y)
\end{equation}
to translate Eq.(5) into the second-order Bessel equation as
\begin{equation}\label{7}
y^2\frac{\pa^2 Z}{\pa y^2}+y\frac{\pa Z}{\pa
y}-(\lambda^2+y^2)Z(y)=-\frac{y^{\delta+1/v}}{sA(s)^{(\mathcal{N}-1)/v}}\delta((\frac{y}{A(s)})^{\frac{1}{v}})
\end{equation}
with parameter $\lambda^2$ under the following conditions:
\begin{equation}
v=\frac{2+\theta}{2},A(s)=\frac{1}{v}[\frac{s^\gamma}{\tilde{D}(s)}]^\frac{1}{2},\lambda=\frac{2+\theta-\mathcal{N}}{2+\theta},\delta=\frac{2+\theta-\mathcal{N}}{2+\theta}.
\end{equation}
Since Eq.(5) should fit the boundary condition $\mathcal
{G}(\pm\infty,t)=0$, i.e. $\tilde{\mathcal {G}}(\pm\infty,s)=0$, we
get the solution of Eq.(5)
\begin{equation}
\tilde{\mathcal {G}}(x,s)=C(s)y^{\delta}K_{\lambda}(y),
\end{equation}
where $K_{n}(x)$ is the modified Bessel function of second kind, and
$C(s)$ can be determined by the normalization of $\mathcal
{G}(x,t)$, i.e. $\int_0^\infty
dxx^{\mathcal{N}-1}\tilde{\mathcal{G}}(x,s)=\frac{1}{2s}$. After
some calculations, we obtain
\begin{equation}
\tilde{\mathcal
{G}}(x,s)=\frac{2+\theta}{\Gamma(\frac{\mathcal{N}}{2+\theta})s}
(\frac{1}{2+\theta}(\frac{s^{\gamma}}{\tilde{D}(s)})^\frac{1}{2})^{\frac{2+\theta+\mathcal{N}}{2+\theta}}|x|^{\frac{2+\theta-\mathcal{N}}{2}}
K_{\frac{2+\theta-\mathcal{N}}{2+\theta}}(\frac{2}{2+\theta}(\frac{s^{\gamma}}{\tilde{D}(s)})^\frac{1}{2}|x|^{\frac{2+\theta}{2}}).
\end{equation}
here,we have used the formula
\begin{equation}
\int_0^\infty dy\cdot y^v
K_\lambda(ay)=2^{v-1}a^{-v-1}\Gamma(\frac{1+v+\lambda}{2})\Gamma(\frac{1+v-\lambda}{2}).
\end{equation}
\textbf{  Case 1. } $D(t)=D\delta(t)$, i.e. $\tilde{D}(s)=D$.

Since $K_\lambda (x)=\frac{1}{2}H_{0,2}^{2,0}\left[\frac{x^2}{4}\left|_{(-\lambda /2,1)(\lambda /2,1)}
\right.\right]$, we can get the Laplace inverse of $\tilde{\mathcal {G}}(x,s)$ by applying the property of the
Laplace inverse of Fox function [25], which yields
\begin{equation}
\mathcal {G}(x,t)=\frac{2+\theta}{2\Gamma(\frac{\mathcal{N}}{2+\theta})}(\frac{1}{(2+\theta)^2
Dt^{\gamma}})^{\frac{\mathcal{N}}{2+\theta}}H_{1,2}^{2,0}\left[\frac{|x|^{2+\theta}}{(2+\theta)^2
Dt^{\gamma}}\left|\begin{array}[c]{cr}(1-\frac{\mathcal{N}\gamma}{2+\theta},\gamma)\\(0,1),&(\frac{2+\theta-\mathcal{N}}{2+\theta},1)\end{array}\right.\right],
\end{equation}
\begin{center}
\begin{figure}[thb]
\scalebox{1.3}[1.1]{\includegraphics{f1.eps}}
 \caption{\label{Fig.1}
\small{The behavior of $\rho(x,t)$ in Eq.(13) is illustrated by considering
$C(\gamma,\theta,\mathcal{N})\mathcal{G}(x,t)$ versus $\frac{x^{2+\theta}}{(2+\theta)^2 Dt^{\gamma}}$ for
typical values of $\gamma$ and $\theta$, where
$C(\gamma,\theta,\mathcal{N})=\frac{2\Gamma(\mathcal{N}/(2+\theta))}{2+\theta}((2+\theta)^2
Dt^{\gamma})^{\mathcal{N}/(2+\theta)}$. Here, for simplicity, we consider $\rho(x,0)=x^{1-\mathcal{N}}\delta(x)$
 }}
 \end{figure}
\end{center}
where $H_{p,q}^{m,n}\left[x\left|_{(b_1,B_1),...,(b_q,B_q)}^{(a_1,A_1),...,(a_p,A_p)}\right.\right]$ is the FOX
function. Thus, we can find the solution by substituting Eq.(12) into Eq.(4), which yields (see fig.1)
\begin{equation}
\rho(x,t)=\frac{2+\theta}{2\Gamma(\frac{\mathcal{N}}{2+\theta})}(\frac{1}{(2+\theta)^2
Dt^{\gamma}})^{\frac{\mathcal{N}}{2+\theta}}\int_{-\infty}^{+\infty}dx'x'^{\mathcal{N}-1}\tilde{\rho}(x')H_{1
,2}^{2,0}\left[\frac{|x-x'|^{2+\theta}}{(2+\theta)^2
Dt^{\gamma}}\left|\begin{array}[c]{cr}(1-\frac{\mathcal{N}\gamma}{2+\theta},\gamma)\\(0,1),&(\frac{2+\theta-\mathcal{N}}{2+\theta},1)\end{array}\right.\right].
\end{equation}
In Fig.1, we show the behavior of the above equation by considering
typical values of $\gamma$, $\mathcal{N}$ and $\theta$ with
$\tilde{\rho}(x)=x^{1-\mathcal{N}}\delta(x)$. Note that the solution
obtained above leads to an anomalous spreading of the initial
condition due to the presence of the time fractional derivative and
the spatial time-dependent diffusion coefficient. The present of Fox
function in Eq.(5) may be associated to the changes produced on the
probability density function for a waiting time that, for this case,
has a long tailed behavior, in contrast to the usual one. At this
point, it is interesting to analyze the asymptotic behavior of
Eq.(13) [25]. For simplicity, we consider
$\tilde{\rho}(x)=x^{1-\mathcal{N}}\delta(x)$, so
$\rho(x,t)=\mathcal{G}(x,t)$ and the asymptotic behavior of
$\rho(x,t)$ is
$$\rho(x,t)\sim
\frac{2+\theta}{2\Gamma(\mathcal{N}/(2+\theta))}(2-\gamma)^{-\frac{1}{2}}\gamma^{\frac{\mathcal{N}\gamma}{(2+\theta)(2-\gamma)}-\frac{1}{2}}(\frac{1}{(2+\theta)^2
Dt^{\gamma}})^{\frac{\mathcal{N}}{(2+\theta)(2-\gamma)}}|x|^{\frac{\mathcal{N}(\gamma-1)}{2-\gamma}}$$
\begin{equation}
\times
exp(-(2-\gamma)\gamma^{\frac{\gamma}{2-\gamma}}(\frac{|x|^{2+\theta}}{(2+\theta)^2
Dt^{\gamma}})^{\frac{1}{2-\gamma}}).
\end{equation}
In this direction, Eq.(14) can be considered as an extension of the
asymptotic behavior of homogeneous and isotropic random walk
models [26].\\
\textbf{  Case 2. } $D(t)=\frac{Dt^{\alpha-1}}{\Gamma(\alpha)}$,
i.e. $\tilde{D}(s)=Ds^{-\alpha}$.

Using the same method as in case 1, we obtain
\begin{equation}
\mathcal {G}(x,t)=\frac{2+\theta}{2\Gamma(\frac{\mathcal{N}}{2+\theta})}(\frac{1}{(2+\theta)^2
Dt^{\gamma+\alpha}})^{\frac{\mathcal{N}}{2+\theta}}H_{1,2}^{2,0}\left[\frac{|x|^{2+\theta}}{(2+\theta)^2
Dt^{\gamma+\alpha}}\left|\begin{array}[c]{cr}(1-\frac{\mathcal{N}(\gamma+\alpha)}{2+\theta},\gamma+\alpha)\\(0,1),&(\frac{2+\theta-\mathcal{N}}{2+\theta},1)\end{array}\right.\right],
\end{equation}
and
$$
\rho(x,t)=\frac{2+\theta}{2\Gamma(\frac{\mathcal{N}}{2+\theta})}(\frac{1}{(2+\theta)^2
Dt^{\gamma+\alpha}})^{\frac{\mathcal{N}}{2+\theta}}\int_{-\infty}^{+\infty}dx'x'^{\mathcal{N}-1}\tilde{\rho}(x')$$
\begin{equation}
\quad \quad \quad \quad \quad \quad \quad \quad \times H_{1,2}^{2,0}\left[\frac{|x-x'|^{2+\theta}}{(2+\theta)^2
Dt^{\gamma+\alpha}}\left|\begin{array}[c]{cr}(1-\frac{\mathcal{N}(\gamma+\alpha)}{2+\theta},\gamma+\alpha)\\(0,1),&(\frac{2+\theta-\mathcal{N}}{2+\theta},1)\end{array}\right.\right].
\end{equation}

Let us go back to Eq.(1), and consider the external force
$F(x)\varpropto x|x|^{\alpha-1}$,$D(x,t)=D\delta(t)|x|^{-\theta}$
and $\mu=2$, $\nu=1$. In this case ,analytical solution can not
easily be obtained for a generic $\alpha$, $\theta$. However, for
$\theta\neq 0$, and $\alpha+\theta+1=0$. By following the same
procedure as the one in the above case, an exact solution can be
obtained and it is given by
\begin{equation}
\rho(x,t)=\frac{2+\theta}{2\Gamma(\frac{\mathcal{N}+\mathcal{K}/D}{2+\theta})\Gamma(\frac{3+\theta-\mathcal{N}}{2+\theta})}(\frac{1}{(2+\theta)^2
Dt^{\gamma}})^{\frac{\mathcal{N}}{2+\theta}}H_{1,2}^{2,0}\left[\frac{|x|^{2+\theta}}{(2+\theta)^2
Dt^{\gamma}}\left|\begin{array}[c]{cr}(1-\frac{\mathcal{N}\gamma}{2+\theta},\gamma)\\(\frac{\mathcal{K}}{(2+\theta)D},1),(\frac{2+\theta-\mathcal{N}}{2+\theta},1)\end{array}\right.\right],
\end{equation}
where, for simplicity, we are considering the initial condition
$\rho(x,0)=\delta(x)$, and the external force (drift)
$F(x)=\mathcal{K}x|x|^{\alpha-1}$. The second moment is given by
$<x^2>\varpropto t^{\frac{2\gamma}{2+\theta}}$, which scales with
the exponent $\frac{2\gamma}{2+\theta}$ and clearly depends only on
$\gamma$ and $\theta$. So, when $\frac{2\gamma}{2+\theta}<1$, $=1$
and $>1$, the system is sub-diffusion, normal diffusion and
supper-diffusion respectively.

Let us now discuss Eq.(1) in the presence of the external force by
employing a mixing between the spatial and time fractional
derivatives. For simplicity, we consider
$(\mathcal{N},\nu,\theta)=(1,1,0)$ and $F(x)=-\mathcal{K}x$. For
this case, Eq.(1) reads
\begin{equation}
\frac{\pa^{\gamma}}{\pa t^{\gamma}}\rho(x,t)=\int_0^t
dt'D(t-t')\frac{\pa^{\mu}}{\pa x^{\mu}}\rho(x,t')-\frac{\pa}{\pa
x}[F(x)\rho(x,t)].
\end{equation}
Applying the Fourier and Laplace transform to Eq.(1) and employing
the Riez representation for the spatial fractional derivatives, we
have
\begin{equation}
s^{\gamma}\hat{\tilde{\rho}}(k,s)-s^{r-1}\hat{\rho}(k,0)=-\tilde{D}(s)|k|^{\mu}\hat{\tilde{\rho}}(k,s)-\mathcal{K}k\frac{d}{dk}\hat{\tilde{\rho}}(k,s),
\end{equation}
where $\hat{\rho}(k,t)=\mathcal{F}\{\rho(x,t)\}=\int_{-\infty}
^{+\infty} \rho(x,t)e^{-ikx}dx$ , so $\hat{\rho}(k,0)=1$. Here, we
consider the diffusion coefficient given by
$D(t)=\frac{Dt^{\alpha-1}}{\Gamma(\alpha)}$, i.e.
$\tilde{D}(s)=Ds^{-\alpha}$ and $\rho(x,0)=\delta(x)$. So we have
\begin{equation}
\hat{\tilde{\rho}}(k,s)=\sum_{n=0}^{\infty}\frac{1}{n!}(\frac{\tilde{D}(s)|k|^{\mu}}{\mathcal{K}\mu})^n
e^{-\frac{\tilde{D}(s)|k|^{\mu}}{\mathcal{K}\mu}}\frac{s^{\gamma-1}}{s^\gamma
+n\mu \mathcal{K}}.
\end{equation}
In order to perform the inverse of Laplace transform on Eq.(20), we express $e^{-z}$ in terms of Fox function,
i.e. $e^{-z}=H_{0,1}^{1,0}\left[z\left|_{(0,1)}\right.\right]$.  Then we can get the Laplace inverse of
$\tilde{\rho}(x,s)$ by applying the property of the Laplace inverse of Fox function [25]
\begin{equation}
\hat{\rho}(k,t)=\sum_{n=0}^{\infty}\frac{1}{n!}\int_0^t dt' \frac{1}{t'}H_{1,1}^{1,0}\left[\frac{D|k|^\mu
t'^\alpha}{\mathcal{K}\mu}\left|^{(0,-\alpha)}_{(n,1)}\right.\right]
E_{\alpha,1}(-n\mu\mathcal{K}(t-t')^{\gamma}),
\end{equation}
where, $E_{\alpha,\beta}(x)$ is the Mittage-Leffler function defined as $E_{\alpha,\beta}(x)=\sum_{n=0}^\infty
\frac{x^n}{\Gamma(n\alpha+\beta)}$ [18]. Here, we used the property of the Laplace transform of convolution
formula, i.e. $\mathcal{L}[f\ast g]=f\cdot g$, where $f\ast g=\int_0^t dt'f(t-t')g(t')$. Note that the solution
of Eq.(18) is a stationary one given in terms of L$\acute{e}$vy distributions. This feature is a characteristic
of the presence of the spatial derivatives in the diffusion equation which changes the probability for a jump
length (see [27] and references therein). In order to get the solution for Eq.(18), we need to get the Fourier
inverse of $\hat{\rho}(k,t)$. Therefore, we only need to perform the inverse of Fourier transform on $H_{1,
1}^{1,0}\left[\frac{D|k|^\mu t'^\alpha}{\mathcal{K}\mu}\left|_{(n,1)}^{(0,-\alpha)}\right.\right]$. By employing
the procedure presented in [28], we have
\begin{equation}
\mathcal{F}^{-1}\{H_{1 \quad 1}^{1 \quad 0}[\frac{D|k|^\mu
t'^\alpha}{\mathcal{K}\mu}|_{(n,1)}^{0,-\alpha}]\}=\frac{1}{\mu\sqrt{\pi}|x|}H_{2,2}^{1,
1}\left[\frac{|x|}{2}(\frac{\mathcal{K}\mu}{Dt'^{\alpha}})^{1/\mu}\left|\begin{array}[c]{cr}(1-n,\frac{1}{\mu}),&(0,\frac{\alpha}{\mu})\\(\frac{1}{2},\frac{1}{2}),&(1,\frac{1}{2})
\end{array}\right.\right],
\end{equation}
and
$$
\rho(x,t)=\sum_{n=0}^{\infty}\int_0^t dt'\frac{1}{n!\mu\sqrt{\pi}|x|t'}H_{2,2}^{1,
1}\left[\frac{|x|}{2}(\frac{\mathcal{K}\mu}{Dt'^{\alpha}})^{1/\mu}\left|\begin{array}[c]{cr}(1-n,\frac{1}{\mu}),&(0,\frac{\alpha}{\mu})\\(\frac{1}{2},\frac{1}{2}),&(1,\frac{1}{2})
\end{array}\right.\right]
$$
\begin{equation}
\times E_{\alpha,1}(-n\mu\mathcal{K}(t-t')^{\gamma}).
\end{equation}



Now, we consider a particular case of Eq.(1) for $\gamma=1$ and
nonzero values of $\mu$ and $\theta$, and consider a linear
drift,i.e. $F(x)=-\mathcal {K}x$. For simplicity, we employ
$D(t)=D\delta(t)$ and the initial condition $\rho(x,0)=\delta(x)$,
then Eq.(1) yields to
\begin{equation}
\frac{\pa}{\pa t}\rho(x,t)=\frac{D}{x^{\mathcal{N}-1}}\frac{\pa}{\pa
x}\{x^{\mathcal{N}-1}|x|^{-\theta}\frac{\pa^{\mu-1}}{\pa
x^{\mu-1}}[\rho(x,t)]^{\nu}\}-\frac{1}{x^{\mathcal{N}-1}}\frac{\pa}{\pa
x}\{x^{\mathcal{N}-1}F(x)\rho(x,t)\}.
\end{equation}
Let us investigated time dependent solutions for Eq.(24). We use
similarity methods to reduce Eq.(24) to ordinary differential
equations. The explicit form for these ordinary differential
equations depends on the boundary conditions or restrictions in the
form of conservation laws. In this direction, we restrict our
analysis to find solution that can be expressed as a scaled function
of the type [29]
\begin{equation}
\rho(x,t)=\frac{1}{\phi(t)^{\mathcal{N}}}\bar{\rho}(z),\quad
z=\frac{|x|}{\phi(t)}.
\end{equation}
Inserting Eq.(25) into Eq.(24), we obtain
\begin{equation}
-(\frac{\dot{\phi(t)}}{\phi(t)^2}+\frac{\mathcal{K}}{\phi(t)})\frac{\pa}{\pa
z}[z^{\mathcal{N}}\bar{\rho}(z)]=\frac{D}{\phi(t)^{\mathcal{N}(\nu-1)+\theta+\mu+1}}\frac{\pa}{\pa
z}[z^{\mathcal{N}-1-\theta}\frac{\pa^{\mu-1}}{\pa
z^{\mu-1}}\bar{\rho}(z)^{\nu}].
\end{equation}
By choosing the ansatz
\begin{equation}
\frac{\dot{\phi(t)}}{\phi(t)^2}+\frac{\mathcal{K}}{\phi(t)}=\frac{kD}{\phi(t)^{\mathcal{N}(\nu-1)+\theta+\mu+1}},
\end{equation}
where $k$ is an arbitrary constant which can be determined by the
normalization condition. By solving Eq.(27), we have that
\begin{equation}
\phi(t)=[(\phi(0))^{\xi}e^{-\xi\mathcal{K}t}+\frac{Dk}{\mathcal{K}}(1-e^{\xi\mathcal{K}t})]^{\frac{1}{\xi}},
\end{equation}
where $\xi=\mathcal{N}(\nu-1)+\theta+\mu$. By substituting Eq.(27)
into Eq.(26), we obtain
\begin{equation}
\frac{\pa}{\pa z}[z^{\mathcal{N}-1-\theta}\frac{\pa^{\mu-1}}{\pa
z^{\mu-1}}\bar{\rho}(z)^{\nu}]=-k\frac{\pa}{\pa
z}[z^{\mathcal{N}}\bar{\rho}(z)].
\end{equation}
Then, we perform an integration and the result is
\begin{equation}
z^{\mathcal{N}-1-\theta}\frac{\pa^{\mu-1}}{\pa
z^{\mu-1}}\bar{\rho}(z)^{\nu}=-kz^{\mathcal{N}}\bar{\rho}(z)+\mathcal{C},
\end{equation}
where $\mathcal{C}$ is another arbitrary constant. Also, we use the
following generic result [30]:
\begin{equation}
D_{x}^{\delta}[x^{\alpha}(a+bx)^{\beta}]=a^{\delta}\frac{\Gamma[\alpha+1]}{\Gamma[\alpha+1-\delta]}x^{\alpha-\delta}(a+bx)^{\beta-\delta}
\end{equation}
with $D_{x}^{\delta}\equiv d^{\delta}/dx^{\delta}$ and
$\delta=\alpha+\beta+1$. By defining $g(x)\equiv
x^{\frac{\alpha}{\nu}}(a+bx)^{\frac{\beta}{\nu}}$ and $\lambda\equiv
\alpha(1-\frac{1}{\nu})-\delta$, and rearranging the indices,
Eq.(31) can be rewritten as follows:
\begin{equation}
D_{x}^{\delta}[g(x)]^{\nu}=a^{\delta}\frac{\Gamma[\alpha+1]}{\Gamma[\alpha+1-\delta]}x^{\lambda}g(x).
\end{equation}
For this case, we consider the ansatz
$\bar{\rho}(z)=\mathcal{A}z^{\frac{\alpha}{\nu}}(1+bz)^{\frac{\beta}{\nu}}$.
By using the property of Eq.(32) in Eq.(30) and ,for simplicity,
choosing $\mathcal{C}=0$, we find
$$\alpha=\frac{(2-\mu)(\mu+\theta)}{1-2\mu-\theta},$$
\begin{equation}
\beta=-\frac{(\mu-1)(\mu-2)}{1-2\mu-\theta},
\end{equation}
$$\nu=\frac{2-\mu}{1+\mu+\theta}.$$
In this case, we have
\begin{equation}
\rho(x,t)=\frac{\mathcal{A}}{\phi(t)^{\mathcal{N}}}[\frac{z^{(\mu+\theta)(1+\mu+\theta)}}{(1+bz)^{(1-\mu)(1+\mu+\theta)}}]^{\frac{1}{1-2\mu-\theta}},
\end{equation}
where $\phi(t)$ is given
above,$\mathcal{A}=[-k\frac{\Gamma(-\beta)}{\Gamma(\alpha+1)}]^{\frac{\mu+\theta+1}{1-2\mu-\theta}}$
and $b$ is an arbitrary constant (to be taken, later on, as $\pm1$
according to the specific solutions that are studied). Several
regions can be analyzed. For simplicity, we illustrate two of them:
$-\infty<\mu<-1-\theta$ with $\theta\geq0$, and $0<\mu<1/2$ with
$0\leq\theta<1/2-\mu$. Let us start by considering the region
$-\infty<\mu<-1-\theta$.  Without loss of generality, we choose
$b=-1$. The normalization condition implies(see Fig. 2)
\begin{equation}
\mathcal{A}\int_{-1}^{1}[\frac{z^{(\mu+\theta)(1+\mu+\theta)}}{(1+bz)^{(1-\mu)(1+\mu+\theta)}}]^{\frac{1}{1-2\mu-\theta}}dz=1.
\end{equation}
So
\begin{equation}
\mathcal{A}=\frac{\Gamma[1-\mu-\theta]}{2\Gamma[\frac{\mu^2+\mu\theta-2\theta-2\mu}{1-2\mu-\theta}]\Gamma[\frac{1-\mu+\mu^2+\theta^2+2\mu\theta}{1-2\mu-\theta}]}.
\end{equation}
\begin{center}
\begin{figure}[thb]
\scalebox{1.3}[1.1]{\includegraphics{f2.eps}}
 \caption{\label{FigII}
\small{Behavior of $\phi(t)^{\mathcal{N}}\rho(x,t)$ versus
$x/\phi(t)$, which illustrates Eq.(32) with typical values for $\mu$
and $\theta$ satisfying $0<\mu<-1-\theta$ and $\theta\geq0$. we
notice that the distribution vanishes at the abcissa equal $\pm1$,
and remains zero outside of this interval.
 }}
 \end{figure}
\end{center}

Let us now analyze the region $0<\mu<1/2$ with
$0\leq\theta<1/2-\mu$. Again without the loss of generality, we
choose $b=1$. The normalization condition implies(see Fig. 3)
\begin{equation}
\mathcal{A}=\frac{\Gamma[\frac{1+\theta-\mu^2-\mu\theta}{1-2\mu-\theta}]}{2\Gamma[\frac{1-\mu+\mu^2+\theta^2+2\mu\theta}{1-2\mu-\theta}]\Gamma[\mu+\theta]}.
\end{equation}
\begin{center}
\begin{figure}[thb]
\scalebox{1.3}[1.1]{\includegraphics{f3.eps}}
 \caption{\label{FigII}
\small{Behavior of $\phi(t)^{\mathcal{N}}\rho(x,t)$ versus
$x/\phi(t)$, which illustrates Eq.(32) with typical values for $\mu$
and $\theta$ satisfying $0<\mu<1/2$ and $0\leq\theta<1/2-\mu$. we
notice that the distribution vanishes at the infinity.
 }}
 \end{figure}
\end{center}

Let us finally mention a connection between the results obtained
here and the solutions that arise from the optimization of the
nonextensive entropy [31]. These distributions do not coincide for
arbitrary value of x. However, the comparison of the
$|x|\rightarrow\infty$ asymptotic behaviors enables us to identify
the type of Tails. By identifying the behavior exhibited in Eq.(34)
with the asymptotic behaviors $1/|x|^{2/(q-1)}$ that appears in [31]
for the entropic problem, we obtain
\begin{equation}
q=\frac{3+\mu+\theta}{1+\mu+\theta}.
\end{equation}
This relation recovers the situation for $\theta=0$.

\section{\label{sec3} Summary and Conclusions}

We have analyzed a generalized fractional diffusion equation which
presents the spatial and time fractional derivatives, includes a
linear external force $F(x)=-\mathcal{K}x$ and a spatial
time-dependent diffusion coefficient $D(x,t)=D(t)|x|^{-\theta}$,
also takes $\mathcal{N}$-dimensional into account. By using Laplace
and Fourier transform, the Green function method and normalized
scaled function, we can find the explicit solutions $\rho(x,t)$
which subjects to the natural boundary condition
$\rho(\pm\infty,t)=0$ and the initial condition
$\rho(x,t)=\tilde{\rho}(x)$. In a word, we have extended the results
previously obtained by the other authors by inviting an linear
external force, a spatial time-dependent diffusion coefficient and
$\mathcal{N}$-dimensional case. We have also discussed the
connection with nonextensive statistics, providing the relation
between our solutions and those obtained within the maximum entropy
principle by using the Tsallis entropy. Finally, we expect that the
results obtained here may be useful to the discission of the
anomalous diffusion systems where fractional diffusion equations
play an important role.


{\bf Acknowledgments:} We would like to thank Hangzhou Dianzi
University for partial financial support and be grateful to the
anonymous referees for useful comments and suggestions.

\thebibliography{00}
 \baselineskip 13pt

\bibitem{mur90} M.Muskat, The Flow of Homogeneous Fluid Through Porous Media, McGraw-Hill, New York, 1937.

\bibitem{pol96} P.Y. Polubarinova-Kochina, Theory of Ground Water Movement, Princeton University Press,Princeton, 1962.

\bibitem{spo93} H. Spohn, J. Phys. 13 (1993) 69.

\bibitem{buc00} J.Buckmaster, J.Fluid Mech. 81 ( 1983) 735.

\bibitem{gro06} P. Grosfils and J.P. Boon, Physica A 362 (2006) 168.

\bibitem{plo04} S.S. Plotkin and P.G. Wolynes, Phys. Rev. Lett. 80 (1998) 5015 .

\bibitem{cro04} D.S.F. Crothers, D. Holland, Y.P. Kalmykov and W.T. Coffey, J. Mol. Liq. 114 (2004) 27.

\bibitem{mez99} R. Metzler, E. Barkai and J. Klafter, Physica A 266 (1999) 343.

\bibitem{cam04} D. Campos, V. Mendez and J. Fort, Phys. Rev. E 69
(2004) 031115.

\bibitem{wan07} M. Wang and S.Y. Chen, J. Colloids Interface Sci. 314
(2007) 264.

\bibitem{wan06} M. Wang, JK. Wang, S.Y. Chen and N. Pan, J. Colloids Interface
Sci. 304 (2006) 246.

\bibitem{met00} R.Metzler and J. Klafter, Phys. Rep. 339 (2000) 1.

\bibitem{met02} R.Metzler and T.F.Nonnenmacher, Chem. Phys. 284 (2002) 67.

\bibitem{wes02} B.J.West, M.Bologna and P.Grigolini, Physics of
Fractal Operators, Springer, New York, 2002.

\bibitem{mai05} F.Mainardi, G.Pagnini and R.K.Saxena,
J.Comput.Appl.Math. 178 (2005) 321.

\bibitem{uch02} V.V.Uchaikin, Chem. Phys. 284 (2002) 507.

\bibitem{ach04} B.N.N.Achar and J.W.Hanneken, J. Mol.
Liq. 114 (2004) 147.

\bibitem{pod99} I. Podlubny, Fractional differential equations, Academic
Press, San Diego, CA, 1999. 54

\bibitem{lez03} E.K. Lenzi, L.C. Malacarne, R.S. Mendes and I.T.
Pedron, Physica A 319 (2003) 245.

\bibitem{len03} E.K. Lenzi, R.S. Mendes, Kwok Sau Fa and
L.C.Malacame, J. Math. Phys. 44 (2003) 2179.

\bibitem{len04} E.K. Lenzi, R.S. Mendes and Kwok Sau Fa, J. Math. Phys.
45 (2004) 3444.

\bibitem{len06} E.K. Lenzi, R.S. Mendes, G. Goncalves, M.K. Lenzi
and L.R. da Silva, Physica A 360 (2006) 215.

\bibitem{mor53} M.P.Morse and H.Feshbach, Methods of Theoretical Physics,
McGraw-Hill, New York, 1953.

\bibitem{ren06} F.Y. Ren, J.R. Liang, W.Y.Qiu and J.B.Xiao,
J.Phys.A:Math.Gen, 39 (2006) 4911.

\bibitem{mat78} A.M. Mathai and R.K.Saxtena, The H-function with
Application in Statistics and Other Disciplines, Wiley Eastern, New
Delhi, 1978.

\bibitem{kla91} J.Klafter, G.Zumoften and A.Blumen, J.phys.A 25 (1991)
4835.

\bibitem{sch07} A.Schot, M.K.Lenzi, L.R.Evangelista, etal, Phys.
Lett. A 366 (2007) 346.

\bibitem{lan06} T.A.M.Langlands, Physica A 367 (2006) 136.


\bibitem{tsa02} C.Tsallis and E.K.Lenzi, Chemical physics, 284 (2002)
341

\bibitem{bol00} M.Bologna, C.Tsallis and P.Grigolini, Phys.Rev.E
62 (2000) 2213.

\bibitem{tsa96} C.Tsallis and D.J.Bukman, Phys.Rev.E
54 (1996) R2197.

\end{document}
