Boundary problems for the fractional and tempered fractional operators

TL;DR

This paper explores boundary conditions for fractional and tempered fractional PDEs modeling anomalous diffusion, emphasizing the importance of boundary data in the complement of the domain for well-posedness and physical relevance.

Contribution

It introduces and justifies new boundary conditions for space fractional PDEs based on stochastic and probabilistic insights, ensuring well-posedness and physical interpretability.

Findings

Boundary conditions depend on the solution outside the domain.

Well-posedness of fractional PDEs with generalized boundary conditions.

Properties of fractional operators discussed.

Abstract

For characterizing the Brownian motion in a bounded domain: $\Omega$, it is well-known that the boundary conditions of the classical diffusion equation just rely on the given information of the solution along the boundary of a domain; on the contrary, for the L\'evy flights or tempered L\'evy flights in a bounded domain, it involves the information of a solution in the complementary set of $\Omega$, i.e., $\mathbb{R}^n\backslash \Omega$, with the potential reason that paths of the corresponding stochastic process are discontinuous. Guided by probability intuitions and the stochastic perspectives of anomalous diffusion, we show the reasonable ways, ensuring the clear physical meaning and well-posedness of the partial differential equations (PDEs), of specifying `boundary' conditions for space fractional PDEs modeling the anomalous diffusion. Some properties of the operators are discussed, and the well-posednesses of the PDEs with generalized boundary conditions are proved.

Equations

$$\hat{p}(\mathbf{k},u)=\frac{1-\varphi (u)}{u}\cdot \frac{{\hat{p}}_0(\mathbf{k})}{1-\Psi (u,\mathbf{k})},$$

$${\int }_{{\mathbb{R}}^n}{e}^{i\mathbf{k}\cdot \mathbf{X}}p(\mathbf{X},t)\mathbf{dX}=\mathbf{E}(e^{i\mathbf{k}\cdot \mathbf{X}})=e^{t\Phi (\mathbf{k})},$$

$$\Phi (\mathbf{k})=i\mathbf{a}\cdot \mathbf{k}-\frac{1}{2}(\mathbf{k}\cdot \mathbf{bk})+{\int }_{{\mathbb{R}}^n\backslash \{0\}}\left[e^{i\mathbf{k}\cdot \mathbf{X}}-1-i(\mathbf{k}\cdot \mathbf{X}){\chi }_{\{|\mathbf{X}|<1\}}\right]\nu (d\mathbf{X});$$

$$\nu (d\mathbf{X})=\frac{\beta \Gamma (\frac{n+\beta }{2})}{2^{1-\beta }{\pi }^{n/2}\Gamma (1-\beta /2)}|\mathbf{X}{|}^{-\beta -n}\mathbf{dX},$$

$$\frac{\partial p(\mathbf{X},t)}{\partial t}={\Delta }^{\beta /2}p(\mathbf{X},t),$$

$$\frac{\partial p(\mathbf{X},t)}{\partial t}=(\Delta +\lambda {)}^{\beta /2}p(\mathbf{X},t),$$

$$\nu (d\mathbf{X})$$

$$+\frac{{\beta }_2\Gamma (\frac{1+{\beta }_2}{2})}{2^{1-{\beta }_2}{\pi }^{1/2}\Gamma (1-{\beta }_2/2)}\delta ({\mathbf{x}}_1)|{\mathbf{x}}_2{|}^{-{\beta }_2-1}d{\mathbf{x}}_{1}d{\mathbf{x}}_2,$$

$$\frac{\partial p({\mathbf{x}}_1,{\mathbf{x}}_2,t)}{\partial t}=\frac{{\partial }^{{\beta }_1}p({\mathbf{x}}_1,{\mathbf{x}}_2,t)}{\partial |{\mathbf{x}}_1{|}^{{\beta }_1}}+\frac{{\partial }^{{\beta }_2}p({\mathbf{x}}_1,{\mathbf{x}}_2,t)}{\partial |{\mathbf{x}}_2{|}^{{\beta }_2}}.$$

$$\frac{5}{2\sqrt{2C}}\frac{\Gamma (3-\beta )}{\Gamma (2-\beta {)}^{3/2}}{\lambda }^{-\frac{1}{2}\beta }\frac{1}{\sqrt{m}},$$

$$m\sim {\lambda }^{-\beta }.$$

$$|Q_m(\mathbf{X})-Q(\mathbf{X})|<\frac{5}{2}\frac{⟨|\mathbf{X}{|}^3⟩}{⟨|\mathbf{X}{|}^2{⟩}^{3/2}}\frac{1}{\sqrt{m}}=\frac{5}{2\sqrt{2C}}\frac{\Gamma (3-\beta )}{\Gamma (2-\beta {)}^{3/2}}{\lambda }^{-\frac{1}{2}\beta }\frac{1}{\sqrt{m}},$$

$$⟨|\mathbf{X}{|}^3⟩=C{\int }_{-\infty }^{\infty }|\mathbf{X}{|}^3{e}^{-\lambda |\mathbf{X}|}|\mathbf{X}{|}^{-\beta -1}d|\mathbf{X}|=2C{\int }_0^{\infty }{e}^{-\lambda |\mathbf{X}|}|\mathbf{X}{|}^{3-\beta -1}d|\mathbf{X}|=2C{\lambda }^{\beta -3}\Gamma (3-\beta )$$

$$⟨|\mathbf{X}{|}^2⟩=C{\int }_{-\infty }^{\infty }|\mathbf{X}{|}^2{e}^{-\lambda |\mathbf{X}|}|\mathbf{X}{|}^{-\beta -1}d|\mathbf{X}|=2C{\int }_0^{\infty }{e}^{-\lambda |\mathbf{X}|}|\mathbf{X}{|}^{2-\beta -1}d|\mathbf{X}|=2C{\lambda }^{\beta -2}\Gamma (2-\beta ).$$

$$\hat{p}(\mathbf{k},t)$$

$$=\sum _{j=0}^{\infty }\mathbf{E}(e^{i\mathbf{k}\cdot \mathbf{X}(t)}|N(t)=j)P(N(t)=j)$$

$$=\sum _{j=0}^{\infty }\mathbf{E}(e^{i\mathbf{k}\cdot ({\mathbf{X}}_0+{\mathbf{X}}_1+\cdots +{\mathbf{X}}_j)}|N(t)=j)P(N(t)=j)$$

$$=\sum _{j=0}^{\infty }{\Phi }_0(\mathbf{k}{)}^j\frac{(\zeta t{)}^j}{j!}{e}^{-\zeta t}$$

$$=e^{\zeta t({\Phi }_0(\mathbf{k})-1)},$$

$$\hat{p}(k,u)=\frac{1}{u+\zeta c^{\beta }|k{|}^{\beta }};$$

$$\hat{p}(k,u)=\frac{1}{u+\zeta C_{\beta ,\lambda }[(\lambda +ik{)}^{\beta }-{\lambda }^{\beta }]+\zeta C_{\beta ,\lambda }[(\lambda -ik{)}^{\beta }-{\lambda }^{\beta }]}.$$

