Method of variable separation for investigating exact solutions and dynamical properties of the time-fractional Fokker-Planck equation ∗†

Based on the idea of variable separation, the time-fractional Fokker-Planck equation with external force ﬁeld is studied by using the property of Mittag-Leﬄer function and some special algorithm skills. In the cases of various external potential functions such as linear potential, harmonic potential, logarithmic potential, exponential potential, and quartic potential, exact solutions and dynamical properties of the above mentioned equation is investigated. The some interesting dynamical behaviors and phenomena are discovered. The proﬁles of some representative exact solutions are il-lustrated by 3D-graphs.


Introduction
Since the original concept of fractional-order derivative appeared in a famous correspondence which L'Hôspital sent to Leibniz in 1695, many mathematicians have further developed this concept. Through more than 300 years of development, the research works on theory and applications of the fractional calculus and the fractional differential equations have become many hot topics in the scientific field. Especially in recent decades, fractional calculus theory and modeling methods have been well applied in the fields of high energy physics, statistical physics, anomalous diffusion, viscoelastic material mechanics, complex control, signal processing, biomedical engineering, geophysics, rheology and economics and so forth in other research areas, and their unique advantages can not be replaced by integer-order calculus and differential equations.
In many scientific research fields mentioned above, the research work on anomalous diffusion phenomenon is undoubtedly a very important research topic. Anomalous diffusion is a kind of common phenomenon in nature, it is present in many physical situations such as transport of fluid in porous media, surface growth, diffusion of plasma, diffusion on liquid surface, two-dimensional rotating flow and so on, the details of introduction can be seen in the references [1][2][3][4] and references cited therein. As an anomalous diffusion model in an external potential field, the following time-fractional Fokker-Planck equation [5][6][7][8][9][10][11] ∂P (x, t) ∂t has a wide range of applications in various scientific are as such as biology, chemistry, physics and finance. As a fractional subdiffusion model of continuous-time random walk, the equation (1.1) is also rewritten as , (1.2) where the function P (x, t) defines the probability density of finding the test particle at a certain position x and a given time t with (t, x) ∈ [0, ∞)×(−∞, +∞), the R 0 D 1−α t is fractional differential operator of Riemann-Liouville type and 0 < α < 1, the function V (x) indicates the potential of over-damped Brownian motion and V ′ (x) = dV (x) dx , m denotes the mass of the test particle, η α is a constant which denotes the generalized friction coefficient with dimension then the equation (1.1) can be easily reduced to the equation (1.2). For these details, please participate in the reference [8]. Especially, when α → 1, equation (1.1) becomes a classical Fokker-Planck equation [12][13][14] (a model in statistical physics) as follows: The equation (1.3) is an normal diffusion model and the corresponding diffusion process is governed by Fick's second law in the force-free limit V (x) = c (constant).
By employing the fractional differential operator R 0 D α−1 t to act on both sides of the By using the following relation of the definitions of Riemann-Liouville fractional derivative and Caputo fractional derivative the equation (1.4) can be reduced to a simple form in the definition of Caputo fractional operator as follows: For details of the above reduction process, see the reference [8,15]. From the above reduction process, the equation ( its analytic solutions in closed form in terms of the Fox H−function were obtained by Wyss [16], Hilfer [17], Metzler and Klafter [8,18], respectively. In the harmonic potential (Ornstein- which studied by Metzler et al [5,8] and its series solution was obtained in [5]. In the quartic double-well potential V (x) = ax 4 +bx 2 and the sextic double-well potential V (x) = ax 6 +bx 4 , by using the method of numerical calculation, the authors investigated the propagators, the auto-correlation function, the survival probality and the lifetime distribution of the equation  [20,21]. Under magnetic field, the study of plasma transport phenomena (diffusion of plasma) can be modeled by classical Fokker-Planck equation on infinite interval [22]. Under the action of various external force fields, the particle motion described by the fractional an exact solution. However, it is very difficult to obtain the exact solution of fractional partial differential equations. Many effective methods in the field of integer order can not be directly applied to solve partial differential equations in the field of fractional order. It is often necessary for us to design a suitable method for solving a specific equation. It is worth mentioning that with the development of the research, on the solving methods of fractional differential equations, people have developed some effective methods to search the exact solutions and the approximate analytical solutions of the fractional differential equation in recent years, these methods contain Adomian decomposition method [23,24], homotopy analysis method [25,26], first integral method [27], invariant analysis method [28,29], fractional variational iteration method [30][31][32], invariant subspace method [33][34][35][36], homogenous balanced principle combined with integral bifurcation method [37][38][39][40], method of separating variables [41][42][43], and so forth. Similar to the methods in references [37][38][39][40][41][42][43], in this paper, by using method of separating variables together with Mittag-Leffler function, we will investigate exact solutions and dynamical properties of the time-fractional Fokker-Planck equation (1.5) in various potential functions such as linear potential, harmonic potential, logarithmic potential, exponential potential and quartic potential.
The organization of this paper is as follows: By using method of separating variables together with Mittag-Leffler function, in various external potential function V (x), we will investigate different kinds of exact solutions of the time-fractional Fokker-Planck equation (1.5) defined by Caputo differential operator. Further, we will discuss the dynamical properties of these exact solutions.

Exact solutions and dynamical properties of the fractional model (1.5) in various external potential fields
In this section, based on the idea of separating variables, by using the property of Mittag-Leffler function and some algorithm skills, we will investigate exact solutions of the equation (1.5) and discuss their dynamical properties. The equation (1.5) can be rewritten as the Obviously, the traditional method is to make Laplace transform for t and Fourier transform for x on both sides of equation (2.1), and then use the corresponding inverse transformation to obtain the solution of equation (2.1). Indeed, anyone who studied fractional differential equations knows that the expression of the exact solution of the equation (2.1) can not be obtained by the two transform because the number of terms of the equation is too many (more than two terms). However, we can use a special variable separation method to obtain the exact solution of the equation. According to the idea of separating variables, similar to the methods in reference [37][38][39][40][41][42][43], we suppose that equation (2.1) has solutions as form in the below where W (x) is an undetermined function which can be determined later. By substituting In equation (2.3), dividing out the Mittag-Leffler function E α (λt α ), we obtain an ODE as follows: By the way, by using the classical variable separation method appeared in [9] and reference cited therein, the form of solutions of equation (2.1) can be supposed as for a given eigenvalue λ n,α , where n ∈ N. Under this hypothesis, equation (2.1) can be reduced the following two eigenequations It is easy to obtain the solution of equation (2.6) as follows: Thus, the solution of equation (2.1) can be expressed as Next, by using equation (2.4) we will search exact solutions of the equation (1.5) and discuss their dynamical properties in various external potential fields.

Case of the linear potential field
When V = mω 2 x is a linear potential function, equation (2.4) can be reduced to a linear and homogeneous ODE of constant coefficient as follows: , then equation (2.10) has a general solution, where ∆ = ω 4 + 4λη 2 α K α , C 1 and C 2 are two arbitrary constants. Plugging (2.11) into (2.2), we obtain a kind of exact solution of the equation (1.5) as follows: In the range of − ω 4 4η 2 α Kα < λ < 0, the solution (2.12) has attenuation property because E α (λt α ) → 0 as t → +∞.
α Kα , then equation (2.10) has a general solution, where C 1 , C 2 are two arbitrary constants. Plugging (2.13) and λ = − ω 4 4η 2 α Kα into (2.2), we obtain a type of exact solution of the equation (1.5) as follows: . (2.14) If λ < − ω 4 4η 2 α Kα , then equation (2.10) has a general solution, Here, we naturally have to ask which of the above three solutions is the unique solution of and the boundary condition is given by then by using (2.16) we obtain an unique solution corresponding to the mixed problem, Similarly, if the initial condition of equation (1.5) is given by and the boundary condition is given by then by using (2.14) we obtain an unique solution corresponding to above mixed problem, It is easy to know that the solutions (2.17) and (2.18) have attenuating property (decay characteristic) according to time increase, which satisfy P (x, t) → 0 as time t → +∞.
In order to intuitively show the dynamical profiles (or properties) of above solutions obtained by us, as examples, the 3D-graphs and 2D-graphs of the solutions (2.17) and (2.18) are illustrated, which are shown in Fig. 1(a), (b), (c) and (d), respectively. phenomenon of plasma transport in a infinite region; this shows that the probability of plasma acoustic particles being distributed at infinity is very low, that is to say, the more distant the particle acoustic wave is transmitted, the weaker the signal is.
In the following discussion, we only discuss the solution of equation (1.5) (i.e. equation (2.1)) under free conditions, and no longer discuss the special solution under a certain boundary condition and an initial value condition in detail. Of course, as long as the certain boundary condition and the initial value condition are given, the special solution of the practical problem can always be obtained by the solution under the free condition.

Case of the harmonic potential field
When V = 1 2 mω 2 x 2 is a harmonic potential function, equation (2.4) can be reduced to a linear and homogeneous ODE of variable coefficient as follows: Especially, when λ = ω 2 ηα , equation (2.19) becomes the following simple equation it has a general solution formed as follows: defines an error function and C 1 , C 2 are two arbitrary constants. Plugging (2.21) and λ = ω 2 ηα into (2.2), we obtain a kind of solution of equation (1.5) formed as π is a new arbitrary constant, the (2.23) is a solution of infinite series type. The 3D and 2D graphs of the profiles for the solution (2.22) are shown in Fig.2(a) and (b), respectively.
It can be seen from the Fig.2 that the probability density function defined by the solution (2.22) does not satisfy the case of normal distribution. In the case of subdiffusion, this is a new type of distribution. Since the solution (2.22) is expressed by the error function, we call this distribution an error-type distribution. This phenomenon seems to be used to explain that the shape of sound wave of the plasma particle is kink-type. But what is being explained in other models is not known at the moment. When λ ̸ = ω 2 ηα , we make a transformation as follows: (2.24) By substituting (2.24) into (2.19), it can be reduced to Letting y(0) = C 1 and dy(0) dx = C 2 as initial conditions of the equation (2.25), then it is easily to obtain a 0 = C 1 , a 1 = C 2 in (2.27). After substituting (2.27) into (2.25), comparing with the same power coefficient of x in the reduced equation, we obtain an iterative formula for coefficients as follows: + λ Kα (n + 2)(n + 1) a n , (n = 0, 1, · · ·). (2.28) According to (2.28) and (2.27), we obtain a general solution of the series type of equation (2.25) as follows: ] . (2.29) Applying (2.29), (2.24) and (2.2), we get an analytic solution of equation (1.5) formed as Of course, directly solving the equation (2.26) and then changing the variables ξ = ω √ 2ηαKα x and ν = − ληα ω 2 back, the following result is obtained.
where C 1 , C 2 are two arbitrary constants, the function M (µ, ν, x) is the first kind of Kummer function which defined by and the Pochhammer symbol denotes (µ) n = Γ(µ+n) Γ(µ) . In the program editing of mathematical software Maple, the function M (µ, ν, x) is usually represented by notation KummerM (µ, ν, x).
According to (2.31), (2.24) and (2.2), we also obtain an analytic solution of equation (1.5) as follows: Although the solutions (2.30) and (2.32) are different in forms, they are in fact equivalent. As in Fig.1 and Fig.2, when the parametric values are properly selected, we plot the 3D and 2D graphs of the solution (2.32), which are shown in Fig.3 under λ < 0 and λ > 0, respectively.
Obviously, no matter how the value of the parameter λ changes, the probability density determined by the solution (2.32) always accords with the normal distribution, as can be seen from Fig. 3.
(a) 3D-graph as λ < 0 (b) 2D-graph at t = 2 as λ < 0 (c) 3D-graph as λ > 0 (d) 2D-graph at t = 2 as λ > 0 From the Fig.3, we know that the random variable defined by the solution (2.32) obeys normal distribution. The corresponding particle motion belongs to Brownian motion. This indicates that in the case of subdiffusion, if the external field is adjusted properly, the particle motion also satisfies Brownian motion and corresponding random variable also obeys normal distribution. It is well known that this phenomenon will not occur in the simple field of subdiffusion. However, with the interference of the external field, that is, the interaction between the subdiffusion and the external field, everything will become possible.
This phenomenon of normal distribution which appeared in the subdiffusion environment is the result of the interaction between the subdiffusion and the external field to achieve a new equilibrium.

Case of the exponential potential field
Kα x is an exponential potential function and λ > 0, equation (2.4) can be reduced to we get a general solution of (2.33) as follows: Kα x is an exponential function and λ > 0, equation (2.4) can be reduced to Similarly, we can obtain an exact solution of the equation (1.5) formed as In order to show profiles of the solutions (2.36) and (2.38) intuitively, we plot their graphs which are shown in Fig.4; under the values of the parameters λ = 0.2, K α = 0.9, α = 0.6, C 1 = 1, C 2 = −1, we plot 3D and 2D graphs of the solution (2.36) which are shown in Fig.4 (a) and (b). Under the values of the parameters λ = 0.2, K α = 0.6, α = 0.75, C 1 = C 2 = 0.2, we plot 3D and 2D graphs of the solution (2.38) which are shown in Fig.4 (c) and (d). Obviously, the probability density determined by the solution (2.36) accords with the normal distribution, but the probability density determined by the solution (2.38) does not accord with the normal distribution, as can be seen from Fig. 4(b) and (d), respectively.

Case of the quartic potential field
When V (x) = 1 4 mω 2 x 4 is a quartic potential and λ = − 3ω 2 ηα , equation (2.4) can be reduced to Under the transformation W = y(x)e − ω 2 x 4 4ηαKα , the equation (2.39) can be reduced to The general solution of equation (2.40) can be expressed by biconfluent Heun function as follows: , (2.41) where the function H(a, b, µ, ν, x) is biconfluent Heun function and the HeunB(a, b, µ, ν, x) is a special notation in the program of mathematical software Maple, the C 1 , C 2 are two arbitrary constants. Thus, the general solution of equation (2.39) is given by . (2.42) Plugging (2.42) and λ = − 3ω 2 ηα into (2.2), we obtain an analytical solution of the equation (1.5) as follows: . (2.43) When V (x) = 1 4 mω 2 x 4 + 1 2 mω 2 x 2 is a quartic double-well potential and λ = − 2ω 2 ηα , equation (2.4) can be reduced to As in the above computational process, we can obtain an analytical solution of the equation (1.5) as follows:

Case of the logarithmic potential field
When V (x) = mω 2 ln |x| is a logarithmic potential function and λ = − ω 2 ηα , equation (2.4) can be reduced to a generalized Bessel equation Under the variable transformation ξ = ω √ ηαKα x, the equation (2.46) can be reduced to a normal Bessel equation of 1-order as follows: its solution is where the J(1, ξ) is the first kind of Bessel function which is defined by ξ) is the second kind of Bessel function which is defined by with ψ(n + 1) = −γ + .
is a logarithmic potential , equation (2.4) can be reduced to a two-order ODE where C 1 , C 2 are two arbitrary constants. Plugging (2.51) into (2.2), we obtain an analytic solution of equation (1.5) as follows: is a logarithmic potential function and , equation (2.4) can be reduced to where C 1 , C 2 are two arbitrary constants. Thus, we obtain an analytic solution of equation (1.5) as follows:  The solution (2.52) implies compacton-type distribution. Obviously, the solution (2.52) is well suited to explain the phenomenon of the separation of large DNA molecules from the gel medium in a finite region, that is to say, the probability of the large DNA molecules separated from the gel medium in the central region of the separator is very high, but the probability near the edge of the separator becomes lower and lower.

Conclusions
In this work, based on the classical variable separation method, by using the properties of According to symmetrical characteristic, this method also can be used to investigate exact solutions and dynamical property of a space-fractional PDE. Of cause, this new approach is also suitable for solving high-dimensional time-fractional PDEs formed as  + ky), b(x, y) = b(x + ky) and c(x, y) = c(x + ky), that is to say, when the external force field is loaded on the plan x + ky, then we can make a new transformation of variable as follows: W (x, y) = W (ξ) with ξ = x + ky. obtained. Here we omit the detailed processes that are solved in the cases mentioned above.