Feynman-Kac Equations for Reaction and Diffusion Processes

TL;DR

This paper develops a comprehensive theoretical framework for deriving Feynman-Kac equations that describe the distribution of functionals of particles undergoing diffusion and chemical reactions, applicable to various diffusion types and reaction rates.

Contribution

It introduces general forms of forward and backward Feynman-Kac equations for reaction-diffusion processes, including normal and anomalous diffusion with linear and nonlinear reactions.

Findings

Derived equations for different diffusion and reaction types

Analyzed occupation time and first passage time using the equations

Included boundary conditions like absorbing and reflecting

Abstract

This paper provides a theoretical framework of deriving the forward and backward Feynman-Kac equations for the distribution of functionals of the path of a particle undergoing both diffusion and chemical reaction. Very general forms of the equations are obtained. Once given the diffusion type and reaction rate, a specific forward or backward Feynman-Kac equation can be obtained. The listed in the paper include the ones for normal/anomalous diffusions and reactions with linear/nonlinear rates. Using the derived equations, we also study the occupation time in half-space, the first passage time to a fixed boundary, and the occupation time in half-space with absorbing or reflecting boundary conditions.

Equations

$$\frac{\partial \rho (x,t)}{\partial t}=K_1\frac{{\partial }^2\rho (x,t)}{\partial x^2}+r(\rho )\rho ,$$

$$\frac{\partial \rho (x,t)}{\partial t}=K_{\alpha }{e}^{-\kappa t}{\mathbb{D}}_t^{1-\alpha }\left(e^{\kappa t}\frac{{\partial }^2\rho (x,t)}{\partial x^2}\right)-\kappa \rho (x,t),$$

$$\rho (x,A,t)={\rho }_0(x)\delta \left(A-{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }\right)e^{{\int }_0^{t}r(\rho (x,u))du}\Psi (t)+{\int }_0^{t}j\left(x,A-{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime },\tau \right)e^{{\int }_{\tau }^{t}r(\rho (x,u))du}\Psi (t-\tau )d\tau ,$$

$$\begin{array}{rl}j(x,A,t)& ={\int }_{\mathbb{R}}{\rho }_0(x-z)\delta \left(A-{\int }_0^{t}U(x-z,t^{\prime })d{t}^{\prime }\right)e^{{\int }_0^{t}r(\rho (x-z,u))du}w(z)\varphi (t)dz\\ & +{\int }_0^t{\int }_{\mathbb{R}}j\left(x-z,A-{\int }_{\tau }^{t}U(x-z,t^{\prime })d{t}^{\prime },\tau \right)e^{{\int }_{\tau }^{t}r(\rho (x-z,u))du}w(z)\varphi (t-\tau )dzd\tau .\end{array}$$

