# First passage time for $g$–subdiffusion process of vanishing particles

Tadeusz Kosztołowicz<EMAIL_ADDRESS>Institute of Physics,
Jan Kochanowski University,
Uniwersytecka 7, 25-406 Kielce, Poland

###### Abstract

Subdiffusion equation and molecule survival equation, both with Caputo
fractional time derivatives with respect to another functions $g_{1}$ and
$g_{2}$, respectively, are used to describe diffusion of a molecule that can
disappear at any time with a constant probability. The process can be
interpreted as “ordinary” subdiffusion and “ordinary” molecule survival
process in which timescales are changed by the functions $g_{1}$ and $g_{2}$.
We derive the first-passage time distribution for the process. The mutual
influence of subdiffusion and molecule vanishing processes can be included in
the model when the functions $g_{1}$ and $g_{2}$ are related to each other. As
an example, we consider the processes in which subdiffusion and molecule
survival are highly related, which corresponds to the case of $g_{1}\equiv
g_{2}$.

Introduction. Subdiffusion occurs in media in which the movement of molecules
is very hindered due to a complex structure of a medium. Within the Continuous
Time Random Walk model the distribution of waiting time for a molecule to jump
has a heavy tail, $\tilde{\psi}(\tau)\sim 1/\tau^{\alpha+1}$, where
$\alpha\in(0,1)$ is a subdiffusion parameter mk ; mk1 ; ks . The parameter is
constant in a homogeneous system in which diffusion properties do not change
with time. The process can be described by the subdiffusion equation with the
“ordinary” Caputo fractional time derivative. Changing time scale may cause
subdiffusion to be accelerated or delayed. An example is the subordinate
method, in which a stochastic process changing this scale is involved sokolov1
; csm ; dybiec ; stan2020 . Recently, it has been proposed the subdiffusion
equation with a fractional time Caputo derivative with respect to another
function $g$ ($g$–subdiffusion equation) kd2022 ; kd2021a ; kd2021b . The
deterministic function $g$ rescales the time variable in the “ordinary”
subdiffusion model. For example, the waiting time for the molecule to jump in
$g$–subdiffusion model is described by the $\psi$ distribution which is
related to $\tilde{\psi}$ as $\psi(\tau)=\tilde{\psi}(g(\tau))$ kd2021b .

We study subdiffusion of a molecule that can be eliminated at any moment from
further diffusion, the molecule can be retained permanently or decayed; we
treat both processes as ”vanishing” of the molecule. We assume that the
probability density of a molecule disappearing does not depend on a molecule
position. Molecule elimination and diffusion processes can be related to each
other. An example is diffusion of an antibiotic through a bacterial biofilm.
Bacteria have different defense mechanisms against the action of the
antibiotic, see aot ; mot  and the references cited in km ; kmwa . One of them
is biofilm compaction, which makes it difficult for antibiotic molecules to
diffuse and even stops them. The permanent arrest of the molecule can be
treated as its ”disappearance” because it is eliminated from further
diffusion. Since the antibiotic mainly attacks rapidly growing bacteria,
another defense mechanism is to temporarily stop the bacteria from
multiplying. This mechanism even facilitates diffusion of the antibiotic
through the biofilm. Another example is subdiffusion of antibiotic molecules
in a medium having a ”plum pudding” structure. This process can occur when a
person infected with cystic fibrosis is also infected with Pseudomonas
aeruginosa bacteria kmwa . The ”pudding” represents the cystic fibrosis mucus,
and the ”plums” represent the Pseudomonas aeruginosa biofilm. The purpose of
the antibiotic is to kill the Pseudomonas aeruginosa bacteria whose defense
mechanisms can destroy the antibiotic molecules. In this case, slowing
subdiffusion in the ”pudding” may increase the chance of avoiding the decay of
the antibiotic molecules. These examples show that the processes of diffusion
and decay of molecules can be related to each other. Slowing down molecules
diffusion can cause an increase as well as a decrease of their disappearance
probability. Thus, relations between diffusion and processes that lead to a
molecules disappearance may be unique to a given process. Later we will pay
attention to the process in which subdiffusion and molecule survival are
highly related, and transport properties of a medium change periodically.

