Fractional compound Poisson processes with multiple internal states

TL;DR

This paper develops mathematical models for particles exhibiting anomalous diffusion with multiple internal states, deriving equations to describe their positions, internal state distributions, and related functionals, with applications in first passage and occupation times.

Contribution

It introduces new fractional Fokker-Planck and Feynman-Kac equations for multi-state anomalous diffusion processes, including internal state dynamics.

Findings

Derived equations for particle positions and internal states

Analyzed stochastic process dynamics

Applied models to first passage and occupation time distributions

Abstract

For the particles undergoing the anomalous diffusion with different waiting time distributions for different internal states, we derive the Fokker-Planck and Feymann-Kac equations, respectively, describing positions of the particles and functional distributions of the trajectories of particles; in particular, the equations governing the functional distribution of internal states are also obtained. The dynamics of the stochastic processes are analyzed and the applications, calculating the distribution of the first passage time and the distribution of the fraction of the occupation time, of the equations are given.

Equations

$$\big{|}G(k,s)\big{>}=\frac{1}{s}\big{[}I-\Phi(s)\big{]}\big{[}I-M^{T}\Phi(s)\Lambda(k)\big{]}^{-1}\big{|}{\rm init}\big{>}.$$

$$\begin{split}M^{T}\frac{\partial}{\partial t}\big{|}G(x,t)\big{>}=&\big{(}M^{T}-I\big{)}{\rm diag}\big{(}B_{\alpha_{1}}^{-1},\ldots,B_{\alpha_{N}}^{-1}\big{)}D^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}_{t}\big{|}G(x,t)\big{>}\\ &+M^{T}{\rm diag}\big{(}K_{\alpha_{1}},\ldots,K_{\alpha_{N}}\big{)}D^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}_{t}\frac{\partial^{2}}{\partial x^{2}}\big{|}G(x,t)\big{>},\end{split}$$

$$A^{-1}(s)\thicksim\frac{\big{|}{\rm eq_{M}}\big{>}\big{<}\Sigma\big{|}}{\big{<}\Sigma\big{|}\Psi(s)\big{|}{\rm eq_{M}}\big{>}+\sigma^{2}k^{2}\big{<}\Sigma\big{|}{\rm eq_{M}}\big{>}}.$$

$$g(k,s)=\big{<}\Sigma\big{|}G(k,s)\big{>}\thicksim\frac{1}{s}\frac{\big{<}\Sigma\big{|}\Psi(s)\big{|}{\rm eq_{M}}\big{>}}{\big{<}\Sigma\big{|}\Psi(s)\big{|}{\rm eq_{M}}\big{>}+\sigma^{2}k^{2}}.$$

$$\big{<}x^{2}(t)\big{>}\sim\mathcal{L}^{-1}\Bigg{\{}\frac{2\sigma^{2}}{s(\varepsilon_{1}B_{\alpha_{1}}s^{\alpha_{1}}+\ldots+\varepsilon_{N}B_{\alpha_{N}}s^{\alpha_{N}})}\Bigg{\}}.$$

$$g(k,s)\sim\sum_{\mathop{i=1}\limits_{i\neq i_{1},i_{2},\cdots,i_{j}}}^{j}\frac{1}{s}\frac{\big{<}\Sigma_{i}\big{|}\Psi_{i}(s)\big{|}{\rm eq_{M}}\big{>}_{i}\big{<}\Sigma_{i}\big{|}{\rm init}\big{>}_{i}}{\big{<}\Sigma_{i}\big{|}\Psi_{i}(s)\big{|}{\rm eq_{M}}\big{>}_{i}+\sigma^{2}k^{2}\big{<}\Sigma_{i}\big{|}{\rm eq_{M}}\big{>}_{i}},$$

$$g(k,s)\sim\sum_{i=1}^{N}\frac{1}{s}\frac{\lambda_{i}s^{\alpha_{i}}}{s^{\alpha_{i}}+K_{\alpha_{i}}k^{2}}\sim\frac{1}{s}-k^{2}\Bigg{[}\sum_{i=1}^{N}K_{\alpha_{i}}\lambda_{i}s^{-\alpha_{i}-1}\Bigg{]},$$