$$\hat{p}(k,u)=\frac{{\tilde{c}}^{\alpha }}{u^{1-\alpha }(1-(1-{\tilde{c}}_{\alpha }{u}^{\alpha })\psi (k))}.$$

$${\hat{p}}_x(\mathbf{k},t)$$

$$={\int }_0^{\infty }{\int }_{{\mathbb{R}}^n}{e}^{i\mathbf{k}\cdot \mathbf{Y}}{p}_y(\mathbf{Y},\tau )d\mathbf{Y}{p}_s(\tau ,t)d\tau$$

$$={\int }_0^{\infty }{e}^{-\tau (-{\Phi }_y(\mathbf{k}))}{p}_s(\tau ,t)d\tau$$

$$=e^{-t{\Phi }_s(-{\Phi }_y(\mathbf{k}))},$$

$$p_e(\tau ,t)=\frac{\partial P(E(t)\le \tau )}{\partial \tau }=\frac{\partial }{\partial \tau }\left[1-P(S(\tau )<t)\right]=-\frac{\partial }{\partial \tau }{\int }_0^t{p}_s(y,\tau )dy.$$

$$p_e(\tau ,u)=-\frac{\partial }{\partial \tau }{u}^{-1}{e}^{-\tau u^{\alpha }}=u^{\alpha -1}{e}^{-\tau u^{\alpha }}.$$

$$p_x(\mathbf{X},t)={\int }_0^{\infty }{p}_y(\mathbf{X},\tau )p_e(\tau ,t)d\tau .$$

$${\hat{p}}_x(\mathbf{k},u)$$

$$={\int }_0^{\infty }{e}^{-\tau {\Phi }_y(\mathbf{k})}{u}^{\alpha -1}{e}^{-\tau u^{\alpha }}d\tau$$

$$=\frac{u^{\alpha -1}}{u^{\alpha }+{\Phi }_y(\mathbf{k})}.$$