We consider the following problem. A molecule is located at point $x_{0}$ at
the initial time $t=0$ in a one–dimensional homogeneous system. We find the
probability distribution of time that a molecule reaches the point $x_{M}$ for
the first time. This problem is potentially of great importance in modeling of
the transport of antibiotics in bacterial biofilms and of epidemic spreading
when the arrival of an infected individual to some point may create a new
source for the spread of the epidemic.

Model. We assume that the disappearance of a molecule is independent of its
position. Let $P(x,t|x_{0})$ be a probability density (Green’s function) that
the molecule is at the point $x$ at time $t$ when particle decay would be
“turned off”, $x_{0}$ is the initial molecule position, and let $\rho(t)$ be a
probability that the molecule still exists at time $t$. Then, the Green’s
function $P_{\rho}$ for the vanishing molecule reads

$P_{\rho}(x,t|x_{0})=\rho(t)P(x,t|x_{0}).$ (1)

We also assume that in some initial time interval a diffusion process is
“ordinary” subdiffusion and the molecule vanishing process is described by the
“ordinary” decay equation. The “ordinary” subdiffusion is described by the
following equation

$\frac{{}^{C}\partial^{\alpha}P(x,t|x_{0})}{\partial
t^{\alpha}}=D\frac{\partial^{2}P(x,t|x_{0})}{\partial x^{2}},$ (2)

where ${}^{C}\partial^{\alpha}f(t)/\partial
t^{\alpha}=\int_{0}^{t}(t-u)^{-\alpha}f^{\prime}(u)du/\Gamma(1-\alpha)$ is the
“ordinary” Caputo fractional derivative, $0<\alpha<1$,
$f^{\prime}(t)=df(t)/dt$, $D$ is a subdiffusion coefficient given in the units
of $m^{2}/s^{\alpha}$.

A frequently used equation to describe the disappearance of a molecule is
$d\phi(t)/dt=-\lambda\phi(t)$, where $\phi$ is a probability density that
molecule decays at time $t$ and $\lambda$ is a decay rate. The solution to
this equation reads $\phi(t)=\lambda{\rm e}^{-\lambda t}$. The probability
that the molecule still exists at time $t$ is
$\rho(t)=1-\int_{0}^{t}\phi(u)du={\rm e}^{-\lambda t}$. This function fulfils
the molecule survival equation (MSE) $d\rho(t)/dt=-\lambda\rho(t)$ with the
initial condition $\rho(0)=1$. However, if the function $\rho$ has a heavy
tail, $\rho(t)\sim 1/t^{\beta}$ with $\beta\in(0,1)$ when
$t\rightarrow\infty$, the MSE involves “ordinary” fractional Caputo derivative
of the order $\beta$ suppl

$\frac{{}^{C}d^{\beta}\rho(t)}{dt^{\beta}}=-\lambda\rho(t),$ (3)

the unit of $\lambda$ is $1/s^{\beta}$.

We assume that at some initial stage the process is described by Eqs. (2) and
(3). In the further course of the process subdiffusion and the molecule
vanishing process may change due to the influence of external factors and/or
due to the interaction of the processes. In order to describe the processes in
the whole time domain, we assume that they are described by equations in which
“ordinary” Caputo derivatives are replaced by $g$–Caputo derivatives, see Eq.
(4) presented below. However, for short time the new equations should take the
form of Eqs. (2) and (3).

The $g$-Caputo fractional derivative ${}^{C}d^{\alpha}_{g}f(t)/dt^{\alpha}$ of
the order $\alpha$ with respect another function $g$ is defined for
$0<\alpha<1$ as almeida

$\frac{{}^{C}d^{\alpha}_{g}f(t)}{dt^{\alpha}}=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(g(t)-g(u))^{-\alpha}f^{\prime}(u)du,$
(4)

where $f^{\prime}(t)=df(t)/dt$, the function $g$ fulfils the conditions
$g(0)=0$, $g(\infty)=\infty$, and $g^{\prime}(t)>0$ for $t>0$. The values of
function $g$ are given in a time unit. When $g(t)=t$, the $g$-Caputo
fractional derivative takes the form of the “ordinary” Caputo derivative. For
$\alpha=1$, the $g$–Caputo derivative takes the form of the Riemann–Stieltjes
derivative ${}^{C}d^{1}_{g}f(t)/dt^{\alpha}\equiv
f^{[1]}_{g}=f^{\prime}(t)/g^{\prime}(t)$.