$$\big{<}x^{2}(t)\big{>}\sim\mathcal{L}^{-1}\Bigg{\{}\sum_{\mathop{i=1}\limits_{i\neq i_{1},i_{2},\cdots,i_{j}}}^{j}\frac{2\sigma^{2}\big{<}\Sigma_{i}\big{|}{\rm init}\big{>}_{i}\big{<}\Sigma_{i}\big{|}{\rm eq_{M}}\big{>}_{i}}{s\big{<}\Sigma_{i}\big{|}\Psi_{i}(s)\big{|}{\rm eq_{M}}\big{>}_{i}}\Bigg{\}},$$

$$\begin{split}M^{T}\frac{\partial}{\partial t}\big{|}G(x,\rho,t)\big{>}=&(M^{T}-I){\rm diag}\big{(}B_{\alpha_{1}}^{-1},\ldots,B_{\alpha_{N}}^{-1}\big{)}\mathcal{D}_{t}^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}\big{|}G(x,\rho,t)\big{>}\\ &+M^{T}\frac{\partial^{2}}{\partial x^{2}}{\rm diag}(K_{\alpha_{1}},\ldots,K_{\alpha_{N}})\mathcal{D}_{t}^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}\big{|}G(x,\rho,t)\big{>}-\rho U(x)M^{T}\big{|}G(x,\rho,t)\big{>},\end{split}$$

$$\begin{split}M\frac{\partial}{\partial t}\big{|}G_{x_{0}}(\rho,t)\big{>}=&{\rm diag}\big{(}B_{\alpha_{1}}^{-1},\ldots,B_{\alpha_{N}}^{-1}\big{)}\mathcal{D}_{t}^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}(M-I)\big{|}G_{x_{0}}(\rho,t)\big{>}\\ &+{\rm diag}(K_{\alpha_{1}},\ldots,K_{\alpha_{N}})\mathcal{D}_{t}^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}M\frac{\partial^{2}}{\partial x_{0}^{2}}\big{|}G_{x_{0}}(\rho,t)\big{>}-\rho U(x_{0})M\big{|}G_{x_{0}}(\rho,t)\big{>}.\end{split}$$

$$\lim_{\rho\rightarrow\infty}g_{0}(\rho,s)=\frac{1}{s}\Bigg{[}1-{\rm exp}\Big{(}-\sqrt{\frac{a_{0}+b_{0}}{2}}B\Big{)}\Bigg{]},$$

$$f(t)=\mathcal{L}^{-1}\Bigg{\{}\exp\Big{(}-\sqrt{\frac{a_{0}+b_{0}}{2}}B\Big{)}\Bigg{\}}.$$

$$M^{T}\frac{\partial}{\partial t}\big{|}G(\rho_{s},t)\big{>}=(M^{T}-I){\rm diag}\Big{\{}\frac{1}{B_{\alpha_{1}}}\mathfrak{D}_{t}^{1-\alpha_{1}},\ldots,\frac{1}{B_{\alpha_{N}}}\mathfrak{D}_{t}^{1-\alpha_{N}}\Big{\}}\big{|}G(\rho_{s},t)\big{>}-\rho_{s}M^{T}{\rm diag}\{U(1),\ldots,U(N)\},$$

$$\begin{split}\mathfrak{D}_{t}^{1-\alpha_{i}}G(\rho_{s},t)=\frac{1}{\Gamma(\alpha_{i})}\Big{[}\frac{\partial}{\partial t}+\rho_{s}U(i)\Big{]}\int_{0}^{t}\frac{\exp(-(t-\tau)\rho_{s}U(i))}{(t-\tau)^{1-\alpha_{i}}}G(\rho_{s},\tau)d\tau.\end{split}$$

$$\underset{t\to \infty }{\lim }{l}_{t^{(1)}/t}(x)=\frac{\sin (\pi \alpha )}{\pi }\frac{{\epsilon }_1({\epsilon }_2+\dots +{\epsilon }_m)(1-x{)}^{\alpha -1}{x}^{\alpha -1}}{({\epsilon }_2+\dots +{\epsilon }_m{)}^2{x}^{2\alpha }+{\epsilon }_1^2(1-x{)}^{2\alpha }+2{\epsilon }_1({\epsilon }_2+\dots +{\epsilon }_m)\cos (\alpha \pi )x^{\alpha }(1-x{)}^{\alpha }},$$

Full text

††preprint: APS/123-QED