$$\{\begin{array}{rl}& \dot{\mathbf{X}}(\tau )=F(\mathbf{X}(\tau ))+\eta (\tau ),\\ & \dot{S}(\tau )=\xi (\tau ),\end{array}$$

$$\frac{d\hat{p}(\mathbf{k},t)}{dt}=-|\mathbf{k}{|}^{\beta }\hat{p}(\mathbf{k},t).$$

$$\frac{\partial p(\mathbf{X},t)}{\partial t}={\Delta }^{\beta /2}p(\mathbf{X},t),$$

Full text

1 Introduction

The phrase ‘anomalous is normal’ says that anomalous diffusion phenomena are ubiquitous in the natural world. It was first used in the title of [24], which reveals that the diffusion of classical particles on a solid surface has rich anomalous behaviour controlled by the friction coefficient. In fact, anomalous diffusion is no longer a young topic. In the review paper [5], the evolution of particles in disordered environments was investigated; the specific effects of a bias on anomalous diffusion were considered; and the generalizations of Einstein’s relation in the presence of disorder were discussed. With the rapid development of the study of anomalous dynamics in diverse field, some deterministic equations are derived, governing the macroscopic behaviour of anomalous diffusion. In 2000, Metzler and Klafter published the survey paper [22] for the equations governing transport dynamics in complex system with anomalous diffusion and non-exponential relaxation patterns, i.e., fractional kinetic equations of the diffusion, advection-diffusion, and Fokker-Planck type, derived asymptotically from basic random walk models and a generalized master equation. Many mathematicians have been involved in the research of fractional partial differential equations (PDEs). For fractional PDEs in a bounded domain $\Omega$, an important question is how to introduce physically meaningful and mathematically well-posed boundary conditions on $\partial\Omega$ or $\mathbb{R}^{n}\backslash\Omega$.

Microscopically, diffusion is the net movement of particles from a region of higher concentration to a region of lower concentration; for the normal diffusion (Brownian motion), the second moment of the particle trajectories is a linear function of the time $t$; naturally, if it is a nonlinear function of $t$, we call the corresponding diffusion process anomalous diffusion or non-Brownian diffusion [22]. The microscopic (stochastic) models describing anomalous diffusion include continuous time random walks (CTRWs), Langevin type equation, Lévy processes, subordinated Lévy processes, and fractional Brownian motion, etc.. The CTRWs contain two important random variables describing the motion of particles [23], i.e., the waiting time ${\xi}$ and jump length $\eta$. If both the first moment of ${\xi}$ and the second moment of $\eta$ are finite in the scaling limit, then the CTRWs approximate Brownian motion. On the contrary, if one of them is divergent, then the CTRWs characterize anomalous diffusion. Two of the most important CTRW models are Lévy flights and Lévy walks. For Lévy flights, the ${\xi}$ with finite first moment and $\eta$ with infinite second moment are independent, leading to infinite propagation speed and the divergent second moments of the distribution of the particles. This causes much difficulty in relating the models to experimental data, especially when analyzing the scaling of the measured moments in time [30]. With coupled distribution of ${\xi}$ and $\eta$ (the infinite speed is penalized by the corresponding waiting times), we get the so-called Lévy walks [30]. Another idea to ensure that the processes have bounded moments is to truncate the long tailed probability distribution of Lévy flights [19]; they still look like a Lévy flight in not too long a time. Currently, the most popular way to do the truncation is to use the exponential tempering, offering the technical advantage of still being an infinitely divisible Lévy process after the operation [21]. The Lévy process to describe anomalous diffusion is the scaling limit of CTRWs with independent ${\xi}$ and $\eta$. It is characterized by its characteristic function. Except Brownian motion with drift, the paths of all other proper Lévy processes are discontinuous. Sometimes, the Lévy flights are conveniently described by the Brownian motion subordinated to a Lévy process [6]. Fractional Brownian motions are often taken as the models to characterize subdiffusion [18].

Macroscopically, fractional (nonlocal) PDEs are the most popular and effective models for anomalous diffusion, derived from the microscopic models. The solution of fractional PDEs is generally the probability density function (PDF) of the position of the particles undergoing anomalous dynamics; with the deepening of research, the fractional PDEs governing the functional distribution of particles’ trajectories are also developed [28, 29]. Two ways are usually used to derive the fractional PDEs. One is based on the Montroll-Weiss equation [23], i.e., in Fourier-Laplace space, the PDF $p({\bf X},t)$ obeys

[TABLE]

where $\hat{p}_{0}({\bf k})$ is the Fourier transform of the initial data; $\phi(u)$ is the Laplace transform of the PDF of waiting times ${\xi}$ and $\Psi(u,{\bf k})$ the Laplace and the Fourier transforms of the joint PDF of waiting times ${\xi}$ and jump length $\eta$. If ${\xi}$ and $\eta$ are independent, then $\Psi(u,{\bf k})=\phi(u)\psi({\bf k})$, where $\psi({\bf k})$ is the Fourier transform of the PDF of $\eta$. Another way is based on the characteristic function of the $\alpha$-stable Lévy motion, being the scaling limit of the CTRW model with power law distribution of jump length $\eta$. In the high dimensional case, it is more convenient to make the derivation by using the characteristic function of the stochastic process. According to the Lévy-Khinchin formula [2], the characteristic function of Lévy process has a specific form

[TABLE]

where

[TABLE]

here $\chi_{I}$ is the indicator function of the set $I$, ${\bf a}\in\mathbb{R}^{n}$, ${\bf b}$ is a positive definite symmetric $n\times n$ matrix and $\nu$ is a sigma-finite Lévy measure on $\mathbb{R}^{n}\backslash\{0\}$. When ${\bf a}$ and ${\bf b}$ are zero and

[TABLE]

the process is a rotationally symmetric $\beta$-stable Lévy motion and its PDF solves

[TABLE]

where $\mathcal{F}(\Delta^{\beta/2}p({\bf X},t))=-|{\bf k}|^{\beta}\mathcal{F}(p({\bf X},t))$ [26]. If replacing (3) by the measure of isotropic tempered power law with the tempering exponent $\lambda$, then we get the corresponding PDF evolution equation

[TABLE]

where $(\Delta+\lambda)^{\beta/2}$ is defined by (32) in physical space and by (34) in Fourier space.

In practice, the choice of $\nu(d{\bf X})$ depends strongly on the concrete physical environment. For example, Figure 1 clearly shows the horizontal and vertical structure. So, we need to take the measure as (if it is superdiffusion)

[TABLE]

where $\beta_{1}$ and $\beta_{2}$ belong to $(0,2)$. If ${\bf a}$ and ${\bf b}$ equal to zero, then it leads to diffusion equation

[TABLE]

Under the guidelines of probability intuitions and stochastic perspectives [15] of Lévy flights or tempered Lévy flights, we discuss the reasonable ways of defining fractional partial differential operators and specifying the ‘boundary’ conditions for their macroscopic descriptions, i.e., the PDEs of the types Eqs. (4), (5), (7), and their extensions, e.g., the fractional Feynman-Kac equations [28, 29]. For the related discussions on the nonlocal diffusion problems from a mathematical point of view, one can see the review paper [10]. The divergence of the second moment and the discontinuity of the paths of Lévy flights predicate that the corresponding diffusion operators should defined on $\mathbb{R}^{n}$, which further signify that if we are solving the equations in a bounded domain $\Omega$, the information in $\mathbb{R}^{n}\backslash\Omega$ should also be involved. We will show that the generalized Dirichlet type boundary conditions should be specified as $p({\bf X},t)|_{\mathbb{R}^{n}\backslash\Omega}=g({\bf X},t)$. If the particles are killed after leaving the domain $\Omega$, then $g({\bf X},t)\equiv 0$, i.e., the so-called absorbing boundary conditions. Because of the discontinuity of the jumps of Lévy flights, a particular concept ‘escape probability’ can be introduced, which means the probability that the particle jumps from the domain $\Omega$ into a domain $H\subset\mathbb{R}^{n}\backslash\Omega$; for solving the escape probability, one just needs to specify $g({\bf X})=1$ for ${\bf X}\in H$ and [math] for ${\bf X}\in({\mathbb{R}}^{n}\backslash\Omega)\backslash H$ for the corresponding time-independent PDEs. As for the generalized Neumann type boundary conditions, our ideas come from the fact that the continuity equation (conservation law) holds for any kinds of diffusion, since the particles can not be created or destroyed. Based on the continuity equation and the governing equation of the PDF of Lévy or tempered Lévy flights, the corresponding flux ${\bf j}$ can be obtained. So the generalized reflecting boundary conditions should be ${\bf j}|_{\mathbb{R}^{n}\backslash\Omega}\equiv 0$, which implies $(\nabla\cdot{\bf j})|_{\mathbb{R}^{n}\backslash\Omega}\equiv 0$. Then, the generalized Neumann type boundary conditions are given as $(\nabla\cdot{\bf j})|_{\mathbb{R}^{n}\backslash\Omega}=g({\bf X},t)$, e.g., for (4), it should be taken as $\left(\Delta^{\beta/2}p({\bf X},t)\right)|_{\mathbb{R}^{n}\backslash\Omega}=g({\bf X},t)$. The well-posednesses of the equations under our specified generalized Dirichlet or Neumann type boundary conditions are well established.

Overall, this paper focuses on introducing physically reasonable boundary constraints for a large class of fractional PDEs, building a bridge between the physical and mathematical communities for studying anomalous diffusion and fractional PDEs. In the next section, we recall the derivation of fractional PDEs. Some new concepts are introduced, such as the tempered fractional Laplacian, and some properties of anomalous diffusion are found. In Sec. 3, we discuss the reasonable ways of specifying the generalized boundary conditions for the fractional PDEs governing the position or functional distributions of Lévy flights and tempered Lévy flights. In Sec. 4, we prove well-posedness of the fractional PDEs under the generalized Dirichlet and Neumann boundary conditions defined on the complement of the bounded domain. Conclusion and remarks are given in the last section.

2 Preliminaries

For well understanding and inspiring the ways of specifying the ‘boundary constrains’ to PDEs governing the PDF of Lévy flights or tempered Lévy flights, we will show the ideas of deriving the microscopic and macroscopic models.

2.1 Microscopic models for anomalous diffusion

For the microscopic description of the anomalous diffusion, we consider the trajectory of a particle or a stochastic process, i.e., ${\bf X}(t)$. If $\left<|{\bf X}(t)|^{2}\right>\sim t$, the process is normal, otherwise it is abnormal. The anomalous diffusions of most often happening in natural world are the cases that $\left<|{\bf X}(t)|^{2}\right>\sim t^{\gamma}$ with $\gamma\in[0,1)\cup(1,2]$. A Lévy flight is a random walk in which the jump length has a heavy tailed (power law) probability distribution, i.e., the PDF of jump length $r$ is like $r^{-\beta-n}$ with $\beta\in(0,2)$, and the distribution in direction is uniform. With the wide applications of Lévy flights in characterizing long-range interactions [3] or a nontrivial “crumpled” topology of a phase (or configuration) space of polymer systems [27], etc, its second and higher moments are divergent, leading to the difficulty in relating models to experimental data. In fact, for Lévy flights $\left<|{\bf X}(t)|^{\delta}\right>\sim t^{\delta/\beta}$ with $0<\delta<\beta\leq 2$. Under the framework of CTRW, the model Lévy walk [25] can circumvent this obstacle by putting a larger time cost to a longer displacement, i.e., using the space-time coupled jump length and waiting time distribution $\Psi(r,t)=\frac{1}{2}\delta(r-vt)\phi(t)$. Another popular model is the so-called tempered Lévy flights [16], in which the extremely long jumps is exponentially cut by using the distribution of jump length $e^{-r\lambda}r^{-\beta-n}$ with $\lambda$ being a small modulation parameter (a smooth exponential regression towards zero). In not too long a time, the tempered Lévy flights display the dynamical behaviors of Lévy flights, ultraslowly converging to the normal diffusion. Figure 2 shows the trajectories of $1000$ steps of Lévy flights, tempered Lévy flights, and Brownian motion in two dimensions; note the presence of rare but large jumps compared to the Brownian motion, playing the dominant role in the dynamics.

Using Berry-Esséen theorem [12], first established in 1941, which applies to the convergence to a Gaussian for a symmetric random walk whose jump probabilities have a finite third moment, we have that for the one dimensional tempered Lévy flights with the distribution of jump length $Ce^{-r\lambda}r^{-\beta-1}$ the convergence speed is

[TABLE]

which means that the scaling law for the number of steps needed for Gaussian behavior to emerge as

[TABLE]

More concretely, letting ${\bf X}_{1}$, ${\bf X}_{2}$, $\cdots$, ${\bf X}_{m}$ be i.i.d. random variables with PDF $Ce^{-r\lambda}r^{-\beta-1}$ and $E(|{\bf X}_{1}|^{2})=\sigma^{2}>0$, then the cumulative distribution function (CDF) $Q_{m}$ of ${\bf Y}_{m}=({\bf X}_{1}+{\bf X}_{2}+\cdots+{\bf X}_{m})/(\sigma\sqrt{m})$ converges to the CDF $Q({\bf X})$ of the standard normal distribution as

[TABLE]

since

[TABLE]

and

[TABLE]

From Eq. (8), it can be seen that with the decrease of $\lambda$, the required $m$ for the crossover between Lévy flight behavior and Gaussian behavior increase rapidly. A little bit counterintuitive observation is that the number of variables required to the crossover increases with the increase of $\beta$.

We have described the distributions of jump length for Lévy flights and tempered Lévy flights, in which Poisson process is taken as the renewal process. We denote the Poisson process with rate $\zeta>0$ as $N(t)$ and its waiting time distribution between two events is $\zeta e^{-\zeta t}$. Then the Lévy flights or tempered Lévy flights are the compound Poisson process defined as ${\bf X}(t)=\sum\limits_{j=0}^{N(t)}{\bf X}_{j},$ where ${\bf X}_{j}$ are i.i.d. random variables with the distribution of power law or tempered power law. The characteristic function of ${\bf X}(t)$ can be calculated as follows. For real ${\bf k}$, we have

[TABLE]

where $\Phi_{0}({\bf k})={\bf E}(e^{i{\bf k}\cdot{\bf X}_{0}})$, being also the characteristic function of ${\bf X_{1}}$, ${\bf X_{2}}$, $\cdots$, ${\bf X_{j}}$ since they are i.i.d.

In the CTRW model describing one dimensional Lévy flights or tempered Lévy flights, the PDF of waiting times is taken as $\zeta e^{-\zeta t}$ with its Laplace transform $\zeta/(u+\zeta)$ and the PDF of jumping length is $c^{-\beta}r^{-\beta-1}$ or $e^{-\lambda r}r^{-\beta-1}$ with its Fourier transform $1-c^{\beta}|k|^{\beta}$ or $1-c_{\beta,\lambda}[(\lambda+ik)^{\beta}-\lambda^{\beta}]-c_{\beta,\lambda}[(\lambda-ik)^{\beta}-\lambda^{\beta}]$. Substituting them into the Montroll-Weiss Eq. (1) with $\hat{p}_{0}(k)=1$ (the initial position of particles is at zero), we get that $\hat{p}(k,u)$ of Lévy flights solves

[TABLE]

and the $\hat{p}(k,u)$ of tempered Lévy flights obeys

[TABLE]

If the subdiffusion is involved, we need to choose the PDF of waiting times as $\tilde{c}^{1+\alpha}t^{-\alpha-1}$ with $\alpha\in(0,1)$ and its Laplace transform $1-\tilde{c}^{\alpha}u^{\alpha}$. Then from (1), we get that

[TABLE]

For high dimensional case, the Lévy flights can also be characterized by Brownian motion subordinated to a Lévy process. Let ${\bf Y}(t)$ be a Brownian motion with Fourier exponent $-|{\bf k}|^{2}$ and $S(t)$ a subordinator with Laplace exponent $u^{\beta/2}$ that is independent of ${\bf Y}(t)$. The process ${\bf X}(t)={\bf Y}({S(t)})$ is describing Lévy flights with Fourier exponent $-|{\bf k}|^{\beta}$, being the subordinate process of ${\bf Y}(t)$. In effect, denote the characteristic function of ${\bf Y}(t)$ as $\Phi_{y}({\bf k})$ and the one of $S(t)$ as $\Phi_{s}(u)$. Then the characteristic function of ${\bf X}(t)$ is as follows:

[TABLE]

where $p_{x}$, $p_{y}$, and $p_{s}$, are respectively the PDFs of the stochastic processes ${\bf X}$, ${\bf Y}$, and $S$. Similarly, in the following, we denote $p$ with subscript (lowercase letter) as the PDF of the corresponding stochastic process (uppercase letter).

This paper mainly focuses on Lévy flights and tempered Lévy flights. If one is interested in subdiffusion, instead of Poisson process, the fractional Poisson process should be taken as the renewal process, in which the time interval between each pair of events follows the power law distribution. Let ${\bf Y}(t)$ be a general Lévy process with Fourier exponent $\Phi_{y}({\bf k})$ and $S(t)$ a strictly increasing subordinator with Laplace exponent $u^{\alpha}$ ($\alpha\in(0,1)$). Define the inverse subordinator $E(t)=\inf\{\tau>0:\,S(\tau)>t\}$. Since $t=S(\tau)$ and $\tau=E(t)$ are inverse processes, we have $P(E(t)\leq\tau)=P(S(\tau)\geq t)$. Hence

[TABLE]

In the above equation, taking Laplace transform w.r.t $t$ leads to

[TABLE]

For the PDF $p_{x}({\bf X},t)$ of ${\bf X}(t)={\bf Y}(E(t))$, there holds

[TABLE]

Performing Fourier transform w.r.t. ${\bf X}$ and Laplace transform w.r.t. $t$ to the above equation results in

[TABLE]

Remark. According to Fogedby [14], the stochastic trajectories of (scale limited) CTRW ${\bf X}(E_{t})$ can also be expressed in terms of the coupled Langevin equation

[TABLE]

where $F({\bf X})$ is a vector field; $E_{t}$ is the inverse process of $S(t)$; the noises $\eta(\tau)$ and $\xi(\tau)$ are statistically independent, corresponding to the distributions of jump length and waiting times.

2.2 Derivation of the macroscopic description from the microscopic models

This section focuses on the derivation of the deterministic equations governing the PDF of position of the particles undergoing anomalous diffusion. It shows that the operators related to (tempered) power law jump lengths should be defined on the whole unbounded domain $\mathbb{R}^{n}$, which can also be inspired by the rare but extremely long jump lengths displayed in Figure 2; the fact that among all proper Lévy processes Brownian motion is the unique one with continuous paths further consolidates the reasonable way of defining the operators. We derive the PDEs based on Eqs. (9), (13), and (16), since they apply for both one and higher dimensional cases. For one dimensional case, sometimes it is convenient to use (10), (11), and (12).

When the diffusion process is rotationally symmetric $\beta$-stable, i.e., it is isotropic with PDF of jump length $c_{\beta,n}r^{-\beta-n}$ and its Fourier transform $1-{|\bf k|}^{\beta}$, where $n$ is the space dimension. In Eq. (9), taking $\zeta$ equal to $1$, we get the Cauchy equation

[TABLE]

Performing inverse Fourier transform to the above equation leads to

[TABLE]

where

[TABLE]

with [8]

[TABLE]

For the more general cases of Eq. (9), there is the Cauchy equation

[TABLE]

so the PDF of the stochastic process ${\bf X}$ solves (taking $\zeta=1$)

[TABLE]

where $\nu(d{\bf Y})$ is the probability measure of the jump length. Sometimes, to overcome the possible divergence of the terms on the right hand side of Eq. (24) because of the possible strong singularity of $\nu(d{\bf Y})$ at zero, the term

[TABLE]

is approximately replaced by

[TABLE]

then the corresponding modification to Eq. (24) is

[TABLE]

where ${\bf y}_{i}$ is the component of ${\bf Y}$, i.e., ${\bf Y}=\{{\bf y}_{1},{\bf y}_{2},\cdots,{\bf y}_{n}\}^{T}$. If $\nu(-d{\bf Y})=\nu(d{\bf Y})$, the integration of the summation term of above equation equals to zero.

If the diffusion is in the environment having a structure like Figure 1, the probability measure should be taken as

[TABLE]

where $\beta_{1},\,\beta_{2},\cdots,\,\beta_{n}$ belong to $(0,2)$. Plugging Eq. (27) into Eq. (24) leads to

[TABLE]

where

[TABLE]

and $\frac{\partial^{\beta_{j}}p({\bf x}_{1},\cdots,{\bf x}_{n},t)}{\partial|{\bf x}_{j}|^{\beta_{j}}}$ in physical space is defined by (21) with $n=1$; in particular, when $\beta_{j}\in(1,2)$, it can also be written as

[TABLE]

It should be emphasized here that when characterizing diffusion processes related with Lévy flights the operators should be defined in the whole space. Another issue that also should be stressed is that when $\beta_{1}=\beta_{2}=\cdots=\beta_{n}=1$, Eq. (28) is still describing the phenomena of anomalous diffusion, including the cases that they belong to $(0,1)$; the corresponding ‘first’ order operator is nonlocal, being different from the classical first order operator, but they have the same energy in the sense that

[TABLE]

[TABLE]

even though $\Delta^{1/2}$ and $\nabla$ are completely different operators, where the notation $\overline{v}$ stands for the complex conjugate of $v$.

If the subdiffusion is involved, the derivation of the macroscopic equation should be based on Eq. (17). For getting the term related to time derivative, the inverse Laplace transform should be performed on $u^{\alpha}\hat{p}({\bf k},u)-u^{\alpha-1}$. Since $\hat{p}({\bf k},t=0)$ is taken as $1$, there exists

[TABLE]

which is usually denoted as ${{}_{0}^{C}D_{t}^{\alpha}}\hat{p}({\bf k},t)$, the so-called Caputo fractional derivative. So, if both the PDFs of the waiting time and jump lengths of the stochastic process ${\bf X}$ are power law, the corresponding models can be obtained by replacing $\frac{\partial}{\partial t}$ with ${{}_{0}^{C}D_{t}^{\alpha}}$ in Eqs. (20), (24), (26), and (28). Furthermore, if there is an external force $F(\bf X)$ in the considered stochastic process $\bf X$, we need to add an additional term $\nabla\cdot(F({\bf X})p({\bf X},t))$ on the right hand side of Eqs. (20), (24), (26), and (28).

Here we turn to another important and interesting topic: tempered Lévy flights. Practically it is not easy to collect the value of a function in the unbounded area ${\mathbb{R}}^{n}\backslash\Omega$. This is one of the achievements of using tempered fractional Laplacian. It is still isotropic but with PDF of jump length $c_{\beta,n,\lambda}e^{-\lambda r}r^{-\beta-n}$. The PDF of tempered Lévy flights solves

[TABLE]

where

[TABLE]

with

[TABLE]

The choice of the constant as the one given in (33) leads to

[TABLE]

However, if $\lambda=0$, one needs to choose the constant as the one given in (22) to make sure $\mathcal{F}\left(\Delta^{\beta/2}p({\bf X},t)\right)=-|{\bf k}|^{\beta}\hat{p}({\bf k},t)$. The reason is as follows.

[TABLE]

For $\beta\in(0,1)\cup(1,2)$, then we have

[TABLE]

where ${{}_{2}}F_{1}$ is the Gaussian hypergeometric function and

[TABLE]

So

[TABLE]

The PDEs for tempered Lévy flights or tempered Lévy flights combined with subdiffusion can be similarly derived, as those done in this section for Lévy flights or Lévy flights combined with subdiffusion. Here, we present the counterpart of Eq. (28),

[TABLE]

where the operator $\frac{\partial^{\beta_{j},\lambda}p({\bf x}_{1},\cdots,{\bf x}_{j},t)}{\partial|{\bf x}_{j}|^{\beta_{j},\lambda}}$ is defined by taking $\beta=\beta_{j}$ and $n=1$ in Eq. (32). Again, even for the tempered Lévy flights, all the related operators should be defined on the whole space, because of the very rare but still possible unbounded jump lengths.

All the above derived PDEs are governing the PDF of the position of particles. If one wants to dig out more deep informations of the corresponding stochastic processes, analyzing the distribution of the functional defined by $A=\int_{0}^{t}U({\bf X(\tau)})d\tau$ is one of the choices, where $U$ is a prespecified function. Denote the PDF of the functional $A$ and position ${\bf X}$ as $G({\bf X},A,t)$ and the counterpart of $A$ in Fourier space as $q$. Then $\hat{G}({\bf X},q,t)$ solves [28]

[TABLE]

for Lévy flights combined with subdiffusion; and [29]

[TABLE]

for tempered Lévy flights combined with subdiffusion, where

[TABLE]

If one is only interested in the functional $A$ (not caring position ${\bf X}$), then $\hat{G}_{{\bf X}_{0}}(q,t)$ is, respectively, governed by [28]

[TABLE]

and [29]

[TABLE]

for Lévy flights and tempered Lévy flights, combined with subdiffusion; the ${\bf X}_{0}$ in $\hat{G}_{{\bf X}_{0}}(q,t)$ means the initial position of particles, being a parameter.

3 Specifying the generalized boundary conditions for the fractional PDEs

After introducing the microscopic models and deriving the macroscopic ones, we have insight into anomalous diffusions, especially Lévy flights and tempered Lévy flights. In Section 2, all the derived equations are time dependent. From the process of derivation, one can see that the issue of initial condition can be easily/reasonably fixed, as classical ones, just specifying the value of $p({\bf X},0)$ in the domain $\Omega$. For Lévy processes, except Brownian motion, all others have discontinuous paths. As a result, the boundary $\partial\Omega$ itself (see Figure 3) can not be hit by the majority of discontinuous sample trajectories. This implies that when solving the PDEs derived in Section 2, the generalized boundary conditions must be introduced, i.e., the information of $p({\bf X},t)$ on the domain ${\mathbb{R}}^{n}\backslash\Omega$ must be properly accounted for.

In the following, we focus on Eqs. (20), (28), (31), (35) to discuss the boundary issues.

3.1 Generalized Dirichlet type boundary conditions

The appropriate initial and boundary value problems for Eq. (20) should be

[TABLE]

In Eq. (40), the term

[TABLE]

According to Eq. (41), $g({\bf X},t)$ should satisfy that there exist positive $M$ and $C$ such that when $|{\bf X}|>M$,

[TABLE]

In particular, when Eq. (42) holds, the function $\int_{\mathbb{R}^{n}\backslash\Omega}\frac{-g({\bf Y},t)}{|{\bf X}-{\bf Y}|^{n+\beta}}d{\bf Y}$ of ${\bf X}$ has any order of derivative if $g({\bf X},t)$ is integrable in any bounded domain. One of the most popular cases is $g({\bf X},t)\equiv 0$, which is the so-called absorbing boundary condition, implying that the particle is killed whenever it leaves the domain $\Omega$. Another interesting case is for the steady state fraction diffusion equation

[TABLE]

Given a domain $H\subset\mathbb{R}^{n}\backslash\Omega$, if taking $g({\bf X})=1$ for ${\bf X}\in H$ and [math] for ${\bf X}\in(\mathbb{R}^{n}\backslash\Omega)\backslash H$, then the solution of (43) means the probability that the particles undergoing Lévy flights lands in $H$ after first escaping the domain $\Omega$ [7]. If $g({\bf X})\equiv 1$ in $\mathbb{R}^{n}\backslash\Omega$, then $p({\bf X})$ equals to $1$ in $\Omega$ because of the probability interpretation. This can also be analytically checked.

For the initial and boundary value problem Eq. (28), it should be written as

[TABLE]

Similar to (41), in (44) the term

[TABLE]

From Eq. (45), for $j=1,\cdots,n$, $g({\bf x}_{1},\cdots,{\bf x}_{j},\cdots,{\bf x}_{n},t)$ should satisfies that there exist positive $M$ and $C$ such that when $|{\bf x}_{j}|>M$,

[TABLE]

The discussions below Eq. (43) still makes sense for Eq. (44). If $g({\bf x}_{1},\cdots,{\bf x}_{j},\cdots,{\bf x}_{n},t)$ satisfies Eq. (46), and it is integrable w.r.t. ${\bf x}_{j}$ in any bounded interval. Then $\int_{\mathbb{R}\backslash(\Omega\cap{\mathbb{R}_{j}})}\frac{-g({\bf x}_{1},\cdots,{\bf y}_{j},\cdots,{\bf x}_{n},t)}{|{\bf x}_{j}-{\bf y}_{j}|^{1+\beta_{j}}}d{\bf y}_{j}$ has any order of partial derivative w.r.t. ${\bf x}_{j}$.

The initial and boundary value problem for Eq. (31) is

[TABLE]

Like the discussions for Eq. (40), $g({\bf X},t)$ should satisfies that there exist positive $M$ and $C$ such that when $|{\bf X}|>M$,

[TABLE]

If Eq. (48) holds and $g({\bf X},t)$ is integrable in any bounded domain, the function $\int_{\mathbb{R}^{n}\backslash\Omega}\frac{-g({\bf Y},t)}{e^{\lambda|{\bf X}-{\bf Y}|}|{\bf X}-{\bf Y}|^{n+\beta}}d{\bf Y}$ of ${\bf X}$ has any order of derivative.

Again, the corresponding tempered steady state fraction diffusion equation is

[TABLE]

For $H\subset\mathbb{R}^{n}\backslash\Omega$, if taking $g({\bf X})=1$ for ${\bf X}\in H$ and [math] for ${\bf X}\in(\mathbb{R}^{n}\backslash\Omega)\backslash H$, then the solution of (49) means the probability that the particles undergoing tempered Lévy flights lands in $H$ after first escaping the domain $\Omega$. If $g({\bf X})\equiv 1$ in $\mathbb{R}^{n}\backslash\Omega$, then $p({\bf X})$ equals to $1$ in $\Omega$.

The initial and boundary value problem (35) should be written as

[TABLE]

For $j=1,\cdots,n$, $g({\bf x}_{1},\cdots,{\bf x}_{j},\cdots,{\bf x}_{n},t)$ should satisfy that there exist positive $M$ and $C$ such that when $|{\bf x}_{j}|>M$,

[TABLE]

If $g({\bf x}_{1},\cdots,{\bf x}_{j},\cdots,{\bf x}_{n},t)$ is integrable w.r.t. ${\bf x}_{j}$ in any bounded interval and satisfies Eq. (51), then $\int_{\mathbb{R}\backslash(\Omega\cap{\mathbb{R}_{j}})}\frac{-g({\bf x}_{1},\cdots,{\bf y}_{j},\cdots,{\bf x}_{n},t)}{e^{\lambda|{\bf x}_{j}-{\bf y}_{j}}|{\bf x}_{j}-{\bf y}_{j}|^{1+\beta_{j}}}d{\bf y}_{j}$ has any order of partial derivative w.r.t. ${\bf x}_{j}$.

The ways of specifying the initial and boundary conditions for Eqs. (36) and (38) are the same as Eq. (40). But for Eq. (36), the corresponding (42) should be changed as

[TABLE]

Similarly, the initial and boundary conditions of Eqs. (37) and (39) should be specified as the ones of Eq. (47). But for Eq. (37), the corresponding (48) needs to be changed as

[TABLE]

For the existence and uniqueness of the corresponding time-independent equations, one may refer to [13].

3.2 Generalized Neumann type boundary conditions

Because of the inherent discontinuity of the trajectories of Lévy flights or tempered Lévy flights, the traditional Neumann type boundary conditions can not be simply extended to the fractional PDEs. For the related discussions, see, e.g., [4, 9]. Based on the models built in Sec. 2 and the law of mass conservation, we derive the reasonable ways of specifying the Neumann type boundary conditions, especially the reflecting ones. Let us first recall the derivation of classical diffusion equation. For normal diffusion (Brownian motion), microscopically the first moment of the distribution of waiting times and the second moment of the distribution of jump length are bounded, i.e., in Laplace and Fourier spaces, they are respectively like $1-c_{1}u$ and $1-c_{2}|{\bf k}|^{2}$; plugging them into Eq. (1) or Eq. (9) and performing integral transformations lead to the classical diffusion equation

[TABLE]

On the other hand, because of mass conservation, the continuity equation states that a change in density in any part of a system is due to inflow and outflow of particles into and out of that part of system, i.e., no particles are created or destroyed:

[TABLE]

where ${\bf j}$ is the flux of diffusing particles. Combining (54) with (55), one may take

[TABLE]

which is exactly Fick’s law, a phenomenological postulation, saying that the flux goes from regions of high concentration to regions of low concentration with a magnitude proportional to the concentration gradient. In fact, for a long history, even up to now, most of the people are more familiar with the process: using the continuity equation (55) and Fick’s law (56) derives the diffusion equation (54). The so-called reflecting boundary condition for (54) is to let the flux ${\bf j}$ be zero along the boundary of considered domain.

Here we want to stress that Eq. (55) holds for any kind of diffusions, including the normal and anomalous ones. For Eqs. (40,44,47,50) governing the PDF of Lévy flights or tempered Lévy flights, using the continuity equation (55), one can get the corresponding fluxes and the counterparts of Fick’s law; may we call it fractional Fick’s law. Combining (40) with (55), one may let

[TABLE]

being the flux for the diffusion operator $\Delta^{\beta/2}$ with $\beta\in(0,2)$, or calling it fractional Fick’s law corresponding to $\Delta^{\beta/2}$. From (44) and (55), one may choose

[TABLE]

where ${\bf\widetilde{Y}}_{i}=\{{\bf x}_{1},\dots,{\bf y}_{i},\cdots,{\bf x}_{n}\}^{T}$, being the flux (fractional Fick’s law) corresponding to the horizontal and vertical type fractional operators. Similarly, we can also get the flux (fractional Fick’s law) corresponding to the tempered fractional Laplacian and tempered horizontal and vertical type fractional operators, being respectively taken as

[TABLE]

and

[TABLE]

with ${\bf\widetilde{Y}}_{i}=\{{\bf x}_{1},\dots,{\bf y}_{i},\cdots,{\bf x}_{n}\}^{T}$.

Naturally, the Neumann type boundary conditions of (40,44,47,50) should be closely related to the values of the fluxes in the domain: $\mathbb{R}^{n}\backslash\Omega$; if the fluxes are zero in it, then one gets the so-called reflecting boundary conditions of the equations. Microscopically the motion of particles undergoing Lévy flights or tempered Lévy flights are much different from the Brownian motion; very rare but extremely long jumps dominate the dynamics, making the trajectories of the particles discontinuous. As shown in Figure 4, the particles may jump into, or jump out of, or even pass through the domain: $\Omega$. But the number of particles inside $\Omega$ is conservative, which can be easily verified by making the integration of (55) in the domain $\Omega$, i.e.,

[TABLE]

where ${\bf n}$ is the outward-pointing unit normal vector on the boundary. If ${\bf j}\,|_{\mathbb{R}^{n}\backslash\Omega}$=0, then for (40) $\Delta^{\frac{\beta}{2}}p({\bf X},t)=\nabla\cdot{\bf j}=0$ in $\mathbb{R}^{n}\backslash\Omega$. So, the Neumann type boundary conditions for (40), (44), (47), and (50) can be, heuristically, defined as

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

respectively. The corresponding reflecting boundary conditions are with $g({\bf X})\equiv 0$.

**Remark:  ** The Neumann type boundary conditions (62)-(65) derived in this section are independent of the choice of the flux ${\bf j}$, provided that it satisfies the condition (55).

4 Well-posedness and regularity of the fractional PDEs with generalized BCs

Here, we show the well-posedesses of the models discussed in the above sections, taking the models with the operator $\Delta^{\frac{\beta}{2}}$ as examples; the other ones can be similarly proved. For any real number $s\in{\mathbb{R}}$, we denote by $H^{s}({\mathbb{R}}^{n})$ the conventional Sobolev space of functions (see [1, 20]), equipped with the norm

[TABLE]

The notation $H^{s}(\Omega)$ denotes the space of functions on $\Omega$ that admit extensions to $H^{s}({\mathbb{R}}^{n})$, equipped with the quotient norm

[TABLE]

where the infimum extends over all possible $\widetilde{u}\in H^{s}({\mathbb{R}}^{n})$ such that $\widetilde{u}=u$ on $\Omega$ (in the sense of distributions). The dual space of $H^{s}(\Omega)$ will be denoted by $H^{s}(\Omega)^{\prime}$. The following inequality will be used below:

[TABLE]

Let $H^{s}_{0}(\Omega)$ be the subspace of $H^{s}({\mathbb{R}}^{n})$ consisting of functions which are zero in ${\mathbb{R}}^{n}\backslash\Omega$. It is isomorphic to the completion of $C^{\infty}_{0}(\Omega)$ in $H^{s}(\Omega)$. The dual space of $H^{s}_{0}(\Omega)$ will be denoted by $H^{-s}(\Omega)$.

For any Banach space $B$, the space $L^{2}(0,T;B)$ consists of functions $u:(0,T)\rightarrow B$ such that

[TABLE]

and $H^{1}(0,T;B)=\{u\in L^{2}(0,T;B):\partial_{t}u\in L^{2}(0,T;B)\}$; see [11].

4.1 Dirichlet problem

For any given $g\in{\mathbb{R}}\cup(L^{2}(0,T;H^{\frac{\beta}{2}}({\mathbb{R}}^{n}))\cap H^{1}(0,T;H^{-\frac{\beta}{2}}({\mathbb{R}}^{n})))\hookrightarrow C([0,T];L^{2}({\mathbb{R}}^{n})$, consider the time-dependent Dirichlet problem

[TABLE]

The weak formulation of (68) is to find $p=g+\phi$ such that

[TABLE]

and

[TABLE]

It is easy to see that $a(\phi,q):=\int_{{\mathbb{R}}^{n}}\Delta^{\frac{\beta}{4}}\phi\,\Delta^{\frac{\beta}{4}}q\,d{\bf X}$ is a coercive bilinear form on $H^{\frac{\beta}{2}}_{0}(\Omega)\times H^{\frac{\beta}{2}}_{0}(\Omega)$ (cf. [31, section 30.2]) and $\ell(q):=\int_{\Omega}(f+\Delta^{\frac{\beta}{2}}g-\partial_{t}g)q\,d{\bf X}$ is a continuous linear functional on $L^{2}(0,T;H^{\frac{\beta}{2}}_{0}(\Omega))$. Such a problem as (70) has a unique weak solution (cf. [31, Theorem 30.A]).

The weak solution actually depends only on the values of $g$ in ${\mathbb{R}}^{n}\backslash\Omega$, independent of the values of $g$ in $\Omega$. To see this, suppose that $g,\widetilde{g}\in{\mathbb{R}}\cup(L^{2}(0,T;H^{\frac{\beta}{2}}({\mathbb{R}}^{n}))\cap H^{1}(0,T;H^{-\frac{\beta}{2}}({\mathbb{R}}^{n})))\hookrightarrow C([0,T];L^{2}({\mathbb{R}}^{n}))$ are two functions such that $g=\widetilde{g}$ in ${\mathbb{R}}^{n}\backslash\Omega$, and $p$ and $\widetilde{p}$ are the weak solutions of

[TABLE]

respectively. Then the function $p-\widetilde{p}\in L^{2}(0,T;H^{\frac{\beta}{2}}_{0}(\Omega))\cap H^{1}(0,T;H^{-\frac{\beta}{2}}(\Omega))$ satisfies

[TABLE]

Substituting $q=p-\widetilde{p}$ into the equation above immediately yields $p-\widetilde{p}=0$ a.e. in ${\mathbb{R}}^{n}\times(0,T)$.

4.2 Neumann problem

Consider the Neumann problem

[TABLE]

Definition 4.1** (Weak solutions).**

The weak formulation of (73) is to find $p\in L^{2}(0,T;H^{\frac{\beta}{2}}({\mathbb{R}}^{n}))\cap C([0,T];L^{2}(\Omega))$ such that

[TABLE]

satisfying the following equation:

[TABLE]

Theorem 4.2** (Existence and uniqueness of weak solutions).**

*If $p_{0}\in L^{2}(\Omega)$, $f\in L^{2}(0,T;H^{\frac{\beta}{2}}(\Omega)^{\prime})$ and $g\in L^{2}(0,T;H^{\frac{\beta}{2}}({\mathbb{R}}^{n}\backslash\Omega)^{\prime})$, then there exists a unique weak solution of (73) in the sense of Definition 4.1. *

*Proof   * Let $t_{k}=k\tau$, $k=0,1,\dots,N$, be a partition of the time interval $[0,T]$, with step size $\tau=T/N$, and define

[TABLE]

Consider the time-discrete problem: for a given $p_{k-1}\in L^{2}({\mathbb{R}}^{n})$, find $p_{k}\in H^{\frac{\beta}{2}}({\mathbb{R}}^{n})$ such that the following equation holds:

[TABLE]

In view of (66), the left-hand side of the equation above is a coercive bilinear form on $H^{\frac{\beta}{2}}({\mathbb{R}}^{n})\times H^{\frac{\beta}{2}}({\mathbb{R}}^{n})$, while the right-hand side is a continuous linear functional on $H^{\frac{\beta}{2}}({\mathbb{R}}^{n})$. Consequently, the Lax–Milgram Lemma implies that there exists a unique solution $p_{k}\in H^{\frac{\beta}{2}}({\mathbb{R}}^{n})$ for (4.2).

Substituting $\displaystyle q=p_{k}$ into (4.2) yields

[TABLE]

Then, summing up the inequality above for $k=1,2,\dots,N$, we have

[TABLE]

By applying Grönwall’s inequality to the last estimate, there exists a positive constant $\tau_{0}$ such that when $\tau<\tau_{0}$ we have

[TABLE]

Since any $q\in H^{\frac{\beta}{2}}(\Omega)$ can be extended to $q\in H^{\frac{\beta}{2}}({\mathbb{R}}^{n})$ with $\|q\|_{H^{\frac{\beta}{2}}({\mathbb{R}}^{n})}\leq 2\|q\|_{H^{\frac{\beta}{2}}(\Omega)}$, choosing such a $q$ in (4.2) yields

[TABLE]

which implies (via duality)

[TABLE]

The last inequality and (4.2) can be combined and written as

[TABLE]

If we define the piecewise constant functions

[TABLE]

and the piecewise linear function

[TABLE]

then (4.2) and (4.2) imply

[TABLE]

and

[TABLE]

respectively, where the constant $C$ is independent of the step size $\tau$. The last inequality implies that $p^{(\tau)}$ is bounded in $H^{1}(0,T;H^{\frac{\beta}{2}}(\Omega)^{\prime})\cap L^{2}(0,T;H^{\frac{\beta}{2}}({\mathbb{R}}^{n}))\hookrightarrow C([0,T];L^{2}(\Omega))$. Consequently, there exists $p\in H^{1}(0,T;H^{\frac{\beta}{2}}(\Omega)^{\prime})\cap L^{2}(0,T;H^{\frac{\beta}{2}}({\mathbb{R}}^{n}))\hookrightarrow C([0,T];L^{2}(\Omega))$ and a subsequence $\tau_{j}\rightarrow 0$ such that

[TABLE]

By taking $\tau=\tau_{j}\rightarrow 0$ in (4.2), we obtain (75). This proves the existence of a weak solution $p$ satisfying (74).

If there are two weak solutions $p$ and $\widetilde{p}$, then their difference $\eta=p-\widetilde{p}$ satisfies the equation

[TABLE]

Substituting $q=p-\widetilde{p}$ into the equation yields

[TABLE]

which implies $p=\widetilde{p}$ a.e. in ${\mathbb{R}}^{n}\times(0,T)$. The uniqueness is proved.

**Remark:  ** From the analysis of this section we see that, although the initial data $p_{0}({\bf X})$ physically exists in the whole space ${\mathbb{R}}^{n}$, one only needs to know its values in $\Omega$ to solve the PDEs (under both Dirichlet and Neumann boundary conditions).

5 Conclusion

In the past decades, fractional PDEs become popular as the effective models of characterizing Lévy flights or tempered Lévy flights. This paper is trying to answer the question: What are the physically meaningful and mathematically reasonable boundary constraints for the models? We physically introduce the process of the derivation of the fractional PDEs based on the microscopic models describing Lévy flights or tempered Lévy flights, and demonstrate that from a physical point of view when solving the fractional PDEs in a bounded domain $\Omega$, the informations of the models in $\mathbb{R}^{n}\backslash\Omega$ should be involved. Inspired by the derivation process, we specify the Dirichlet type boundary constraint of the fractional PDEs as $p({\bf X},t)|_{\mathbb{R}^{n}\backslash\Omega}=g({\bf X},t)$ and Neumann type boundary constraints as, e.g., $(\Delta^{\beta/2}p({\bf X},t))|_{\mathbb{R}^{n}\backslash\Omega}=g({\bf X},t)$ for the fractional Laplacian operator.

The tempered fractional Laplacian operator $(\Delta+\lambda)^{\beta/2}$ is physically introduced and mathematically defined. For the four specific fractional PDEs given in this paper, we prove their well-posedness with the specified Dirichlet or Neumann type boundary constraints. In fact, it can be easily checked that these fractional PDEs are not well-posed if their boundary constraints are (locally) given in the traditional way; the potential reason is that locally dealing with the boundary contradicts with the principles that the Lévy or tempered Lévy flights follow.