We involve the $g$–Caputo fractional time derivative in both the subdiffusion
equation and the molecular survival equation. We assume that $g$–subdiffusion
is controlled by a function $g_{1}$, and the $g$–molecule survival process by
a function $g_{2}$,

$\frac{{}^{C}\partial^{\alpha}_{g_{1}}P(x,t|x_{0})}{\partial
t^{\alpha}}=D\frac{\partial^{2}P(x,t|x_{0})}{\partial x^{2}},$ (5)
$\frac{{}^{C}d^{\beta}_{g_{2}}\rho(t)}{dt^{\beta}}=-\lambda\rho(t),$ (6)

$0<\alpha,\beta\leq 1$. As we have assumed, Eqs. (5) and (6) describe the
process over the entire time domain, but in the short time these equations
take the form of Eqs. (2) and (3), respectively. The latter condition is
fulfilled if $g_{i}(t)\rightarrow t$ when $t\rightarrow 0$. We assume that
$g_{i}(t)=t+\xi_{i}(t)$ with $\xi_{i}$ satisfying $\xi_{i}(t)\rightarrow 0$
when $t\rightarrow 0$, $i=1,2$. The functions $g_{1}$ and $g_{2}$ satisfy the
general conditions $g_{1}(0)=g_{2}(0)=0$,
$g_{1}(\infty)=g_{2}(\infty)=\infty$, and
$g_{1}(t),g_{2}(t),g^{\prime}_{1}(t),g^{\prime}_{2}(t)>0$ for $t>0$.

Solutions to Eqs. (5) and (6). The Green’s function $P$ is a solution to Eq.
(5) for the initial condition $P(x,0|x_{0})=\delta(x-x_{0})$, where $\delta$
is the delta–Dirac function, and for the boundary conditions
$P(\pm\infty,t|x_{0})=0$,

$P(x,t|x_{0})=\frac{1}{2\sqrt{D}}f_{-1+\alpha/2,\alpha/2}\left(g_{1}(t);\frac{|x-x_{0}|}{\sqrt{D}}\right),$
(7)

where

$\displaystyle
f_{\nu,\mu}(z;a)=\frac{1}{z^{1+\nu}}\sum_{k=0}^{\infty}\frac{1}{k!\Gamma(-\nu-\mu
k)}\left(-\frac{a}{z^{\mu}}\right)^{k},$ (8)

$a,\mu>0$. The solution to Eq. (6) for the initial condition $\rho(0)=1$ is

$\rho(t)={\rm E}_{\beta}\big{(}-\lambda g^{\beta}_{2}(t)\big{)},$ (9)

where ${\rm E}_{\beta}(z)=\sum_{n=0}^{\infty}z^{n}/\Gamma(\beta n+1)$ is the
Mittag–Leffler function, $0<\beta\leq 1$. We note that for $\beta=1$ we have
$\rho(t)={\rm e}^{-\lambda g_{2}(t)}$. From Eqs. (1), (7), and (9) we get

$\displaystyle P_{\rho}(x,t|x_{0})=\frac{{\rm E}_{\beta}\big{(}-\lambda
g^{\beta}_{2}(t)\big{)}}{2\sqrt{D}}$ (10) $\displaystyle\times
f_{-1+\alpha/2,\alpha/2}\left(g_{1}(t);\frac{|x-x_{0}|}{\sqrt{D}}\right).$

The calculations provided the above solutions are shown in suppl .

First passage time. The process of molecule vanishing is independent of its
position. We assume that the correlation of the functions $g_{1}$ and $g_{2}$,
if any, is not local and does not depend on the variable $x$. Let the system
consist of two parts $M_{1}$ and $M_{2}$ separated by a point $x_{M}$,
$M_{1}=(-\infty,x_{M})$ and $M_{2}=(x_{M},\infty)$. At the initial moment
$t=0$ the molecule is at the point $x_{0}$ located in the region $M_{1}$.