Introduction.— Poisson process is the most fundamental stochastic process of renewal theory. The application of Poisson process naturally coming to our minds is in queueing theory to model the random events: the arrival of customers at a store or phone calls at an exchange Kleinrock:76 . Renewal process generalizes Poisson process for arbitrary holding times, being still independent identically distributed (i.i.d.) godreche:00 . If the probability density function (PDF) of the holding/waiting times between two subsequent events has the asymptotic behavior $\phi(\tau)\sim 1/\tau^{\alpha+1},0<\alpha<1$ when time is long enough laskin:00 , it is called fractional Poisson process. The residence time statistics for $N$ fractional Poisson processes are discussed in burov:00 , which reveal sharp transitions for a critical number of degrees of freedom $N$ and give application of detecting nonergodic kinetics from the measurements of many blinking chromophores.

We further generalize the renewal processes to have multiple internal states, where the holding times for different internal states are drawn from different distributions. The case of two internal states is considered in lowen:00 ; Godec:17 with applications, including trapping in amorphous semiconductors, electronic burst noise, movement in systems with fractal boundaries, the digital generation of $1/f$ noise, and ionic currents in cell membranes. The fractional Poisson processes with multiple internal states have a lot of potential applications, e.g., the particles moving in multiphase viscous liquid composed of materials with different chemical properties; Niemann, Barkai, and Kantz niemann:00 detailedly investigate a stochastic signal with multiple states, in which each state has an associated joint distribution for the signal’s intensity and its holding time. This letter adopts the notations of niemann:00 . If introducing jumps to the fractional Poisson processes, then we reach the fractional compound Poisson (fcP) processes meerschaert:00 . The jumps follow the fractional Poisson processes to arrive and, most of the time, the size of the jumps is random, with specified PDF. More specifically, let $\mathcal{N}(t)$ be the fractional Poisson process and $\{\xi_{i},i=1,2,\ldots\}$ a sequence of i.i.d. random variables (jump lengths). Then $X(t)=\sum\limits_{i=1}^{\mathcal{N}(t)}\xi_{i}$ is the fcP process. The continuous time random walk (CTRW) with i.i.d. power law waiting time cartea:00 ; metzler:00 ; Metzler:16 ; klafter:00 is a specific fcP process, being widely used to model various anomalous diffusions, e.g., mRNA molecules in living cells golding:00 , price fluctuation in financial market scalas:00 .

In this letter, we focus on the fcP processes $X(t)$ with multiple internal states, i.e., $\mathcal{N}(t)$ of $X(t)$ is a fractional Poisson process with multiple internal states. Each internal state has an associated distribution of waiting time, but the distributions of jump lengths are all simply taken as normal distribution. We derive the Fokker-Planck equations metzler:00 , governing the PDF of positions of $X(t)$, the Feynman-Kac equations turgeman:00 ; carmi:00 ; cairoli:00 , describing the distribution of the functional satya:00 of the paths of $X(t)$, and the equations, characterizing the functional distribution of the internal states. From the Fokker-Planck equations, we obtain the evolution of the mean square displacement (MSD) for the process, showing that it strongly depends on the properties of the internal transition matrix feller:00 and sometimes the distribution of the initial position of the particles. The applications of the Feynman-Kac equations and the equations governing the distribution of the functional are given to calculate the distribution of first passage time carmi:00 ; redner:00 ; Deng:17 and the distribution of the fraction of the occupation time godreche:00 . Some properties of the distributions are obtained.

Model.— We consider the fcP processes with finite internal states, denoting their number as $N$. The internal states determine the distributions of waiting times and the transition of the internal states is described by a Markov chain with its transition matrix $M$; the dimension of $M$ is $N\times N$. The element $m_{ij}$ of the matrix $M$ represents the transition probability from state $i$ to state $j$. The bras \big{<}\cdots\big{|} and kets \big{|}\cdots\big{>} denote the row and column vectors, respectively. From Ref. feller:00 one can see that for the ergodic or periodic chain there always exists an equilibrium distribution, denoted by \big{<}\rm eq_{M}\big{|}, with the property of \big{<}{\rm eq_{M}}\big{|}M=\big{<}{\rm eq_{M}}\big{|}. For the transition matrix $M$, as is well known, its largest eigenvalue is $1$ and the right eigenvector for this eigenvalue is the vector of 1s Walker:11 , denoted as \big{|}\Sigma\big{>}. That is M\big{|}\Sigma\big{>}=\big{|}\Sigma\big{>}. The second largest eigenvalue is strictly less than $1$ if $M$ is irreducible. It also obviously holds that M^{T}\big{|}{\rm eq_{M}}\big{>}=\big{|}{\rm eq_{M}}\big{>} and \big{<}\Sigma\big{|}M^{T}=\big{<}\Sigma\big{|}. We use the notation \big{<}{\rm init}\big{|} to represent the initial distribution of the internal states. Based on the CTRW, we define the waiting time distribution matrix $\Phi(t)={\rm diag}(\phi^{(1)}(t),\phi^{(2)}(t),\ldots,\phi^{(N)}(t))$ and the jump length one $\Lambda(x)={\rm diag}(\lambda^{(1)}(x),\lambda^{(2)}(x),\ldots,\lambda^{(N)}(x))$.

