fractional kinetic Fokker-Planck

# An operator splitting scheme for the fractional kinetic Fokker-Planck equation

Manh Hong Duong Department of Mathematics, Imperial College London, London SW7 2AZ, UK  and  Yulong Lu Department of Mathematics, Duke University, Durham NC 27708, USA
###### Abstract.

In this paper, we develop an operator splitting scheme for the fractional kinetic Fokker-Planck equation (FKFPE). The scheme consists of two phases: a fractional diffusion phase and a kinetic transport phase. The first phase is solved exactly using the convolution operator while the second one is solved approximately using a variational scheme that minimizes an energy functional with respect to a certain Kantorovich optimal transport cost functional. We prove the convergence of the scheme to a weak solution to FKFPE. As a by-product of our analysis, we also establish a variational formulation for a kinetic transport equation that is relevant in the second phase. Finally, we discuss some extensions of our analysis to more complex systems.

###### Key words and phrases:
Operator splitting methods, variational methods, fractional kinetic Fokker-Planck equation, kinetic transport equation, optimal transportation.
###### 2010 Mathematics Subject Classification:
Primary: 49S05, 35Q84; Secondary: 49J40.

## 1. Introduction

In this paper, we study the existence of solutions to the following fractional kinetic Fokker-Planck equation (FKFPE)

 (1.1) {∂tf+v⋅∇xf=divv(∇Ψ(v)f)−(−△v)sfinRd×Rd×(0,∞),f(x,v,0)=f0(x,v)inRd×Rd,

with . In the above, denotes the divergence operator; the differential operators and with subscripts and indicate that these operators act only on the corresponding variables; the operator is the fractional Laplacian operator on the variable , where the fractional Laplacian , is defined by

 −(−△)sf(x):=−F−1(|ξ|2sF[f](ξ))(x).

Here denotes the Fourier transform on , i.e. . Note that the fractional Laplacian operator with is a non-local operator since it can also be expressed as the singular integral

 −(−△)sf(x)=−Cd,s∫Rdf(x)−f(y)|x−y|d+2sdy,

where the normalisation constant is given by and is the Gamma function. See [33] for more equivalent definitions of fractional Laplacian operator.

The equation (1.1) is interesting to us because it can be viewed as the Fokker-Planck (forward Kolmogorov) equation of the following generalized Langevin equation

 (1.2) dXtdt=Vt, dVtdt=−∇Ψ(Vt)+Lst,

where is the Lévy stable process with exponent . The stochastic differential equation (SDE) (1.2) describes the motion of a particle moving under the influence of a (generalized) frictional force and a stochastic noise and in the absence of an external force field. FKFPE (1.1) is the evolution of the probability distribution of . In particular, the fractional operator is the Markov generator of the process . When and , equation (1.1) becomes the classical kinetic Fokker-Planck (or Kramers) equation (without external force field) which is a local PDE and has been used widely in chemistry as a simplified model for chemical reactions  [32, 26] and in statistical mechanics [35, 38]. The non-local Lévy process plays an important role in modelling systems that include jumps and long-distance interactions such as anomalous diffusion or transport in confined plasma [5]. Singular limits of Equation (1.1) with was studied in [12], see also [11] for a similar result for the same equation but on a spatially bounded domain. In a recent work  [1], the authors have extended [12] to a system that contains an additional external force field and they have also proved its well-posedness by the means of the Lax-Milgram theorem. We will prove the existence of solutions of (1.1) for a general based on the trick of operator splitting. For more recent developments on PDEs involving the fractional Laplacian operator, we refer the interested reader to expository surveys [42, 43, 41].

The aim of this paper is to develop a variational formulation for approximating solutions to equation (1.1). The theory of variational formulation for PDEs took off with the introduction of Wasserstein gradient flows by the seminal work of Jordan, Kinderlehrer and Otto [30]. Such a variational structure has important applications for the analysis of an evolution equation such as providing general methods for proving well-posedness [4] and characterizing large time behaviour (e.g., [10]), giving rise to natural numerical discretizations (e.g., [22]), and offering techniques for the analysis of singular limits (e.g., [39, 40, 6, 18]). There are now a significantly large number of papers in exploring variational structures for local PDEs, see the aforementioned papers and references therein as well as the monographs [4, 44] for more details. However, variational formulations for non-local PDEs are less understood. Erbar [24] showed that the fractional heat equation is a gradient flow of the Boltzmann entropy with respect to a new modified Wasserstein distance that is built from the Lévy measure and based on the Benamou-Brenier variant of the Wasserstein distance. Bowles and Agueh  [8] proved the existence of the fractional Fokker-Planck equation

 (1.3) {∂tf=divv(∇Ψ(v)f)−(−△v)sfinRd×(0,∞),f(v,0)=f0(v)inRd,

which can be viewed as the spatially homogeneous version of equation (1.1) or the fractional heat equation with a drift. Erbar’s proof is variational based on the so-called “evolution variational inequality” concept introduced in [4]. However, it seems that his method can not be extended to the fractional Fokker-Planck equation since the distance that he introduced was particularly tailored for the Boltzmann entropy. Instead, Bowles and Agueh’s proof is “semi-variational” based on a novel splitting argument which we sketch now. They split up the original dynamics (1.3) into two processes: a fractional diffusion process, namely , and a transport process in the field of the potential , namely , and then alternatively run these processes on a small time interval. Furthermore, the transport process can be understood as a Wasserstein gradient flow of the potential energy. By adopting a suitable interpolation of the individual processes, they were able to show that the constructed splitting scheme converges to a weak solution of (1.3). In the literature, the technique of operator splitting is often used to construct numerical methods for solving PDEs, see [27]. On the theoretical side, the idea of splitting had also been used to study the well-posedness of PDEs, see [9, 2] on kinetic equations and [3, 16] on fractional PDEs.

In the present work, we adopt the same splitting argument in [8] to construct a weak solution to the fractional kinetic equation (1.1). More specifically, we split the dynamics described in (1.1) by two phases:

1. Fractional diffusion phase. At every fixed position , the probability density , as a function of velocity , evolves according to the fractional heat equation

 (1.4) ∂tf=−(−△v)sf.
2. Kinetic transport phase. The density evolves according to the following equation

 (1.5) ∂tf+v⋅∇xf=divv(∇Ψ(v)f).

We expect that successive alternative iterating the above two phases with vanishing period of time would give an approximation to the dynamics (1.1). The key difference between our splitting scheme above and the scheme in [8] is that the transport process here is not only driven by the potential energy but also the kinetic energy. In [8], the transport process is approximated by a discrete Wasserstein gradient flow based on the work [31]. However, due to the presence of the kinetic term, the kinetic transport equation is not a Wasserstein gradient flow; thus one can no longer use the Wasserstein distance. To overcome this obstacle, we employ instead the minimal acceleration cost function and the associated Kantorovich optimal transportation cost functional that has been used in [28, 19] for the kinetic Fokker-Planck equation and in [25] for the isentropic Euler system, see Section 3.

### 1.1. Main result

Throughout the paper, we make the following important assumption on the potential .

###### Assumptions 1.1.

is non-negative and .

We adopt the following notion of weak solution to KFPE (1.1).

###### Definition 1.2.

Let be a non-negative function such that for some and . We say that is a weak solution to (1.1) if it satisfies the following:

1. for any .

2. for a.e. .

3. For any test function ,

 ∫T0∫R2df(x,v,t)(∂tφ+v⋅∂xφ+∇vΨ⋅∇vφ−(−△v)sφ)dtdxdv +∫R2df0(x,v)φ(0,x,v)=0.

The main result of the paper is the following theorem.

###### Theorem 1.3.

Suppose that Assumption 1.1 holds. Given a for some and , there exists a weak solution to (1.1) in the sense of Definition 1.2.

The proof of Theorem 1.3 is constructive, that is we will build a converging sequence to a solution of (1.1) from the splitting scheme discussed above that will be rigorously formulated in Section 4. The proof is based on a series of lemmas and is postponed to Section 5. As a by-product of the analysis, we also construct a discrete variational scheme and obtain its convergence for the kinetic transport equation, see Theorem 3.3 in Section 3; thus extending the work [31] to include the kinetic feature. Furthermore, some possible extensions to more complex systems are discussed in Section 6. It is not clear to us how to obtain the uniqueness and regularity result. The bootstrap argument in [30] to prove smoothness of weak solutions (and hence also uniqueness) seems not working for the fractional Laplacian operator due to the lack of a product rule. It should be mentioned that in the recent paper [34], the author has proved the existence and uniqueness of a solution to the fractional Fokker-Planck equation (1.3) in some weighted Lebesgue spaces. It would be an interesting problem to generalize [34] to FKFPE. This is to be investigated in future work.

### 1.2. Organization of the paper

The rest of the paper is organized as follows. Section 2 summarizes some basic results about the fractional heat equation. Section 3 studies the kinetic transport equation and its variational formulation. The splitting scheme of the paper is formulated explicitly in Section 4 and some a priori estimates are established for the discrete sequences as well as their time-interpolation. The proof of the main result is presented in Section 5. Finally, in Section 6 we discuss several possible extensions of the analysis to more complex systems.

### 1.3. Notation

Let be the collection of probability measures on with finite second moments. Let be the subset of probability measures in that are absolutely continuous with respect to the Lebesgue measure on . For , the 2-Wasserstein distance is defined by

 W2(μ,ν):=(inf{∫R2d|x−y|2p(dx,dy):p∈P(μ,ν)})12

where is the set of probability measures on with marginals , i.e. if and only if

 p(A×Rd)=μ(A),p(Rd×A)=ν(A)

hold every Borel set . In the case that with densities , we may write instead of .

We use the notation to denote the push-forward of a probability measure on under map , that is a probability measure on satisfying for all smooth test function ,

 ∫R2dφ(x,v)dF#μ=∫R2dφ(F(x,v))dμ.

## 2. The fractional heat equation

This section collects some basic results on the fractional heat equation. We start by defining the fractional heat kernel

 (2.1) Φs(v,t):=F−1(e−t|⋅|2s)(v).

Remember that the fractional Laplacian operator in (1.1) is only an operator in -variable. With the fractional heat kernel, the solution to the fractional heat equation (1.4) with initial condition can be expressed as

 (2.2) f(x,v,t)=Φs(⋅,t)∗vf0(x,v)

where is the convolution operator in -variable. The following elementary result is immediate from the definition of the kernel; see also [8].

###### Lemma 2.1.

For any , .

For any and , .

for all and .

Lemma 2.1 (3) demonstrates a significant difference between the fractional heat kernel and standard Gaussian kernel, i.e. the former has infinite second moment. The loss of second moment bound may lead to infinite potential energy for example when the potential . To overcome this issue, it is more convenient to make a renormalisation on the fractional heat kernel. To be more precise, for any , let us denote and set where is an indicator function of a centred ball of radius . Given a function , we can define the renormalised convolution

 (2.3) ¯fh,R:=Φhs,R∗vf∥Φhs,R∥L1(Rd)

It is clear that the new defined convolution satisfies pointwise. Moreover, we have the following lemma.

###### Lemma 2.2.

Let be a function such that . Suppose that and with . Then

.

 (2.4) ∫R2d¯fh,R(x,v)F(v)dxdv ≤∫R2df(x,v)F(v)dxdv +12∥D2F∥L∞(Rd)∫BR|w|2Φhs(w)dw∫BRΦhs(w)dw.
###### Proof.

Notice that it suffices to prove part (2) since part (1) follows directly from part (2) by setting . The proof is similar to that of [8, Lemma 4.1], but for completeness we give the proof below. First from the definition of , one sees that

 ∫R2d¯fh,R(x,v)F(v)dxdv=∫R2dF(v)∫BRΦhs(w)f(x,v−w)dwdxdv∫BRΦhs(w)dw.

Using change of variable and Taylor’s expansion, we can write the numerator as

 ∫R2dF(v)∫BRΦhs(w)f(x,v−w)dwdxdv =∫R2dF(w+z)∫BRΦhs(w)f(x,z)dwdxdz =∫BRΦhs(w)∫R2dF(w+z)f(x,z)dxdzdw =∫BRΦhs(w)∫R2d[F(z)+w⋅∇F(z)+12wTD2F(ξw,z)w]f(x,z)dxdzdw ≤∫BRΦhs(w)dw∫R2dF(z)f(x,z)dxdz+∣∣∣∫BRΦhs(w)∫R2dw⋅∇F(z)f(x,z)dxdzdw∣∣∣ +12∥D2F∥∞∫BR|w|2Φhs(w)∫R2df(x,z)dxdz =∫BRΦhs(w)dw∫R2dF(z)f(x,z)dxdz+12∥D2F∥∞∫BR|w|2Φhs(w).

Note that in the above is an intermediate point between and and the term with the modulus vanishes since the kernel is symmetric with respect to the origin.

The following lemma provides an upper bound for the ratio on the right side of (2.4).

###### Lemma 2.3.

For any , there exists a constant such that

 (2.5) ∫BR|w|2Φhs(w)dw∫BRΦhs(w)dw≤C(h1s+hR2−2s)

holds for all .

###### Proof.

This lemma follows directly from a two-sided point-wise estimate on as shown in [8, Proposition 2.1]. See also equation (16) in [8]. ∎

## 3. The kinetic transport equation and its variational formulation

### 3.1. The minimum acceleration cost

Consider the kinetic transport equation with initial value

 (3.1) ∂tf(x,v,t)+v⋅∇xf=divv(∇Ψ(v)f(x,v,t)), f(x,v,0)=f0(x,v).

We are interested in the variational structure of (3.1) which is an interesting problem on its own right. In [31], Kinderlehrer and Tudorascu proved that the transport equation , which is the spatially homogeneous version of (3.1), is a Wasserstein gradient flow of the energy . Their proof is via constructing a discrete variational scheme as in [30]. However, due to the absence of the entropy term, which is super-linear, several non-trivial technicalities were introduced to obtain the compactness of the discrete approximations thus establishing the convergence of the scheme. For the kinetic transport equation (3.1), due to the presence of the kinetic term, it is not a Wasserstein gradient flow in the phase space thus the Wasserstein distance can no longer be used. Therefore to construct a discrete variational scheme for this equation, we need a different Kantorovich optimal transportation cost functional. To this end, we will employ the Kantorovich optimal transportation cost functional that is associated to the minimal acceleration cost. This cost functional has been used before in [28, 19] for the kinetic Fokker-Planck equation and in [25] for the isentropic Euler system. We follow the heuristics of defining the minimal acceleration cost as in [25]. Consider the motion of particle going from position with velocity to a new position with velocity , within a time interval of length . Suppose that the particle follows a curve such that

 (ξ,˙ξ)|t=0=(x,v)  and  (ξ,˙ξ)|t=h=(x′,v′)

and such that the average acceleration cost along the curve, that is is minimized. Then the curve is actually a cubic polynomial and the minimal average acceleration cost is given by where

 (3.2) Ch(x,v;x′,v′):=|v′−v|2+12∣∣x′−xh−v′+v2∣∣.

The Kantorovich functional associated with the cost function , is defined by, for any ,

 (3.3) Wh(μ,ν)2=infp∈P(μ,ν)∫R4dCh(x,v;x′,v′)p(dxdvdx′dv′),

where is the set of all couplings between and . It is important to notice that is not a distance. In fact, is not symmetric in the arguments , due to the asymmetry of the cost function . In addition, does not vanish when . Instead, we have that

 Wh(μ,ν)=0⟺ν=(Fh)#μ,

where is the free transport map defined by

 Fh:Rd×Rd →Rd×Rd (3.4) (x,v) ↦Fh(x,v)=(x+hv,v).

It is also useful to define the map

 Gh:Rd×Rd →Rd×Rd (3.5) (x,v) ↦Gh(x,v)=(√3(2xh−v),v).

The composition is then given by

 (Gh∘Fh)(x,v)=(√3(2xh+v),v).

Although the Kantorovich functional is not a distance, the next lemma shows that can be expressed in terms of the usual Wasserstein distance .

###### Lemma 3.1.

[25, Proposition 4.4] Let and be given by (3.1) and (3.1) respectively. The Kantorovich functional can be expressed in terms of the 2-Wasserstein distance as

 (3.6) Wh(μ,ν)=W2((Gh∘Fh)#μ,G#hν)for all  μ,ν∈P2(R2d).

As a consequence, the infimum in (3.3) is attained and thus is a minimum.

### 3.2. Variational formulation

With being defined, we want to interpret (3.1) as a generalized gradient flow of the potential energy with respect to . For doing so, we consider the variational problem

 (3.7) inff∈P2aA(f):=12hWh(f0,f)2+∫R2dΨ(v)f(x,v)dxdv.

Here is an initial probability density with and is the time step. The next lemma establishes some properties about the minimizer to (3.7).

###### Lemma 3.2.

For being sufficiently small, the variational problem (3.7) has a unique minimizer .

Let be small enough such that for some fixed . If for , then

 (3.8) ∥f∥pLp(R2d)≤(1−αh)p−1∥f0∥pLp(R2d).

satisfies the following Euler-Lagrange equation: for any ,

 (3.9) 1h∫R4d[(x′−x)⋅∇x′φ(x′,v′)+(v′−v)⋅∇v′φ(x′,v′)]P∗(dxdvdx′dv′) −∫R2dv′⋅∇x′φ(x′,v′)f(x′,v′)dx′dv′ +∫R2d∇v′Ψ(v′)⋅∇v′φ(x′,v′)f(x′,v′)dx′dv′=R,

where is the optimal coupling in and

 (3.10) R =−h2∫R4d∇v′Ψ(v′)⋅∇x′φ(x′,v′)P∗(dxdvdx′dv′) =−h2∫R2d∇v′Ψ(v′)⋅∇x′φ(x′,v′)f(x′,v′)dx′dv′.
###### Proof.

Thanks to Lemma 3.1, we can rewrite the functional as

 A(f) =12hW2((Gh∘Fh)#f0,(Gh)#f)2+∫R2dΨ(v)(Gh)#f(dxdv) =12hW2(~f0,~f)2+∫R2dΨ(v)~f(dxdv)=:~A(~f),

where and . According to [31, Proposition 1] (see also [8, Proposition 3.1 ]), the functional has a unique minimizer, denoted by . Therefore, the problem (3.7) has a unique minimizer .

This follows directly from [31, Proposition 1] and the fact that if then

 ∥~f∥pLp(R2d)=(2h√3)d(p−1)∥f∥pLp(R2d).

The derivation of the Euler-Langrange equation for the minimizer of the variational problem (3.7) follows the now well-established procedure (see e.g. [30, 28]). For the reader’s convenience, we sketch the main steps here. First, we consider the perturbation of defined by push-forwarding under the flows :

 ∂ψs∂s=ζ(ψs,ϕs), ∂ϕs∂s=η(ψs,ϕs), ψ0(x,v)=x, ϕ0(x,v)=v,

where will be chosen later. Let us denote to be the push forward of under the flow . Since , it follows that , and an explicit calculation gives

 (3.11) ∂sγs∣∣s=0=−divx(fζ)−divv(fη)

in the sense of distributions. Second, thanks to the optimality of , we have that for all defined via the flow above. Then the standard variational arguments as in [30, 28] leads to the following stationary equation on :

 12h∫R4d[∇x′Ch(x,v;x′,v′)⋅ζ(x′,v′)+∇v′Ch(x,v;x′,v′)⋅η(x′,v′)]P∗(dxdvdx′dv′) (3.12) +∫R2df(x,v)∇Ψ(v)⋅η(x,v)dxdv=0,

where is the optimal coupling in the definition of . Third, we choose and with a given as follows

 (3.13) ζ(x′,v′) =−h26∇x′φ(x′,v′)+12h∇v′φ(x′,v′), η(x′,v′) =−12h∇x′φ(x′,v′)+∇v′φ(x′,v′).

Now from the definition of the cost functional in (3.2), we have that

 ∇x′Ch =24h(x′−xh−v′+v2), ∇v′Ch =2(v′−v)−12(x′−xh−v′+v2).

Therefore, together with (3.13), we calculate

 ∇x′Ch⋅ζ+∇v′Ch⋅η =24h(x′−xh−v′+v2)⋅(−h26∇x′φ(x′,v′)+12h∇v′φ(x′,v′)) +(2(v′−v)−12(x′−xh−v′+v2))⋅(−12h∇x′φ(x′,v′)+∇v′φ(x′,v′)) =2((x′−x)−hv′)⋅∇x′φ+2(v′−v)⋅∇v′φ.

The Euler-Lagrange equation (3.9) for the minimizer follows directly by substituting the equation above back into (3.12). ∎

We now can build up a discrete variational scheme for the kinetic transport equation as follows. Given with and is the time step. For every integer , we define as the minimizer of the minimization problem

 (3.14) inff∈P2a{12hWh(fk−1,f)2+∫R2dΨ(v)f(x,v)dxdv}.

The following theorem extends the work [31] to the kinetic transport equation.

###### Theorem 3.3.

Suppose that Assumption 1.1 holds. Given a for some and , there exists a weak solution to equation  (3.1) in the sense of Definition 1.2 but with the fractional Laplacian term removed.

###### Proof.

The proof of this theorem follows the same lines as that of the Theorem 1.3, that is to show that the discrete variational scheme (3.14) above converges to a weak solution of the kinetic transport equation. Since the proof of Theorem 1.3 will be carried out in details in Section 5, we omit this proof here. ∎

## 4. A splitting scheme for FKFPE

### 4.1. Definition of splitting scheme

As we mentioned in the introduction section, we object to construct an operator splitting scheme for equation (1.1) by continuously alternating processes (1.4) and (1.5), where the later is approximated by the generalized gradient flow of the potential energy, or equivalently, the density after a short time step is approximately given by the solution to the variational problem (3.7). However, there is an issue associated with iterating (1.4) and (3.7). That is, the solution of the fractional heat equation may not have a finite second moment (see Lemma 2.1 (3)). Hence it can not be used as the initial condition in the variational problem (3.7) since the potential energy might be infinite. To around this issue, we define an approximate fractional diffusion process by using the renormalised convolution (2.3) based on the truncted fractional heat kernel. To be more precise, given a fixed , let us consider a uniform partition of the time interval with and . With an initial condition , for we iteratively compute the following:

• Given a trunction parameter , compute the renormalised convolution

 (4.1) ¯fnh,R:=Φhs,R∗vfn−1h,R∥Φhs,R∥L1(Rd).
• Solve for the minimizer of the problem

 (4.2) fn+1h,R:=argminf∈P2a(Rd)12hWh(¯fnh,R,f)2+∫R2dΨ(v)f(x,v)dxdv.

Note that thanks to Lemma 3.2 (1) the minimizer in (4.2) is well-defined and unique. Moreover, it follows from Lemma 3.2 (3) that satisfies the following equation

 (4.3) ∫R4d[(x′−x)⋅∇x′φ(x′,v′)+(v′−v)⋅∇v′φ(x′,v′)]Pn+1h,R(dxdvdx′dv′) =h∫R2d(v′⋅∇x′φ(x′,v′)−∇v′Ψ(v′)⋅∇v′φ(x′,v′))fn+1h,R(x′,v′)dx′dv′+Rn+1h,R,

where is the optimal coupling in and

 (4.4) Rn+1h,R=h22∫R2d∇vΨ(v)⋅∇xφ(x,v)fn+1h,R(dxdv).

With the scheme being defined above, we obtain a discrete approximating sequence . Below we define a time-interpolation based on and our ultimate goal is to prove that this sequence converges to a weak solution of (1.1).

Time-interpolation: We define by setting

 (4.5) fh,R(t):=Φs(t−tn)∗vfnh,R for t∈[tn,tn+1).

It is clear that by definition solves the fractional heat equation on every with initial condition . Notice also that is only right-continuous in general. For convenience, we also define

 (4.6) ~fn+1h,R=limt↑tn+1fh,R(t).

### 4.2. A priori estimates

In this section, we prove some useful a priori estimates for the discrete-time sequence as well as for the time-interpolation sequence . We start by proving an upper bound for the sum of the Kantorovich functionals .

###### Lemma 4.1.

Let and be the sequences constructed from the splitting scheme. Then there exists a constant , independent of and , such that

 (4.7) N∑n=1Wh(¯¯¯fnh,R,fnh,R)2≤C(h∫R2dΨ(v)f0(x,v)dxdv+T∥D2Ψ∥∞(h1/2+hR2−2s)).

Since