the heat flow can be disrupted by the response of the material. It has been proved (see e.g., [@V; @LS]) that in a material with thermal memory $$\label{2014.5}
e(t,x)=\bar{\beta} u(t,x)+\int^t_0 n(t-s)u(s,x)ds$$ holds with some appropriate constant $\bar{\beta}$ and kernel $n$. Typically, $n(t)$ is a positive decreasing function which blows up near $t=0$, for instance $n(t)=t^{-\alpha}$ for $\alpha\in (0,1)$. In this case the convolution implies that nearer past affects the present more. If in addition the internal energy $e(t,x)$ depends also on past random effects, then becomes $$\label{2014.6}
e(t,x)=\bar{\beta} u(t,x)+\int^t_0 n(t-s)u(s,x)ds+\int^t_0
\ell(t-s)h(s, u(s,x))dW_s,$$ where $W$ is a random process, such as Brownian motion, modeling the random effects. If $u(0, x)=0$, $\bar{\beta}=0$, $n(t)=\Gamma ( 1-\beta_1 )^{-1} t^{-\beta_1}$ and $\ell(t)=\Gamma (2- \beta_2 )^{-1}t^{1-\beta_2}$ for some constants $\beta_i \in (0,1)$, then (\[2014.6\]) (after differentiation in $t$) becomes $$\begin{aligned}
\label{e:1.5}
-{\rm div} \, \vec  F  = \frac{\partial e}{\partial t}(t,x)&=&\frac{1}{\Gamma (1- \beta_1 )}\frac{\partial }{\partial t} \int^t_0  (t-s)^{-\beta_1}    u(s,x)ds
\nonumber \\
&&+\frac1{\Gamma (2- \beta_2 )}\frac{\partial }{\partial t}\int^t_0
(t-s)^{1-\beta_2}h(s, u(s,x))dW_s.\end{aligned}$$ Since $$\begin{aligned}
\int_0^t (t-s)^{-{\beta_2}} \int_0^s
h(a, u(a,x))
dW_ads&= \int_0^t
 \int_a^t  (t-s)^{-{\beta_2}}ds
h(a, u(a,x))
dW_a
\\&=\frac{1}{1-{\beta_2}} \int_0^t
 (t-a)^{1-{\beta_2}}
h(a, u(a,x))
dW_a,\end{aligned}$$ by the definition of Caputo derivative (\[heat 1\]) we have $$\begin{aligned}
\partial^{\beta_2}_t \int_0^t
h(s, u(s,x))
dW_s&=
\frac{1}{\Gamma (1-{\beta_2})}\frac{\partial 