The initial state of the stochastic process $X(t)$ is drawn from the initial distribution. After confirming the initial state, e.g., state $i$, the waiting time and jump length will be obtained from the distributions $\phi^{(i)}(t)$ and $\lambda^{(i)}(x)$, respectively. The new internal state is drawn from the distribution M\big{|}{\rm init}\big{>}. Then repeat the procedure. For the Markov chain, there are plenty of practically or theoretically important transition matrices feller:00 . The simple and representative one should be $\begin{pmatrix}0&1\\ 1&0\\ \end{pmatrix},$ signifying $2$ alternating internal states. The model will be, respectively, discussed for the case that $M$ is irreducible and the case that $M$ is reducible, i.e., the digraph of $M$ is strongly connected or not.

Fokker-Planck equations.—We use the notation $g^{(i)}(x,t)$, $i=1,2,\ldots,N$ to represent the PDF of finding the particle, at time $t$, position $x$ and internal state $i$. Let \big{|}G(x,t)\big{>} be the column vector $\{g^{(i)}(x,t),\,i=1,\ldots,N\}$. The initial condition for \big{|}G(x,t)\big{>} is taken as \delta(x)\delta(t)\big{|}{\rm init}\big{>}. Similarly to the derivation of fractional Fokker-Planck equation from CTRW model with i.i.d. waiting times klafter:00 ; metzler:00 , we can obtain the equation in the Fourier-Laplace space

[TABLE]

In this letter, we take the waiting time distributions as asymptotical power laws, i.e., in the Laplace space $\Phi(s)=I-\Psi(s)$, where $\Psi(s)={\rm diag}(B_{\alpha_{1}}s^{\alpha_{1}},\ldots,B_{\alpha_{N}}s^{\alpha_{N}})$, $0<\alpha_{1},\ldots,\alpha_{N}<1$. As for the jump length distributions we choose all of them as Gaussian distribution, so in the Fourier space $\Lambda(k)=(1-\sigma^{2}k^{2})I$. Thus from the above equations, taking the inverse Laplace and Fourier transformations leads to the Fokker-Planck equation with $N$ internal states

[TABLE]

where D^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}_{t}={\rm diag}\Big{(}D_{t}^{1-\alpha_{1}},\ldots,D_{t}^{1-\alpha_{N}}\Big{)}, and $D_{t}^{1-\alpha_{i}},i=1,\ldots,N$ are the Riemann-Liouville derivatives; the factors $K_{\alpha_{i}}=\sigma^{2}/B_{\alpha_{i}}$ of the diagonal matrix represent diffusion coefficients with the dimension ${\rm cm}^{2}/{\rm sec}^{\alpha}$; it can be noted that if $N=1$ the usual fractional Fokker-Planck equation metzler:00 is recovered.

Next our aim is to calculate the PDF of finding the particle at position $x$ at time $t$, denoted as $g(x,t)$, and the mean squared displacement (MSD) of the process. Without loss of generality, we assume $0<\alpha_{1}\leqslant\alpha_{2}\leqslant\ldots\leqslant\alpha_{N}<1$, and take the equilibrium and initial distributions, respectively, as \big{|}\rm eq_{M}\big{>}=(\varepsilon_{1},\ldots,\varepsilon_{N})^{T} and \big{<}{\rm init}\big{|}=(\lambda_{1},\ldots,\lambda_{N}). We denote the matrix $I-M^{T}\Phi(s)\Lambda(k)$ in Eq. (1) as $A(s)$; it is irreducible or not if and only if $M$ is or not. Plugging the distributions of waiting time and jump length into $A(s)$, for the irreducible transition matrix $M$ the asymptotic expression of the inverse matrix of $A(s)$ can be expressed as