The probability that the molecule leaves the region $M_{1}=(-\infty,x_{M})$
first time in the time interval $(t,t+\Delta t)$, where $\Delta t$ is assumed
to be small, is $F(t;x_{0},x_{M})\Delta t=\rho(t)[R(t;x_{0},x_{M})-R(t+\Delta
t;x_{0},x_{M})]$. $R$ is a probability that the molecule did not leave the
$M_{1}$ region by time $t$,
$R(t;x_{0},x_{M})=\int_{-\infty}^{x_{M}}P_{abs}(x,t|x_{0})dx$, where
$P_{abs}(x,t|x_{0})$ is the probability of finding the molecule in the region
$M_{1}$ in a system with a fully absorbing wall located at $x_{M}$ redner .
The commonly used boundary condition at the absorbing wall is
$P_{abs}(x_{M},t|x_{0})=0$. The Green’s function for a system with a fully
absorbing wall can be found by means of the method of images, which for
$x,x_{0}<x_{M}$ gives $P_{abs}(x,t|x_{0})=P(x,t|x_{0})-P(x,t|2x_{M}-x_{0})$
feller ; chandra . After calculations we get

$\displaystyle
R(t;x_{0},x_{M})=1-f_{-1,\alpha/2}\left(g_{1}(t);\frac{x_{M}-x_{0}}{\sqrt{D}}\right).$
(11)

Taking the limit of $\Delta t\rightarrow 0$ we obtain

$F(t;x_{0},x_{M})=-\rho(t)\frac{dR(t;x_{0},x_{M})}{dt}.$ (12)

Eqs. (8), (11), and (12) provide

$\displaystyle F(t;x_{0},x_{M})=g_{1}^{\prime}(t){\rm
E}_{\beta}\big{(}-\lambda g^{\beta}_{2}(t)\big{)}$ (13) $\displaystyle\times
f_{0,\alpha/2}\left(g_{1}(t);\frac{x_{M}-x_{0}}{\sqrt{D}}\right).$

For $t\leq 0$ we put $F(t;x_{0},x_{M})\equiv 0$.

The probability $Q(t;x_{0},x_{M})$ that the molecule has passed the point
$x_{M}$ until time $t$ is calculated by means of the formula

$Q(t;x_{0},x_{M})=\int_{0}^{t}F(u;x_{0},x_{M})du.$ (14)

The function $Q$ will be calculated numerically. Another function describing
spreading of vanishing molecules is the mean square displacement of a particle
$\sigma_{\rho}^{2}(t)=\int_{-\infty}^{\infty}(x-\bar{x})^{2}P(x,t|x_{0})dx$,
where $\bar{x}$ is a mean value of $x$. For the Green’s function Eq. (10) we
get

$\sigma_{\rho}^{2}(t)=\rho(g_{2}(t))\frac{2D}{\Gamma(1+\alpha)}g_{1}^{\alpha}(t).$
(15)

When $g_{1}(t)\equiv t$ and $\rho(t)\equiv 1$, , the function
$\sigma_{\rho=1}^{2}(t)\sim t^{\alpha}$ is frequently used to define a kind of
diffusion; subdiffusion is when $0<\alpha<1$, normal diffusion is for
$\alpha=1$, and superdiffusion for $\alpha>1$.