$$\rho (x,p,t)={\rho }_0(x)e^{ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_0^{t}r(\rho (x,u))du}\Psi (t)+{\int }_0^{t}j(x,p,\tau )e^{ip{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(\rho (x,u))du}\Psi (t-\tau )d\tau ,$$

$$\begin{array}{rl}j(x,p,t)& ={\int }_{\mathbb{R}}{\rho }_0(x-z)e^{ip{\int }_0^{t}U(x-z,t^{\prime })d{t}^{\prime }}{e}^{{\int }_0^{t}r(\rho (x-z,u))du}w(z)\varphi (t)dz\\ & +{\int }_0^t{\int }_{\mathbb{R}}j(x-z,p,\tau )e^{ip{\int }_{\tau }^{t}U(x-z,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(\rho (x-z,u))du}w(z)\varphi (t-\tau )dzd\tau .\end{array}$$

$$\frac{\partial \rho (x,p,t)}{\partial t}=[ipU(x,t)+r(\rho (x,t))]\rho (x,p,t)+j(x,p,t)-i(x,p,t),$$

$$i(x,p,t)={\rho }_0(x)e^{ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_0^{t}r(\rho (x,u))du}\varphi (t)+{\int }_0^{t}j(x,p,\tau )e^{ip{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(\rho (x,u))du}\varphi (t-\tau )d\tau .$$

$$j(x,p,t)={\int }_{\mathbb{R}}i(x-z,p,t)w(z)dz.$$

$$\mathcal{L}\{\rho (x,p,t)e^{-ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{t}r(\rho (x,u))du}\}=[{\rho }_0(x)+\mathcal{L}\{j(x,p,t)e^{-ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{t}r(\rho (x,u))du}\}]\tilde{\Psi }(s),$$

$$\mathcal{L}\{i(x,p,t)e^{-ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{t}r(\rho (x,u))du}\}=[{\rho }_0(x)+\mathcal{L}\{j(x,p,t)e^{-ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{t}r(\rho (x,u))du}\}]\tilde{\varphi }(s),$$

$$\mathcal{L}\{i(x,p,t)e^{-ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{t}r(\rho (x,u))du}\}=\frac{\tilde{\varphi }(s)}{\tilde{\Psi }(s)}\mathcal{L}\{\rho (x,p,t)e^{-ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{t}r(\rho (x,u))du}\}.$$

$$i(x,p,t)={\int }_0^{t}K(t-\tau )\rho (x,p,\tau )e^{ip{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(\rho (x,u))du}d\tau ,$$

$$\tilde{K}(s)=\frac{\tilde{\varphi }(s)}{\tilde{\Psi }(s)}=\frac{s\tilde{\varphi }(s)}{1-\tilde{\varphi }(s)}.$$

$$\begin{array}{rl}\frac{\partial \rho (x,p,t)}{\partial t}& ={\int }_{\mathbb{R}}{\int }_0^{t}K(t-\tau )\rho (x-z,p,\tau )e^{ip{\int }_{\tau }^{t}U(x-z,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(\rho (x-z,u))du}w(z)d\tau dz\\ & -{\int }_0^{t}K(t-\tau )\rho (x,p,\tau )e^{ip{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(\rho (x,u))du}d\tau +[ipU(x,t)+r(\rho (x,t))]\rho (x,p,t),\end{array}$$

$$\begin{array}{rl}\frac{\partial n(x,p,t)}{\partial t}& ={\int }_{\mathbb{R}}{\int }_0^{t}K(t-\tau )n(x-z,p,\tau )e^{ip{\int }_{\tau }^{t}U(x-z,t^{\prime })d{t}^{\prime }}w(z)d\tau dz\\ & -{\int }_0^{t}K(t-\tau )n(x,p,\tau )e^{ip{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime }}d\tau +ipU(x,t)n(x,p,t).\end{array}$$

$$\frac{\partial \rho (x,p,t)}{\partial t}=\frac{{\sigma }^2}{2}\frac{{\partial }^2}{\partial x^2}{\int }_0^{t}K(t-\tau )\rho (x,p,\tau )e^{ip{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(\rho (x,u))du}d\tau +[ipU(x,t)+r(\rho (x,t))]\rho (x,p,t).$$

$$\begin{array}{rl}\frac{\partial \rho (x,p,t)}{\partial t}& =\frac{{\sigma }^2}{2\tau }\frac{{\partial }^2\rho (x,p,t)}{\partial x^2}+ipU(x,t)\rho (x,p,t)\\ & +r(\rho (x,t))\rho (x,p,t),\end{array}$$

$$\begin{array}{rl}\frac{\partial \rho (x,p,t)}{\partial t}& =\frac{{\sigma }^2}{2{\tau }^{\alpha }\Gamma (1-\alpha )}\frac{{\partial }^2}{\partial x^2}\left[e^{ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_0^{t}r(\rho (x,u))du}{\mathbb{D}}_t^{1-\alpha }\{\rho (x,p,t)e^{-ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{t}r(\rho (x,u))du}\}\right]\\ & +[ipU(x,t)+r(\rho (x,t))]\rho (x,p,t),\end{array}$$

$$\begin{array}{rl}\frac{\partial \rho (x,p,t)}{\partial t}& =\frac{{\sigma }^2\alpha C_{\tau }}{2\Gamma (1-\alpha )}\frac{{\partial }^2}{\partial x^2}\left[e^{ip{\int }_0^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_0^{t}r(\rho (x,u))du}\frac{\partial }{\partial t}{\int }_0^{t}H(t-\tau )\rho (x,p,\tau )e^{-ip{\int }_0^{\tau }U(x,t^{\prime })d{t}^{\prime }}{e}^{-{\int }_0^{\tau }r(\rho (x,u))du}d\tau \right]\\ & +[ipU(x,t)+r(\rho (x,t))]\rho (x,p,t),\end{array}$$

$$\begin{array}{rl}\frac{\partial \rho (x,p,t)}{\partial t}& ={\int }_{\mathbb{R}}{\int }_0^{t}K(t-\tau )\rho (x-z,p,\tau )e^{ip{\int }_{\tau }^{t}U(x-z,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(u)du}w(z)d\tau dz\\ & -{\int }_0^{t}K(t-\tau )\rho (x,p,\tau )e^{ip{\int }_{\tau }^{t}U(x,t^{\prime })d{t}^{\prime }}{e}^{{\int }_{\tau }^{t}r(u)du}d\tau +[ipU(x,t)+r(t)]\rho (x,p,t),\end{array}$$

$$\rho (x,A,t)={\rho }_0(x)\delta (A-tU(x))e^{r(x)t}\Psi (t)+{\int }_0^{t}j(x,A-(t-\tau )U(x),\tau )e^{(t-\tau )r(x)}\Psi (t-\tau )d\tau ,$$

$$\begin{array}{rl}j(x,A,t)& ={\int }_{\mathbb{R}}{\rho }_0(x-z)\delta (A-tU(x-z))e^{r(x-z)t}w(z)\varphi (t)dz\\ & +{\int }_0^t{\int }_{\mathbb{R}}j(x-z,A-(t-\tau )U(x-z),\tau )e^{(t-\tau )r(x-z)}w(z)\varphi (t-\tau )dzd\tau .\end{array}$$

$$\rho (x,p,t)={\rho }_0(x)e^{iptU(x)}{e}^{r(x)t}\Psi (t)+{\int }_0^{t}j(x,p,\tau )e^{ip(t-\tau )U(x)}{e}^{(t-\tau )r(x)}\Psi (t-\tau )d\tau ,$$

$$\begin{array}{rl}j(x,p,t)& ={\int }_{\mathbb{R}}{\rho }_0(x-z)e^{iptU(x-z)}{e}^{r(x-z)t}w(z)\varphi (t)dz\\ & +{\int }_0^t{\int }_{\mathbb{R}}j(x-z,p,\tau )e^{ip(t-\tau )U(x-z)}{e}^{(t-\tau )r(x-z)}w(z)\varphi (t-\tau )dzd\tau .\end{array}$$

$$\rho (x,p,s)={\rho }_0(x)\tilde{\Psi }(s-ipU(x)-r(x))+j(x,p,s)\tilde{\Psi }(s-ipU(x)-r(x)),$$

$$j(x,p,s)={\int }_{\mathbb{R}}{\rho }_0(x-z)w(z)\tilde{\varphi }(s-ipU(x-z)-r(x-z))dz+{\int }_{\mathbb{R}}j(x-z,p,s)w(z)\tilde{\varphi }(s-ipU(x-z)-r(x-z))dz.$$

$$\rho (k,p,s)=\tilde{\Psi }\left(s-ipU\left(-i\frac{\partial }{\partial k}\right)-r\left(-i\frac{\partial }{\partial k}\right)\right)[{\rho }_0(k)+j(k,p,s)],$$

$$j(k,p,s)=\hat{w}(k)\tilde{\varphi }\left(s-ipU\left(-i\frac{\partial }{\partial k}\right)-r\left(-i\frac{\partial }{\partial k}\right)\right)[{\rho }_0(k)+j(k,p,s)],$$

$$\rho (k,p,s)=\frac{\tilde{\Psi }\left(s-ipU\left(-i\frac{\partial }{\partial k}\right)-r\left(-i\frac{\partial }{\partial k}\right)\right){\rho }_0(k)}{1-\hat{w}(k)\tilde{\varphi }\left(s-ipU\left(-i\frac{\partial }{\partial k}\right)-r\left(-i\frac{\partial }{\partial k}\right)\right)}.$$

Full text

††preprint: APS/123-QED

I Introduction

Diffusion is the net movement of particles from a region of high concentration to a region of low concentration. It is called normal diffusion if the ensemble mean squared displacement of the diffusing particles scales linearly, otherwise it is called anomalous diffusion, i.e., $\langle x^{2}(t)\rangle\sim t^{\mu}$ with $\mu\neq 1$, being subdiffusion for $\mu\in(0,1)$ and superdiffusion for $\mu\in(1,2)$ Magdziarz2007 ; Magdziarz2008 ; Cairoli2015 . For pure diffusion processes (with external forces), the equations governing the distribution of functionals of the path of particle are well established Cairoli2015a ; Cairoli2017 . Influenced by Feynman’s thesis of the derivation of Schrodinger’s equation, in 1949, Kac derives the classical Feynman-Kac equation for normal diffusion Kac1949 . Turgeman et al. present the fractional Feynman-Kac equation for non-Brownian functionals Turgeman2009 and tempered fractional Feynman-Kac equation for tempered non-Brownian functionals is provided in Wu2016 .

Chemical reaction leads to the transformation of one type of particles to another type. The main concern of this paper is on the distribution of the functional of the paths of a particle performing both diffusion and chemical reaction. We first recall the governing equation for the probability density function (PDF) $\rho(x,t)$ of the reaction diffusion process. Assume that the reaction rate law is $\dot{\rho}=r(\rho)\rho$. For normal diffusion, the well-known reaction diffusion equation is

[TABLE]

with $K_{1}$ being the diffusion constant. For anomalous diffusion, the situation becomes complicated, and the resulted equation can not be obtained by simply adding the reaction term into the non-Markovian transport equation. For example, in Sokolov2006 ; Henry2006 , the authors derive the subdiffusive reaction diffusion equation:

[TABLE]

where $K_{\alpha}$ is the generalized diffusion coefficient and $\mathbb{D}_{t}^{1-\alpha}$ is the Riemann-Liouville fractional derivative. For the multispecies system with anomalous diffusion and reaction dynamics, refer to Langlands2008 ; Fedotov2010a ; Froemberg2011 , and for the front propagation problems, see Yadav2007 ; Iomin2012 ; Iomin2013 .

Functionals of paths of particles have applications in rather diverse fields, including probability theory, the KPZ varieties, finance, economics, even in characterizing the stochastic behaviour of daily temperature records Majumdar2005 . So far, it seems there are no research works on the functional distribution of the reaction diffusion processes. Moreover, the functional in all published works is defined as the form $A=\int_{0}^{t}U(x(\tau))d\tau$, where $x(t)$ is the particle trajectory and $U(x)$ is a specified function. Here we extend it as $A=\int_{0}^{t}U(x(\tau),\tau)d\tau$, which may be used to model the first passage time to a moving boundary Tuckwell1984 . As for the chemical reactions, both nonlinear reaction rate $r(\rho(x,t))$ and linear reaction rates $r(t)$ and $r(x)$ are considered, i.e., the equations governing reaction process are, respectively, $\dot{\rho}=r(\rho)\rho$, $\dot{\rho}=r(t)\rho$, or $\dot{\rho}=r(x)\rho$. The framework Turgeman2009 for deriving the Feynman-Kac equations for pure diffusion process does not work for the reaction diffusion process. So, we present new theoretical frameworks of deriving the forward and backward Feynman-Kac equations for the distribution of the more general functionals of the path of a particle undergoing both diffusion and chemical reaction with linear/nonlinear rates; the new framework also applies to the pure diffusion process and the functional $A=\int_{0}^{t}U(x(\tau))d\tau$. Based on the very general forms of the derived forward and backward Feynman-Kac equations, we more specifically provide the equations for normal/anomalous diffusions and reactions with linear/nonlinear rate. And the derived equations are used to calculate the occupation time in half-space, the first passage time to a fixed boundary, and the occupation time in half-space with absorbing or reflecting boundary conditions.

The remainder of this paper is organized as follows. In Section II, we derive the general forward Feynman-Kac equations for the process with nonlinear reaction rate $r(\rho(x,t))$, and discuss some important properties and examples for the case of linear reaction rates $r(t)$ and $r(x)$. The backward Feynman-Kac equations are presented in Section III. Following the derived equations, in Section IV, some specific functional distributions are analytically solved, including the occupation time in half-space, the first passage time to a fixed boundary and the occupation time in half-interval with absorbing or reflecting boundary conditions, and some numerical simulations are performed to verify the analytical results. We conclude the paper with some discussions in Section V.

II Derivation of the Forward Feynman-Kac Equations

Now, we start to derive the forward Feynman-Kac equations for the reaction diffusion processes, composed of three subsections, respectively, for nonlinear reaction rate $r(\rho(x,t))$, linear reaction rate $r(t)$, and linear reaction rate $r(x)$. We limit to the case that all the reaction rates are nonpositive.

II.1 Nonlinear reaction rate $r(\rho(x,t))$

Let $\phi(t)$ and $w(x)$ be the waiting time and symmetric jump length PDF, respectively. We define a stochastic process $\{x(t),A(t)\}$ describing the time-varying positions and functional values that makes physical sense for the reaction diffusion processes. Denote $\rho(x,A,t)$ as the joint PDF of finding the particle at position $x$ with the functional value $A$ at time $t$. We extend the method in Fedotov2010 to formulate the balance equations for the density $\rho(x,A,t)$ and the auxiliary density $j(x,A,t)$. The latter is the joint PDF of particles exactly arriving at position $x$ with the functional value $A$ at time $t$. The balance equation for $\rho(x,A,t)$ can be written as

[TABLE]

where $\Psi(t)=\int_{t}^{+\infty}\phi(\tau)d\tau$ is the survival probability function, i.e., the probability that the particles keep waiting up to time $t$. The first term on the right hand side (RHS) of Eq. (3) represents the particles that stay at their initial position $x$ up to time $t$. The corresponding functional value of these motionless particles can be calculated as $A=\int_{0}^{t}U(x,t^{\prime})dt^{\prime}$. Hence the Dirac delta function $\delta(\cdot)$ is introduced. $r(\rho(x,t))$ denotes the nonlinear reaction rate with $\rho(x,t)=\int_{-\infty}^{+\infty}\rho(x,A,t)dA$. The second term on the RHS of Eq. (3) describes the particles that arrive at position $x$ at some previous time $\tau<t$ and don’t move during the time interval $[\tau,t]$. The difference of the functional values between the time interval $[0,\tau]$ and $[0,t]$ is $\int_{\tau}^{t}U(x,t^{\prime})dt^{\prime}$, since no jumps take place during $[\tau,t]$. If the usual functional $A=\int_{0}^{t}U(x(t^{\prime}))dt^{\prime}$ is used, one just needs to replace $\int_{0}^{t}U(x,t^{\prime})dt^{\prime}$ and $\int_{\tau}^{t}U(x,t^{\prime})dt^{\prime}$ in Eq. (3) with $tU(x)$ and $(t-\tau)U(x)$, respectively.

Moreover, the balance equation for $j(x,A,t)$ can be written in the following form:

[TABLE]

The first term on the RHS of Eq. (4) expresses the particles that stay at the initial position $x-z$ and jump to $x$ instantaneously at time $t$. The second term demonstrates the particles that arrive at the position $x-z$ at some previous time $\tau<t$, and react up to time $t$ and jump to $x$ exactly at time $t$.

In order to derive the forward Feynman-Kac equation, we need to conduct the Fourier transforms $A\rightarrow p$ (denoted by $\mathcal{F}\{f(A)\}=\int_{-\infty}^{+\infty}e^{ipA}f(A)dA$) for both Eqs. (3) and (4), since $A$ could be negative. Then, we obtain

[TABLE]

and

[TABLE]

Since Eqs. (5) and (6) have nonlinear terms, the standard Fourier-Laplace transform techniques used in Turgeman2009 ; Carmi2010 ; Carmi2011 ; Wu2016 cannot be applied directly. Instead, we differentiate the density $\rho(x,p,t)$ given by Eq. (5) with respect to $t$, which leads to

[TABLE]

where

[TABLE]

It follows from Eqs. (6) and (8) that

[TABLE]

In order to express $i(x,p,t)$ in terms of $\rho(x,p,t)$, we divide Eqs. (5) and (8) by the same factor $e^{ip\int_{0}^{t}U(x,t^{\prime})dt^{\prime}}e^{\int_{0}^{t}r(\rho(x,u))du}$ and then take the Laplace transform $t\rightarrow s$ (denoted by $\mathcal{L}\{f(t)\}=\int_{0}^{+\infty}e^{-st}f(t)dt$) for both equations

[TABLE]

[TABLE]

where $\tilde{\Psi}(s)=\mathcal{L}\{\Psi(t)\}$ and $\tilde{\phi}(s)=\mathcal{L}\{\phi(t)\}$. From Eqs. (10) and (11), we obtain

[TABLE]

Performing inverse Laplace transform on Eq. (12) gives

[TABLE]

where $K(t)$ is the memory kernel defined by its Laplace transform Fedotov2010

[TABLE]

Substituting Eq. (13) into Eqs. (7) and (9) gives the general forward Feynman-Kac equation with nonlinear reaction rate $r(\rho(x,t))$, i.e.,

[TABLE]

where $K(t)$ is defined by Eq. (14). This generalized Feynman-Kac equation has few restrictions on the forms of the waiting time and jump length densities, which allows us to conveniently model standard transport, anomalous diffusion or tempered dynamics together with nonlinear reaction terms.

Obviously, plugging $r(\rho(x,t))\equiv 0$ into Eq. (15) gives the forward Feynman-Kac equation for the pure diffusion process without reactions. To avoid confusion, we denote this equation in terms of $n(x,p,t)$, namely,

[TABLE]

Assume that the symmetric jump length distribution $w(z)$ satisfies $\int_{\mathbb{R}}z^{2}w(z)dz=\sigma^{2}<+\infty$. Expanding the first term on the RHS of Eq. (15) for small $z$ and truncating the Taylor series at the second order (see Appendix A for more details) gives

[TABLE]

Now we consider several special cases of Eqs. (15) and (17). To start with, if $p=0$, then $\rho(x,0,t)=\int_{-\infty}^{+\infty}e^{ipA}\rho(x,A,t)dA=\rho(x,t)$. Thus, when $p=0$, Eqs. (15) and (17) reduce to Eqs. (14) and (22) in Fedotov2010 respectively, which describe the nonlinear master equations for reaction diffusion processes. Furthermore, we can specify the waiting time distributions and write down the corresponding Feynman-Kac equations in the continuum limit of large scales and long time as follows.

(a) Exponential waiting time distribution

Suppose $\phi(t)=\frac{1}{\tau}\exp(-\frac{t}{\tau})$. Thus, $\tilde{\phi}(s)=\int_{0}^{+\infty}e^{-st}\phi(t)dt=\frac{1}{\tau s+1}$ and $\tilde{K}(s)=\frac{s\tilde{\phi}(s)}{1-\tilde{\phi}(s)}=\frac{1}{\tau}$. Inserting $K(t)=\frac{1}{\tau}\delta(t)$ into Eq. (17) gives

[TABLE]

which is the standard Feynman-Kac equation for the normal diffusion (Markovian) process with the reaction rate $r(\rho(x,t))$.

(b) Power-law waiting time distribution

When $\phi(t)=\alpha\tau^{\alpha}t^{-1-\alpha}\mathbf{1}_{[\tau,+\infty)}(t)$ with $0<\alpha<1$, we can calculate (see Appendix B) $\tilde{\phi}(s)$, which has the asymptotic form $\tilde{\phi}(s)\sim 1-\tau^{\alpha}\Gamma(1-\alpha)s^{\alpha}$ for $s\rightarrow 0$. Meanwhile, according to Eq. (14), $\tilde{K}(s)\approx\frac{s^{1-\alpha}}{\tau^{\alpha}\Gamma(1-\alpha)}$. Thus, Eq. (17) reduces to

[TABLE]

where $\mathbb{D}_{t}^{1-\alpha}$ is the Riemann-Liouville fractional derivative defined as $\mathbb{D}_{t}^{1-\alpha}\rho(x,p,t)=\frac{1}{\Gamma(\alpha)}\frac{\partial}{\partial t}\int_{0}^{t}\frac{\rho(x,p,t^{\prime})}{(t-t^{\prime})^{1-\alpha}}dt^{\prime}$.

(c) Tempered power-law waiting time distribution

Assume $\phi(t)=C_{\tau}^{-1}e^{-\lambda t}t^{-1-\alpha}\mathbf{1}_{[\tau,+\infty)}(t)$, where $0<\alpha<1$ and $C_{\tau}=\int_{\tau}^{+\infty}e^{-\lambda t}t^{-1-\alpha}dt$. Then it can be calculated (see also Appendix B) that $\tilde{\phi}(s)$ has the asymptotic form $\tilde{\phi}(s)\sim 1-\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}(s+\lambda)^{\alpha}+\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}\lambda^{\alpha}$ for $s\rightarrow 0$. Similarly, $\tilde{K}(s)\approx\frac{\alpha C_{\tau}}{\Gamma(1-\alpha)}\frac{s}{(s+\lambda)^{\alpha}-\lambda^{\alpha}}$. Under this condition, Eq. (17) can be rewritten as

[TABLE]

where $H(t)=e^{-\lambda t}t^{\alpha-1}E_{\alpha,\alpha}(\lambda^{\alpha}t^{\alpha})$ and its Laplace transform $\tilde{H}(s)=\frac{1}{(s+\lambda)^{\alpha}-\lambda^{\alpha}}$. Here $E_{\alpha,\alpha}(\cdot)$ is the two-parameter Mittag-Leffler function Podlubny1999 , which is defined as $E_{\alpha,\beta}(z)=\sum\limits_{k=0}^{\infty}\frac{z^{k}}{\Gamma(\alpha k+\beta)}$ with $\alpha>0$ and $\beta>0$ .

II.2 Linear reaction rate $r(t)$

Following the similar derivation procedure to the previous subsection, we obtain the forward Feynman-Kac equation with linear reaction rate $r(t)$, i.e.,

[TABLE]

where $K(t)$ is defined by Eq. (14).

By simple calculations, it can be easily checked that $\rho(x,p,t)=n(x,p,t)e^{\int_{0}^{t}r(u)du}$ satisfies Eq. (21) if $n(x,p,t)$ is the solution Eq. (16), which means that for the reaction diffusion process with linear reaction rate $r(t)$, the effect of transport with memory and the linear reaction rate dependening on time $t$ can be decoupled. In fact, this also implies that the conclusion of Eqs. (15-17) in Fedotov2010 for the case of constant reaction rate still holds for the linear reaction rate, and even for describing the distribution of functionals.

II.3 Linear reaction rate $r(x)$

We further consider the linear reaction rate $r(x)$ and the traditional functional definition $A=\int_{0}^{t}U(x(\tau))d\tau$, in order to give another derivation of the corresponding forward Feynman-Kac equation. To start with, under these conditions, we have

[TABLE]

and

[TABLE]

Performing Fourier transforms $A\rightarrow p$ on both sides of Eq. (22) and Eq. (23) gives that

[TABLE]

and

[TABLE]

Conducting Laplace transform $t\rightarrow s$ to Eqs. (24) and (25), respectively, we obtain

[TABLE]

and

[TABLE]

Conducting Fourier transforms $x\rightarrow k$ to Eqs. (26) and (27) gives

[TABLE]

and

[TABLE]

since we have the identity Carmi2010 $\mathcal{F}\{xf(x)\}=-i\frac{\partial}{\partial k}\hat{f}(k)$. Substituting $j(k,p,s)=\frac{\hat{w}(k)\tilde{\phi}\left(s-ipU\left(-i\frac{\partial}{\partial k}\right)-r\left(-i\frac{\partial}{\partial k}\right)\right)\rho_{0}(k)}{1-\hat{w}(k)\tilde{\phi}\left(s-ipU\left(-i\frac{\partial}{\partial k}\right)-r\left(-i\frac{\partial}{\partial k}\right)\right)}$ into Eq. (28), we get

[TABLE]

Furthermore, we specify three kinds of typical waiting time and jump length distributions and then calculate their asymptotic forms in the Laplace or Fourier domain as $s\rightarrow 0$ or $k\rightarrow 0$. Naturally, the inverse Fourier-Laplace transform technique is applied to write down the respective forward Feynman-Kac equations in the long time and large scales limits.

(a) Exponential waiting time and Gaussian jump length distributions

Suppose $\phi(t)=\frac{1}{\tau}\exp(-\frac{t}{\tau})$ and $w(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{x^{2}}{2\sigma^{2}}}$. Then $\tilde{\phi}(s)=\frac{1}{1+\tau s}\sim 1-\tau s$ and $\hat{w}(k)=\int_{-\infty}^{+\infty}e^{ikx}w(x)dx=e^{-\frac{1}{2}\sigma^{2}k^{2}}\sim 1-\frac{\sigma^{2}k^{2}}{2}$. Substituting the asymptotic forms of $\tilde{\phi}(s)$ and $\hat{w}(k)$ into Eq. (30) and conducting inverse Fourier-Laplace transform, we have

[TABLE]

(b) Tempered power-law waiting time and Gaussian jump length distributions

Plugging $\tilde{\phi}(s)\sim 1-\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}(s+\lambda)^{\alpha}+\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}\lambda^{\alpha}$ and $\hat{w}(k)\sim 1-\frac{\sigma^{2}k^{2}}{2}$ into Eq. (30) and performing inverse Laplace-Fourier transform give that

[TABLE]

where $H(t;x)$ is defined by its Laplace transform $\tilde{H}(s;x)=\frac{1}{(s+\lambda-r(x)-ipU(x))^{\alpha}-\lambda^{\alpha}}$. It can also be written as

[TABLE]

where the tempered fractional substantial derivative $\mathbb{D}_{t}^{1-\alpha,\lambda,x}$ is defined in the Laplace domain as $\mathcal{L}\{\mathbb{D}_{t}^{1-\alpha,\lambda,x}\rho(x,p,t)\}=(s+\lambda-r(x)-ipU(x))^{1-\alpha}\rho(x,p,s)$. And in the time domain,

[TABLE]

which is equivalent to

[TABLE]

It should be noted that Eqs. (32) and (33) are the generalizations of Eqs. (11) and (5) in Wu2016 respectively to the reaction diffusion cases.

(c) Tempered power-law waiting time and jump length distributions

In this case, we assume that the waiting time and jump length obey different tempered power-law distributions. Let $w(x)=C_{\varepsilon}^{-1}e^{-\gamma|x|}|x|^{-1-\beta}\mathbf{1}_{[\varepsilon,+\infty)}(|x|)$ for $0<\beta<2$. The normalization factor $C_{\varepsilon}$ is defined as $C_{\varepsilon}=\int_{\varepsilon}^{+\infty}e^{-\gamma x}x^{-1-\beta}dx+\int_{-\infty}^{-\varepsilon}e^{-\gamma|x|}|x|^{-1-\beta}dx$ to make sure $\int_{-\infty}^{+\infty}w(x)dx=1$. Then it could be calculated (see Appendix C) that $\hat{w}(k)$ has the asymptotic form $\hat{w}(k)\sim 1-\frac{2\Gamma(1-\beta)}{\beta C_{\varepsilon}}(k^{2}+\gamma^{2})^{\beta/2}+\frac{2\Gamma(1-\beta)}{\beta C_{\varepsilon}}\gamma^{\beta}$ for $k\rightarrow 0$. Substituting the asymptotic $\tilde{\phi}(s)$ and $\hat{w}(k)$ into Eq. (30) and conducting inverse Laplace-Fourier transform, we obtain

[TABLE]

where the tempered fractional Riesz derivative $\nabla_{x}^{\beta,\gamma}$ (see Wu2016 for more details) is defined in the Fourier domain as $\mathcal{F}\{\nabla_{x}^{\beta,\gamma}\rho(x,p,t)\}=-(k^{2}+\gamma^{2})^{\beta/2}\rho(k,p,t)$. Equivalently,

[TABLE]

where the operator $\mathbb{D}_{t}^{1-\alpha,\lambda,x}$ is defined as Eq. (34) or Eq. (35). When $r(x)\equiv 0$, Eq. (37) reduces to Eq. (36) in Wu2016 .

III Derivation of the Backward Feynman-Kac Equations

Let $\rho_{x_{0}}(A,t)$ be the PDF of the functional $A$ at time $t$ with the initial position $x_{0}$. The backward Feynman-Kac equation regarding $\rho_{x_{0}}(A,t)$ with the nonlinear reaction rate $r(\rho(x,t))$ can be written as follows

[TABLE]

It should be noted that $\rho(x_{0},t)\neq\int_{-\infty}^{+\infty}\rho_{x_{0}}(A,t)dA$. While the first term on the RHS of Eq. (38) indicates the motionless particles remaining at their initial position $x_{0}$ up to time $t$, the second term alternatively represents the particles that jump to the location $x_{0}+z$ at time $\tau$ ($\tau<t$). Similarly, conducting the Fourier transform $A\rightarrow p$, we obtain

[TABLE]

In what follows, we consider a special case of Eq. (39), which supposes the reaction rate is $r(x)$ and the functional is defined as $A=\int_{0}^{t}U[x(t^{\prime})]dt^{\prime}$. Thus, under these assumptions, Eq. (39) can be rewritten as

[TABLE]

Consequently, the standard Laplace-Fourier transform technique is applicable to further simplify the concerned equation. Conducting to Eq. (40) the Laplace transform $t\rightarrow s$ and the Fourier transform $x_{0}\rightarrow k$, we get

[TABLE]

Here we list three kinds of typical waiting time and jump length distributions and write down the respective backward equations in the continuum limit. We omit the power-law waiting time distribution cases, since they could be recovered when the exponential tempering exponents are set to zeros.

(a) Exponential waiting time and Gaussian jump length distributions

Plugging $\tilde{\phi}(s)\sim 1-\tau s$ and $\hat{w}(k)\sim 1-\frac{\sigma^{2}k^{2}}{2}$ into Eq. (41) and performing inverse Laplace-Fourier transform give that

[TABLE]

(b) Tempered power-law waiting time and Gaussian jump length distributions

With the substitution of $\tilde{\phi}(s)\sim 1-\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}(s+\lambda)^{\alpha}+\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}\lambda^{\alpha}$ and $\hat{w}(k)\sim 1-\frac{\sigma^{2}k^{2}}{2}$ into Eq. (41) and inverse Laplace-Fourier transform, we obtain

[TABLE]

where $H(t;x_{0})$ is defined by its Laplace transform $\tilde{H}(s;x_{0})=\frac{1}{(s+\lambda-r(x_{0})-ipU(x_{0}))^{\alpha}-\lambda^{\alpha}}$. Equivalently,

[TABLE]

where the tempered fractional substantial derivative $\mathbb{D}_{t}^{1-\alpha,\lambda,x_{0}}$ is defined in the Laplace domain as $\mathcal{L}\{\mathbb{D}_{t}^{1-\alpha,\lambda,x_{0}}\rho_{x_{0}}(p,t)\}=(s+\lambda-r(x_{0})-ipU(x_{0}))^{1-\alpha}\rho_{x_{0}}(p,s)$. And in the time domain,

[TABLE]

which is equivalent to

[TABLE]

(c) Tempered power-law waiting time and jump length distributions

Substituting $\tilde{\phi}(s)\sim 1-\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}(s+\lambda)^{\alpha}+\frac{\Gamma(1-\alpha)}{\alpha C_{\tau}}\lambda^{\alpha}$ and $\hat{w}(k)\sim 1-\frac{2\Gamma(1-\beta)}{\beta C_{\varepsilon}}(k^{2}+\gamma^{2})^{\beta/2}+\frac{2\Gamma(1-\beta)}{\beta C_{\varepsilon}}\gamma^{\beta}$ into Eq. (41) and doing Laplace-Fourier transform lead to

[TABLE]

where the operator $\nabla_{x_{0}}^{\beta,\gamma}$ is defined in the Fourier domain as $\mathcal{F}\{\nabla_{x_{0}}^{\beta,\gamma}\rho_{x_{0}}(p,t)\}=-(k^{2}+\gamma^{2})^{\beta/2}\rho_{k}(p,t)$. Equivalently,

[TABLE]

where the operator $\mathbb{D}_{t}^{1-\alpha,\lambda,x_{0}}$ is defined as Eq. (45) or Eq. (46).

IV Applications of the derived equations

In this section, we present the distributions of specific functionals of the paths of particles performing temporal tempered anomalous dynamics with piecewise constant reaction rate. We analytically solve the corresponding backward Feynman-Kac equations for obtaining the distributions, the moments, and other properties of interest. As pointed out in Wu2016 , our analysis is mainly based on the derived backward Feynman-Kac equations, since here we are more concerned with the functional distributions than the particles’ positions. From this perspective, backward equations are more convenient, although the forward equations could lead to the same conclusions with the extra integration over $x$ step. In what follows, we assume the reaction rate function $r(x)$ satisfying

[TABLE]

where $\kappa_{1}$ and $\kappa_{2}$ are negative constants. Here we confine our analysis to negative reaction rates in order to model the spontaneous evanescent process in which the particles are destroyed or removed at different constant rates depending on their positions. As for the reproduction process with positive reaction rates, one must specify the rules regarding the waiting time of the newborn particles Abad2010 ; Fedotov2010 , which is a problem to be explored further.

IV.1 Occupation time in half-space and its fluctuations

Define the occupation time of a particle in the positive half-space as $T^{+}=\int_{0}^{t}U(x(\tau))d\tau$, where $U(x)=1$ for $x\geq 0$ and $U(x)=0$ for $x<0$. Since obviously $T^{+}\geq 0$, we rely on the Laplace transform $T^{+}\rightarrow p$, instead of the Fourier transform. In order to find the PDF of $T^{+}$, we consider the backward equation Eq. (43) (or Eq. (44)), conduct Laplace transform $t\rightarrow s$, and substitute the assumed $U(x_{0})$ and $r(x_{0})$. Denote $K_{\alpha}=\frac{\alpha C_{\tau}\sigma^{2}}{2\Gamma(1-\alpha)}$. Thus,

[TABLE]

Solving the equations in each half-space individually and demanding $\rho_{x_{0}}(p,s)<\infty$ for $|x_{0}|\rightarrow\infty$, we have

[TABLE]

To determine the constants $C_{1}$ and $C_{2}$, we require that $\rho_{x_{0}}(p,s)$ and its first derivative $\frac{\partial\rho_{x_{0}}(p,s)}{\partial x_{0}}$ are continuous at $x_{0}=0$. Consequently, we have

[TABLE]

Suppose the particle starts from $x_{0}=0$. Then

[TABLE]

which describes the PDF of $T^{+}$ in the Laplace domain ($T^{+}\rightarrow p$ and $t\rightarrow s$). When both $\kappa_{1}$ and $\kappa_{2}$ are set to zeros in Eq. (53), the Eq. (47) in Wu2016 is recovered, describing the distribution of $T^{+}$ for tempered anomalous motions without reactions. However, it seems difficult to invert Eq. (53) analytically, even for the $\kappa_{1}=\kappa_{2}=0$ case as reported in Wu2016 , while the occupation fraction $T^{+}/t$ obeys Lamperti distribution for anomalous motions without exponential tempering ($\lambda=0$) and reactions ($\kappa_{1}=\kappa_{2}=0$) Carmi2010 . Specially, when $\alpha=1$, $\rho_{0}(p,s)=(s-\kappa_{2})^{-1/2}(s+p-\kappa_{1})^{-1/2}$, apparently different from the arcsine law of the occupation fraction for Brownian motion Majumdar2005 ; Barkai2006 .

Furthermore, in order to evaluate the expectation and fluctuation of the occupation time, we calculate the moments of $T^{+}$ in the Laplace domain ($t\rightarrow s$) from Eq. (53) as follows:

[TABLE]

According to Eq. (54), when $\kappa_{1}=\kappa_{2}$, $\langle T^{+}\rangle_{s}=\frac{1}{2(s-\kappa_{1})^{2}}$, and correspondingly in the time domain $\langle T^{+}\rangle=\frac{te^{\kappa_{1}t}}{2}$, or $\langle\frac{T^{+}}{t}\rangle=\frac{e^{\kappa_{1}t}}{2}$. Here $e^{\kappa_{1}t}$ can be interpreted as the survival probability. This result is further confirmed by the simulation results in Figs. 1-2 and it is remarkably different from the previous studies, for example Wu2016 ; CarmiBarkai2011 , which presents $\langle T^{+}\rangle=\frac{t}{2}$, or $\langle\frac{T^{+}}{t}\rangle=\frac{1}{2}$ for (tempered) anomalous diffusion processes.

[TABLE]

where $B_{1}(s)=2(\alpha-1)\lambda^{\alpha}+(2+\alpha)(s+\lambda-\kappa_{1})^{\alpha}$ and $B_{2}(s)=2(\alpha-1)\lambda^{\alpha}+(2-\alpha)(s+\lambda-\kappa_{1})^{\alpha}$. Substitution of $\kappa_{1}=\kappa_{2}=0$ into Eq. (55) gives the following special case of our derivation

[TABLE]

which is exactly in agreement with the previous work (Eq. (54) in Wu2016 ).

Particularly, assume $\kappa_{1}=\kappa_{2}$, we obtain from Eq. (55) that

[TABLE]

and inversely in the time domain,

[TABLE]

As $t\rightarrow 0$, namely $s\rightarrow\infty$, both $\lambda$ and $\kappa_{1}$ can be ignored in Eq. (57). Thus, $\langle(T^{+})^{2}\rangle\sim\frac{4-\alpha}{8}t^{2}$, which is the expected result Wu2016 since initially tempering and reaction terms have negligible influence on the process.

IV.2 First passage time

The time $T_{f}$ when a particle starting at $x_{0}=-b$ $(b>0)$ passes $x=0$ for the first time is called the first passage time Redner2001 . According to Kac1951 ; Wu2016 , the relationship between the distribution of $T_{f}$ and the PDF of the occupation time in half-space $T^{+}$ (in the Laplace domain $T^{+}\rightarrow p$) satisfies

[TABLE]

Denote the PDF of $T_{f}$ as $f(t)$, which satisfies $f(t)=-\frac{\partial}{\partial t}\lim_{p\rightarrow\infty}\rho_{-b}(p,t)$ from Eq. (59). According to Eqs. (51) and (52),

[TABLE]

Hence, in the Laplace domain $t\rightarrow s$, we have

[TABLE]

When $t\rightarrow 0$, corresponding to $s\rightarrow\infty$, we have $\tilde{f}(s)\sim\exp\left(-\frac{b}{\sqrt{K_{\alpha}}}s^{\frac{\alpha}{2}}\right)$ from Eq. (61), which is equivalent to the one-sided Lévy law in the time domain. Furthermore, we obtain $\tilde{f}(s)\sim\frac{-\kappa_{2}}{s-\kappa_{2}}$ as $s\rightarrow 0$ and equivalently $f(t)\sim-\kappa_{2}e^{\kappa_{2}t}$ as $t\rightarrow\infty$. Consequently, we assert that $f(t)$ decays very fast to zero when $t\rightarrow 0$, behaves as $t^{-1-\alpha/2}$ for short but not too short time scales and finally decays to zero according to the exponential law $-\kappa_{2}e^{\kappa_{2}t}$ for adequately large time scales. This conclusion is further confirmed by simulations presented in Figs. 3 and 4. Note that the exponential law is distinct from the celebrated $t^{-1.5}$ decay law Wu2016 ; Redner2001 , due to the influence of additional reaction terms.

IV.3 Occupation time in half-interval

In this subsection, we suppose that the particle motion is restricted to the interval $(-L,L)$ with the absorbing or reflecting boundaries. The occupation time $T^{+}_{a}$ or $T^{+}_{r}$ (corresponding to the absorbing or reflecting boundary conditions, respectively) in the positive half-interval is defined as $T^{+}_{a}=T^{+}_{r}=\int_{0}^{t}U[x(\tau)]d\tau$, where $U(x)=1$ for $0\leq x<L$ and $U(x)=0$ for $-L<x<0$. Thus, we have

[TABLE]

(a) Absorbing boundary conditions

Solving Eq. (62) in each interval, respectively, we obtain

[TABLE]

With the absorbing boundary condition $\rho_{x_{0}}(p,s)|_{x_{0}=\pm L=0}$ and the continuity of $\rho_{x_{0}}(p,s)$ and its first derivative at $x_{0}=0$, we have

[TABLE]

where we denote for simplicity $E_{1}(p,s)=\sqrt{\frac{(s+\lambda+p-\kappa_{1})^{\alpha}-\lambda^{\alpha}}{K_{\alpha}}}$, $E_{2}(s)=\sqrt{\frac{(s+\lambda-\kappa_{2})^{\alpha}-\lambda^{\alpha}}{K_{\alpha}}}$ and $E(p,s)=\frac{E_{2}(s)}{E_{1}(p,s)}$. Suppose the particle departs at $x=0$. Then

[TABLE]

(b) Reflecting boundary conditions

Solving Eq. (62) with $\frac{\partial\rho_{x_{0}}(p,s)}{\partial x_{0}}|_{x_{0}=\pm L}=0$ gives

[TABLE]

Similarly, in order to determine the constants $C_{1}$ and $C_{2}$, we assume the continuity of $\rho_{x_{0}}(p,s)$ and its first derivative at $x_{0}=0$ in Eq. (66). Hence,

[TABLE]

where we denote $F(p,s)=\frac{E_{2}(s)\sinh[LE_{2}(s)]}{E_{1}(p,s)\sinh[LE_{1}(p,s)]}$ for simplicity. If the particle starts from $x_{0}=0$, then

[TABLE]

Especially, when $\lambda=\kappa_{1}=\kappa_{2}=0$, we recover

[TABLE]

which was previously derived in Carmi2011 using the similar method.

V Summary

The functional distributions of the trajectories of diffusion processes have been well developed, including the ones for normal and anomalous diffusions. Chemical reaction is another important process that results in the transformation of chemical substances. What happens for the functional distribution if we have both reaction and diffusion? This paper is answering this question. The previous ideas of deriving the governing equations for the functional distributions of pure (normal/anomalous) diffusion processes do not work for the general reaction diffusion processes. We provide a theoretical framework of deriving the governing equations for the functional distributions of the trajectories of the stochastic processes with both reaction and diffusion.

It is well known that the net effect of normal diffusion and chemical reaction processes is just the sum of the individual rates of change. However, the net effect becomes complicated for the two processes: anomalous diffusion and chemical reaction. What about the functional distributions of the net effect of the diffusion and reaction processes? We show that for any type of diffusion, if the reactional rate is a constant or only time dependent $r(t)$, then the functional distribution of the diffusion and reaction processes is the one of pure diffusion process multiplied by $e^{\int_{0}^{t}r(u)du}$. For the case with general reaction rates, we derive a series of specific Feynman-Kac equations with various diffusion types.

Several functional distributions of interest are investigated by means of analytically solving the derived backward Feynman-Kac equations, including the occupation time in half-space, the first-passage time to a fixed boundary and the occupation time in half-interval with absorbing or reflecting boundary conditions. Compared with the previous work, we are more concerned with the influence of the reaction terms to the temporal-tempered anomalous dynamics and validate some new theoretical conclusions by simulations. It is first found that the mean occupation time in half-space of the reaction diffusion process obeys the $\frac{te^{\kappa_{1}t}}{2}$ law, where $\kappa_{1}$ is the assumed reaction rate in the whole space. The PDF of the first passage time to the fixed boundary $x=0$ starting from $x=-b$ has the asymptotic form $-\kappa_{2}e^{\kappa_{2}t}$ as $t\rightarrow\infty$, and $\kappa_{2}$ is the prescribed reaction rate in the negative half-space.

The coming research topics that need to explore should include functional distributions of particle creation processes, the first-passage time to a moving boundary, reaction diffusion dynamics under the influence of external potentials and the corresponding functional distributions.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant No.11671182. We thank Eli Barkai for the helpful discussions.

APPENDIX

V.1 Taylor expansion of the first term on the RHS of Eq. (15)

We express the second order Taylor expansion of the first term on the RHS of Eq. (15) for small $z$ as follows:

[TABLE]

since we have the assumptions $\int_{\mathbb{R}}w(z)dz=1$, $\int_{\mathbb{R}}zw(z)dz=0$ and $\int_{\mathbb{R}}z^{2}w(z)dz=\sigma^{2}$.

V.2 Asymptotic forms of the Laplace transforms of (tempered) power-law waiting time PDF

In the first place, consider the power-law waiting time PDF $\phi(t)=\alpha\tau^{\alpha}t^{-1-\alpha}\mathbf{1}_{[\tau,+\infty)}(t)$ with $0<\alpha<1$, which obviously satisfies the normalization $\int_{0}^{+\infty}\phi(t)dt=1$. According to the formula $\int_{0}^{+\infty}(e^{iky}-1)\alpha y^{-1-\alpha}dy=-\Gamma(1-\alpha)(-ik)^{\alpha}$ for $0<\alpha<1$ in Meerschaert2012 , we have

[TABLE]

Usually, $O(s)$ is omitted in practice since we are mainly attentive to the long-time asymptotic behaviour of $\phi(t)$, which corresponds to $\tilde{\phi}(s)$ in the limit of $s\rightarrow 0$.

Next, we proceed to calculate the Laplace transform of the tempered power-law PDF $\phi(t)=C_{\tau}^{-1}e^{-\lambda t}t^{-1-\alpha}\mathbf{1}_{[\tau,+\infty)}(t)$ ($0<\alpha<1$ and $C_{\tau}$ is the normalization factor) as follows:

[TABLE]

V.3 Asymptotic forms of the Fourier transform of tempered power-law jump length PDF

As for the tempered power-law distribution $w(x)=C_{\varepsilon}^{-1}e^{-\gamma|x|}|x|^{-1-\beta}\mathbf{1}_{[\varepsilon,+\infty)}(|x|)$ with $0<\beta<2$ and the normalization factor $C_{\varepsilon}$, we calculate its Fourier transform $\hat{w}(k)=\int_{-\infty}^{+\infty}e^{ikx}w(x)dx$ from two aspects. First of all, when $0<\beta<1$, we have

[TABLE]

where $\theta=\arg(\gamma+ik)=\arctan(k/\gamma)$. Thus, as $k\rightarrow 0$, $\cos\theta\beta\rightarrow 1$. Secondly, when $1<\beta<2$, we have

[TABLE]

according to the recurrence relation $\Gamma(1-\beta)=\frac{\Gamma(2-\beta)}{1-\beta}$.