[TABLE]

Then we obtain

[TABLE]

The PDF in the Fourier-Laplace space given in Eq. (3) is the same as the PDF of (natural-form) distributed-order diffusion discussed in sandev:00 ; chechkin:00 with $p(\alpha)=\sum^{N}_{i=1}\varepsilon_{i}\delta(\alpha-\alpha_{i})$, even though their backgrounds are completely different. From Eq. (3), the MSD can be got as

[TABLE]

Therefore, when $t$ is large enough ($s$ tends to [math]), the MSD of the process behaves as \big{<}x^{2}(t)\big{>}\sim\mathcal{L}^{-1}\Big{\{}\frac{2K_{\alpha_{1}}}{\varepsilon_{1}s^{1+\alpha_{1}}}\Big{\}}\sim\frac{2K_{\alpha_{1}}}{\varepsilon_{1}\Gamma(1+\alpha_{1})}t^{\alpha_{1}}; i.e., the MSD of the process with the irreducible transition matrix behaves asymptotically as $t^{\alpha}$, where $\alpha$ is the smallest one among all exponents of the power-law waiting time distributions of the internal states. When the transition matrix of the Markov chain of the internal states is irreducible, the equilibrium distribution \big{|}{\rm eq_{M}}\big{>} does not depend on the initial distribution, and the initial distribution has no influence on the final PDF or MSD.

On the other hand, if the transition matrix is reducible, then the initial distribution of the internal states often makes an important impact on the final results. First, we consider the case that the transition matrix has the form of {\rm diag}\big{\{}M_{1}(s),\ldots M_{j}(s)\big{\}}, where the matrices $M_{i},i=1,\ldots,j$ are irreducible with the dimension of $n_{i}\times n_{i}$, and $n_{1}+n_{2}+\ldots+n_{j}=N$. From the form of the transition matrix, one can see that the internal states of the process actually consist of several independent Markov chains with the transition matrices $M_{1},M_{2},\ldots,M_{j}$. Following the structure of $M$, since $\Psi(s)$ and $I$ are diagonal, we can rewrite $\Psi(s)$ and $I$ as the form consistng of $\Psi_{i}(s)$ and $I_{i}$, $i=1,2,\ldots,j$. Then the matrix $A(s)$ has the form of {\rm diag}\big{\{}A_{1}(s),\ldots,A_{j}(s)\big{\}}, where $A_{i}(s)=I_{i}-M_{i}^{T}\Psi_{i}(s)+\sigma^{2}k^{2}M_{i}^{T}$. The vectors can also be rewritten as \big{|}{\rm init}\big{>}=\big{(}\big{|}{\rm init}\big{>}_{1},\ldots,\big{|}{\rm init}\big{>}_{j}\big{)}^{T}, \big{|}{\rm eq_{M}}\big{>}=\big{(}\big{|}{\rm eq_{M}}\big{>}_{1},\ldots,\big{|}{\rm eq_{M}}\big{>}_{j}\big{)}^{T}, and \big{<}\Sigma\big{|}=\big{(}\big{<}\Sigma\big{|}_{1},\ldots,\big{<}\Sigma\big{|}_{j}\big{)}. For the convenience of statement, we redefine the subscripts, i.e., let $\Psi_{i}(s)$, \big{|}{\rm init}\big{>}_{i}, and \big{|}{\rm eq_{M}}\big{>}_{i}, etc be consist of $\{B_{\alpha_{ir}}s^{\alpha_{ir}}\}$, $\{\lambda_{i,r}\}$, and $\{\varepsilon_{i,r}\}$ etc, respectively, where $r=1,2,\ldots,n_{i}$. After obtaining the inverse matrix of $A(s)$, we have the PDF $g(x,t)$ in the Fourier-Laplace space

[TABLE]

where $i_{m}=m$ if \big{|}{\rm init}\big{>}_{m}=0 otherwise $i_{m}=0$, $m=1,2,\cdots,j$. In a very particular case, i.e., the transition matrix is an identity matrix and $i_{m}=0$ for $m=1,2,\cdots,N$, the PDF of the process has the form