Examples. As examples, we consider the case of fully correlated subdiffusion
and molecule vanishing processes. We assume that $g_{1}(t)\equiv
g_{2}(t)\equiv g(t)$ and that the function $\xi_{i}$ depends on a periodic
function. As $g$ we consider the following functions: $g_{0}(t)=t$,
$g_{A}(t)=t(1+a{\rm sin}(\tilde{\omega}(t)t+\phi))$, $g_{B}(t)=t+a((\omega
t)(\omega t-{\rm sin}(2\omega t))+(1-{\rm cos}(2\omega t))/2$, and
$g_{C}(t)=t+a({\rm sin}(\omega t+\phi)-{\rm sin}(\phi))$, $a$ and $\omega$ are
positive parameters. Since $g(t),g^{\prime}(t)>0$ we assume $a<1$ in the cases
of $A$ and $B$, and $a<1$, $a\omega<1$ in the case of $C$. In the case of $A$
the function $\tilde{\omega}$ depends on $t$, for the numerical calculations
we put $\tilde{\omega}(t)=t/(\tau_{0}+t)$, where $\tau_{0}$ is a parameter
given in a time unit.

Fig. 1 shows the plots of the first-passage time distribution Eq. (13). In
Fig. 2 there is presented the time evolution of the probability that the
molecule will reach the region $M_{2}$ by time $t$, which means that it has
passed the point $x_{M}$ at least once. The function $Q$ has been numerically
calculated from Eq. (14). Examples of the functions $\sigma_{\rho}^{2}$ Eq.
(15) and $\rho$ Eq. (9) are shown in Figs. 3 and 4.

Figure 1: Time evolution of the first-passage time distribution for
$g$–subdiffusion with molecule vanishing processes controlled by the functions
$g_{1}(t)\equiv g_{2}(t)\equiv g(t)$, where $g=g_{0}$ (symbol $O$ in the
legend), $g=g_{A}$ (symbol $A$), $g=g_{B}$ (symbol $B$), and $g=g_{C}$ (symbol
$C$). The parameters are $\alpha=0.6$, $\beta=0.5$, $\lambda=0.01$,
$\omega=1.5$, $\phi=0$, $a=0.4$, $x_{0}=0$, and $x_{M}=10$,
$\tilde{\omega}=1.5t/(1+t)$ for $g_{B}$, $\tau_{0}=1$, all quantities are
given in arbitrarily chosen units. Figure 2: Time evolution of the function
$Q$, the description is as in Fig. 1. Figure 3: The plots of function
$\sigma_{\rho}^{2}$ Eq. (15) for $\rho(t)\equiv 1$ (filled symbols) and for
$\rho(t)$ varying in time (open symbols) are compared, the further description
is as in Fig. 1. Figure 4: Time evolution of the probability that a particle
still exists $\rho^{2}(t)$, the description is as in Fig. 1.

The function $g_{0}$ does not change the timescale. In the cases $A$, $B$, and
$C$, the timescale is changed and periodic functions are involved in the
function $g$. In the case of $A$, $g_{A}$ contains ${\rm
sin}(\tilde{\omega}(t))$, the period of this function decreases with time and
goes to zero when $t\rightarrow\infty$. For the parameters given in the
caption of Fig. 1, $g_{A}$ does not generate noticeable oscillating effects,
see Fig. 1–4. For processes $B$ and $C$ oscillation effects are observed in
the figures. In process $B$, the amplitude of timescale oscillations increases
with time and is proportional to $t$. The change in the timescale is
relatively large, the rate of change is described by the function
$g^{\prime}_{B}(t)=1+at{\rm sin}^{2}(\omega t)$, which results in $g_{B}\sim
t^{2\alpha}$ in the long time limit. Thus, assuming $\alpha>0.5$, the
oscillation effect is superimposed here on the superdiffusion effect. In the
case of $C$, the amplitude of oscillations are constant, the oscillation
effect decreases with time. Here we have $\sigma_{\rho}^{2}(t)\sim t^{\alpha}$
when $t\rightarrow\infty$, there is the “ordinary” subdiffusion effect in the
long time limit.

Final remarks. Subdiffusion and a molecule vanishing process are described by
Eqs. (5) and (6). In general, involving $g$–Caputo fractional time derivatives
in a model, instead of “ordinary” Caputo derivatives, we can describe changes
of the course of the processes caused by external factors (e.g. temperature
change causing the temporal evolution of parameters), and by the influence of
different processes on each other. It has already been shown that the
$g$–subdiffusion equation describe subdiffusion with a change of parameters
$\alpha$ and $D$ kd2022  and a transition from ”ordinary” subdiffusion to slow
subdiffusion (ultraslow diffusion) kd2021a . In this paper, we use the
equations with $g$–Caputo fractional derivatives to model subdiffusion and
molecule vanishing process which can be related to each other.

We have assumed the molecule vanishing has a constant rate $\lambda$ and does
not depend on the molecule position. This process can be interpreted as
subdiffusion with $A\rightarrow\emptyset$ reaction. In kl  it was shown that
the form of a subdiffusion–reaction equation depends on whether a diffusing
molecule $A$ can decay with a probability independent of its position, or
whether it must meet a $B$ molecule, with a certain probability, for a
reaction to occur. Subdiffusion-reaction equations are qualitatively different
in both cases kl ; sokolov . If the function $\rho$ also depended on the
variable $x$ then the process could be described by a following
$g$–subdiffusion–reaction equation with a reaction $nA+mB\rightarrow\emptyset$
($n$ and $m$ are positive parameters),

$\displaystyle\frac{{}^{C}\partial^{\alpha}_{g_{1}}P_{\rho}(x,t|x_{0})}{\partial
t^{\alpha}}=D\frac{\partial^{2}P_{\rho}(x,t|x_{0})}{\partial x^{2}}$ (16)
$\displaystyle-\Xi(P_{\rho},C_{A}(x,t),C_{B}(x,t)),$

where $P_{\rho}$ describes subdiffusion of $A$ molecule, $\Xi$ is a reaction
term which can be controlled by both functions $g_{1}$ and $g_{2}$, $C_{A}$,
$C_{B}$ denote concentrations of $A$ and $B$ molecules respectively. We note
that if $n\neq 1$, a probability of the reaction depends on the concentration
$C_{A}$. The reaction term may be non-linear, the analysis of a general
$g$–subdiffusion–reaction equation will be presented elsewhere.
$G$–subdiffusion with $A\rightarrow\emptyset$ reaction can be described by the
subdiffusion equation and the molecule survival equation using the model
presented in this paper.

## References

  * (1) R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339, 1-77 (2000).
  * (2) R. Metzler and J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37, R161-R208 (2004).
  * (3) J. Klafter and I.M. Sokolov, First step in random walks. From tools to applications (Oxford UP, New York, 2011).
  * (4) I.M. Sokolov, Thermodynamics and fractional Fokker-Planck equations, Phys. Rev. E 63, 05611 (2001).
  * (5) A.V. Chechkin, F. Seno, R. Metzler, and I.M. Sokolov, Brownian yet non-Gaussian diffusion: From superstatistics to subordination of diffusing diffusivities, Phys. Rev. X 7, 021002 (2017).
  * (6) B. Dybiec and E. Gudowska–Nowak, Subordinated diffusion and continuous time random walk asymptotics, Chaos 20, 043129 (2010).
  * (7) A. Stanislavsky and A. Weron, Accelerating and retarding anomalous diffusion: A Bernstein function approach, Phys Rev. E 101, 052119 (2020).
  * (8) T. Kosztołowicz and A. Dutkiewicz, $G$-subdiffusion equation that describes transient subdiffusion, arXiv:2205.11625 (2022).
  * (9) T. Kosztołowicz and A. Dutkiewicz, Subdiffusion equation with Caputo fractional derivative with respect to another function, Phys. Rev. E 104, 014118 (2021).
  * (10) T. Kosztołowicz and A. Dutkiewicz, Stochastic interpretation of $g$-subdiffusion process, Phys. Rev. E 104, L042101 (2021).
  * (11) G.G. Anderson and G.A. O’Toole, Bacterial Biofilms, Current Topics in Microbiology and Immunology 322, p. 85 (Berlin, Springer, 2008).
  * (12) T.F.C. Mah and G.A. O’Toole, Mechanisms of biofilm resistance to antimicrobial agents, Trends Microbiol. 9, 34-39 (2001).
  * (13) T. Kosztołowicz and R. Metzler, Diffusion of antibiotics through a biofilm in the presence of diffusion and absorption barriers, Phys. Rev. E 102, 032408 (2020).
  * (14) T. Kosztołowicz, R. Metzler, S. Wa̧sik, and M. Arabski, Modelling experimentally measured of ciprofloxacin antibiotic diffusion in Pseudomonas aeruginosa biofilm formed in artificial sputum medium, PLoS One 15(12), e0243003 (2020).
  * (15) T. Kosztołowicz, How to solve equations with the $g$–Caputo fractional derivative, Supplemental material.
  * (16) R. Almeida, A Caputo fractional derivative of a function with respect to another function, Commun. Nonlinear Sci. Numer. Simulat. 44, 460-481 (2017).
  * (17) S. Redner, A guide to first passage processes (Cambridge UP, Cambridge, 2001).
  * (18) W. Feller, An introduction to probability theory and its applications, Vol. 2 (Wiley, New York, 1968).
  * (19) S. Chandrasekhar, Stochastic problems in physics and astronomy, Rev. Mod. Phys. 15, 1-89 (1943).
  * (20) T. Kosztołowicz and K.D. Lewandowska, Subdiffusion-reaction processes with $A\rightarrow B$ reactions versus subdiffusion-reaction processes with $A+B\rightarrow B$ reactions, Phys. Rev. E 90, 032136 (2014).
  * (21) I.M. Sokolov, M.G.W. Schmidt, and F. Sagués, Reaction–subdiffusion equation, Phys. Rev. E 73, 031102 (2001).