[TABLE]

which is the same as the PDF of (modified-form) distributed-order diffusion sandev:00 ; chechkin:00 with $p(\alpha)=\sum^{N}_{i=1}\varepsilon_{i}\delta(\alpha-\alpha_{i})$ (their physical backgrounds are completely different). From (4), there exists

[TABLE]

which behaves as \big{<}x^{2}(t)\big{>}\sim t^{\alpha^{\ast}} for large $t$ with \alpha^{\ast}=\max\limits_{\mathop{1\leqslant i\leqslant j}\limits_{i\neq i_{1},i_{2},\cdots,i_{j}}}\Big{\{}\min\limits_{1\leqslant r\leqslant n_{i}}\{\alpha_{ir}\}\Big{\}}, being confirmed by simulating the stochastic process (see Fig. 1); the influence of the distribution of the initial states is also observed.

Next, we consider the case that the transition matrix is not strictly the form of block diagonal matrix. Without loss of generality, we assume that the first $i$ rows of the transition matrix still keep the form of block diagonal matrix while the others not. If the elements of the initial distribution from $(i+1)$-th to $N$-th are [math] (denoting this part of the vector as \big{|}{\rm init}\big{>}_{i+1,N}), then the results illustrated by Fig. 1 still hold, just neglecting the last $(N-i)$ internal states; on the other hand, if \big{|}{\rm init}\big{>}_{i+1,N} does not equal to zero and, at the same time, the elements of \big{|}{\rm init}\big{>} from first to $i$-th are nonzero, the results for large $t$ can still be displayed by Fig. 1, since all the particles in the states ($(i+1),\cdots,N$) finally go to some of the first $i$ states. A little bit complex case is that as least one of the first $i$ elements of —init¿ is zero and \big{|}{\rm init}\big{>}_{i+1,N}\neq 0; in this case, taking \big{<}{\rm init}\big{|}M^{N-1} as the new initial distribution, the results are still depicted by Fig. 1, ignoring the last $(N-i)$ internal states.

Equations governing the distribution of the functionals of the paths and internal states of the process.—Two types of functionals will be considered. One is still defined as $A=\int_{0}^{t}U(x(\tau))d\tau$ Kac:49 , being widely discussed carmi:00 ; turgeman:00 ; cairoli:00 ; satya:00 ; Wu:16 , where $U(x)$ is a prescribed function and $x(t)$ is a trajectory of a particle. The other functional is first introduced here, defined as $A_{s}=\int^{t}_{0}U(i(\tau))d\tau$, where $i(\tau)$ represents that the particle is in the state $i$ at time $\tau$, naturally its values belonging to $\{1,2,\ldots,N\}$.

We use the notation $g^{(i)}(x,A,t),i=1,2,\ldots,N$ to represent the joint PDF of finding the particle at position $x$ with the functional $A$ and in the internal state of $E_{i}$ at time t. Gather all $g^{(i)}(x,A,t)$ to form a column vector denote by \big{|}G(x,A,t)\big{>}. Following the process of the derivation of the fractional Feynman-Kac equation carmi:00 , we have

[TABLE]

with the initial condition \big{|}G(x,A,t)\big{>}\,\big{|}_{t=0}=\delta(A)\delta(x)\delta(t)\big{|}{\rm init}\big{>}, and $\mathcal{D}_{t}^{{\rm diag}(1-\alpha_{1},\ldots,1-\alpha_{N})}={\rm diag}(\mathcal{D}_{t}^{1-\alpha_{1}},\ldots,\mathcal{D}_{t}^{1-\alpha_{N}})$ with $\mathcal{D}_{t}^{1-\alpha_{i}}$ being the fractional substantial derivative carmi:00 ; friedrich:00 . Next we derive the backward version of Eq. (5). We use the notation \big{|}{\rm init}\big{>}_{x_{0}}=\big{(}\lambda^{(1)}_{x_{0}},\lambda^{(2)}_{x_{0}},\ldots,\lambda^{(N)}_{x_{0}}\big{)} to represent the initial distribution of the internal states of the process starting at $x_{0}$, and $g_{x_{0}}^{(i)}(A,t)$ the PDF of the functional $A$ of the process at $t$, starting at $x_{0}$ with the internal state $E_{i}$. After some calculations, one can get

[TABLE]

If one is only interested in the PDF of $A$ at $t$ of the process starting at $x_{0}$, just calculate $g_{x_{0}}(A,t)=\sum_{i=1}^{N}\lambda_{x_{0}}^{(i)}g^{(i)}_{x_{0}}(A,t)$. Next, we give a specific application of (6) for calculating the distribution of the first passage time $t_{f}$, being the time that the particle starting at $x_{0}$ ($<B$) first reaches $x=B$. Define $A_{f}=\int_{0}^{t}U(x(\tau))d\tau$, where $U(x)=0$ if $x<B$ otherwise it equals to $1$. According to carmi:00 ; Wu:16 , there exists the relation Pr\{t_{f}>t\}=Pr\Big{\{}\max_{0\leq\tau<t}x(\tau)<B\Big{\}}=\lim_{\rho\rightarrow\infty}g_{x_{0}}(\rho,t). For the process of two internal states with the alternating transition matrix, and the coefficients $K_{\alpha_{1}}=K_{\alpha_{2}}=B_{\alpha_{1}}=B_{\alpha_{2}}=1$, we have

[TABLE]

where $a_{0}=-2+s^{\alpha_{1}}+s^{\alpha_{2}}$ and $b_{0}=\sqrt{4+s^{2\alpha_{1}}+s^{2\alpha_{2}}-2s^{\alpha_{1}+\alpha_{2}}}$; here $x_{0}$ is taken as [math] but not essential. Thus one can obtain the PDF of the first passage time

[TABLE]

When $t$ is big and $\alpha_{1}>\alpha_{2}$, then $a_{0}\sim-2+s^{\alpha_{2}}$ and $b_{0}\sim\sqrt{4+s^{2\alpha_{2}}}\sim 2+\frac{1}{4}s^{2\alpha_{2}}$. Thus, \exp\Big{(}-\sqrt{\frac{a_{0}+b_{0}}{2}}B\Big{)}\sim\exp\Big{(}-\sqrt{\frac{s^{\alpha_{2}}}{2}}B\Big{)}\sim 1-\frac{B}{\sqrt{2}}s^{\frac{s^{\alpha_{2}}}{2}}, i.e., $f(t)\sim\frac{B}{\sqrt{2}\mid\Gamma(-\alpha_{2}/2)\mid}t^{-\alpha_{2}/2-1}$. The result is confirmed by simulations given in Fig. 2.

We then turn to the distribution of $A_{s}=\int^{t}_{0}U(i(\tau))\tau$, ignoring the position $x$ of the particle. Denote $G^{(i)}(A_{s},t)$ as the joint PDF of finding the particle with the functional $A_{s}$ in the state $i$ at time $t$. Let \big{|}G(A_{s},t)\big{>}=\begin{pmatrix}G^{(1)}(A_{s},t),\ldots,G^{(N)}(A_{s},t)\end{pmatrix}^{T}. Then we get the governing equation

[TABLE]

where

[TABLE]

A direct application of (7) is to calculate the distribution of the fraction of the occupation time, i.e., the distribution of $t^{(i)}/t$, denoted as $l_{t^{(i)}/t}(x)$, where $t^{(i)}$ represents the occupation time of state $i$. Without loss of generality, we only consider the occupation time of the first state by letting $U(i(\tau))=1$ if $i(\tau)=1$, otherwise $U(i(\tau))=0$. Here we just present three results (in the case that the transition matrix is irreducible reducible:17 ): 1) if $\alpha_{1}<\alpha_{2}\leq\alpha_{3}\leq\cdots\leq\alpha_{N}$, then $l_{t^{(i)}/t}(x)\sim\delta(x-1)$; 2) if $\alpha_{1}$ is not the strictly smallest one, then $l_{t^{(i)}/t}(x)\sim\delta(x)$; 3) if $\alpha_{1}=\alpha_{2}=\cdots=\alpha_{m}$ ($2\leq m\leq N$) but smaller than other exponents, then

[TABLE]

from which the classic arcsine law can be recovered stefani:00 .

Conclusion.—We derive the Fokker-Planck equations as well as equations of the functionals of the paths and internal states of the fcP process with multiple internal states. Based on the equations, the MSD is analyzed, and the first passage time and fraction of occupation time are calculated. If the transition matrix of the internal states is reducible, the initial distribution of the internal states significantly influence the final results.