Nonlinear Kinetics on Lattices Based on the Kinetic Interaction Principle

Abstract

Master equations define the dynamics that govern the time evolution of various physical processes on lattices. In the continuum limit, master equations lead to Fokker–Planck partial differential equations that represent the dynamics of physical systems in continuous spaces. Over the last few decades, nonlinear Fokker–Planck equations have become very popular in condensed matter physics and in statistical physics. Numerical solutions of these equations require the use of discretization schemes. However, the discrete evolution equation obtained by the discretization of a Fokker–Planck partial differential equation depends on the specific discretization scheme. In general, the discretized form is different from the master equation that has generated the respective Fokker–Planck equation in the continuum limit. Therefore, the knowledge of the master equation associated with a given Fokker–Planck equation is extremely important for the correct numerical integration of the latter, since it provides a unique, physically motivated discretization scheme. This paper shows that the Kinetic Interaction Principle (KIP) that governs the particle kinetics of many body systems, introduced in G. Kaniadakis, Physica A 296, 405 (2001), univocally defines a very simple master equation that in the continuum limit yields the nonlinear Fokker–Planck equation in its most general form.

1. Introduction

Nonlinear kinetics in the coordinate or velocity space or in the phase space has been used systematically over the last few decades. A popular approach involves studying the temporal evolution of a particle system while it interacts with a medium in thermodynamic equilibrium. In Ref. [1], the most general nonlinear kinetics in phase space was considered in the limit of the diffusive approximation. In this limit, the interaction between the particles and the medium is viewed as a diffusive process in the velocity space. It was shown that the particle probability density function $f=f(t,r,v)$ (henceforward particle density function) obeys the following evolution equation

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}+\frac{1}{m}\frac{\partial U}{\partial v}\frac{\partial f}{\partial r}-\frac{1}{m}\frac{\partial V}{\partial r}\frac{\partial f}{\partial v}=\frac{\partial }{\partial v}\left\{Da\left(f\right)b\left(f\right)\frac{\partial }{\partial v}\left[\beta V+\beta U+ln\frac{a\left(f\right)}{b\left(f\right)}\right]\right\},\end{array}$$

(1)

where $D=D(r,v)$ is the diffusion coefficient, while $V=V\left(r\right)$ and $U=U\left(v\right)$ represent the external potential and the particle kinetic energy respectively. In addition, $\partial f/\partial r$, $\partial f/\partial v$, denote the gradients of the particle density with respect to velocity and space, and m is the particle mass. The constant $\beta$ is given by $\beta =1/k_{\mathrm{B}}T$ where $k_{\mathrm{B}}$ is the Boltzmann constant and T the temperature of the system.

The above equation can be viewed as a generalized Kramers equation that describes a nonlinear kinetics in the phase space. The nonlinearity is due to the presence of the two arbitrary functions $a\left(f\right)$ and $b\left(f\right)$ that depend on the particle density function f. In the case where $a\left(f\right)=f$ and $b\left(f\right)=1$, Equation (1) reduces to the ordinary Kramers equation that describes the standard linear kinetics, e.g., [2].

Equation (1) has been employed to study the kinetic foundations of nonlinear statistical systems governed by $\kappa$-entropy. Generalized statistical mechanics, based on $\kappa$-entropy [1,3,4], preserves the main features of ordinary Boltzmann-Gibbs statistical mechanics. For this reason, it has attracted the interest of many researchers over the last 16 years, who have studied its foundations and mathematical aspects [5,6,7,8,9,10,11,12], the underlying thermodynamics [13,14,15,16,17], and specific applications of the theory in various scientific and engineering fields. A non-exhaustive list of application areas includes quantum statistics [18,19,20], quantum entanglement [21,22], plasma physics [23,24,25,26,27], nuclear fission [28], astrophysics [29,30,31,32,33,34,35], geomechanics [36], genomics [37], complex networks [38,39], economy [40,41,42,43] and finance [44,45,46,47].

It is important to note that the right hand side of Equation (1), which introduces nonlinear effects into kinetics, describes an unconventional diffusion process in the velocity space. This process can clearly also be studied in the coordinate space, where it describes an anomalous diffusion process. Such processes have long been observed and studied in various fields of condensed and soft matter physics [48].

Knowledge of the master equation that is associated with a given Fokker–Planck partial differential equation is crucial for two reasons. First, as stated above the master equation provides a uniquely defined discretization scheme for numerical integration. Second, the master equation is physically meaningful, since it defines the process dynamics on lattices. Hence, the study of the lattice dynamics allows an improved understanding of the dynamics in continuous space.

For example, consider the Fokker–Planck equation in the form (1) that involves two arbitrary functions $a\left(f\right)$ and $b\left(f\right)$. Numerical integration of this equation requires discretization in the spatial variable. The resulting master equation cannot contain any other functions besides $a\left(f\right)$ and $b\left(f\right)$, that are already present in the original partial differential equation. Hence, any lattice (finite-size) dynamical effects that are lost in the continuum limit cannot be re-introduced at this level. In addition, the master equation derived from the discretization of the Fokker–Planck equation depends on the particular differencing scheme (e.g., forward differences versus central differences). To avoid these shortcomings, it is worth investigating the use of physical principles at the lattice scale to univocally define the discrete kinetics and directly derive the master equation. The master equation thus obtained is clearly useful, not only for deriving the nonlinear Fokker–Planck Equation (1) in the continuum limit, but also as the starting point for the direct numerical integration of the Fokker–Planck equation.

The present paper shows that the Kinetic Interaction Principle (KIP) which governs the particle kinetics [1] univocally defines a very simple lattice-based master equation, which yields the nonlinear Fokker–Planck equation in the continuous limit. Both the master equation and the ensuing Fokker–Planck equation involve physically motivated terms.

The paper is organized as follows. In Section 2, we introduce the specific form of the Fokker–Planck equation that we will investigate. We then recall the KIP that underlies the kinetics of a particle system, and we apply it to derive the most general nonlinear master equation in the nearest-neighbour approximation. Furthermore we express the master equation in the form of a discrete continuity equation and introduce the nonlinear particle current. In Section 3, we calculate the particle current in the continuous limit and derive the respective Fokker–Planck Equation (2). In Section 4 we focus on standard finite-difference-based discretization schemes of the Fokker–Planck equation, and we show that in the nonlinear case the KIP-based master equation cannot be derived from such discretization schemes. In Section 5, we focus on an important subclass of nonlinear kinetics, described by Equation (2), which is related to ordinary Fickian diffusion. Finally, in Section 6 we present our concluding remarks and suggest future extensions of this research to lattices in higher dimensions.

2. Kinetic Interaction Principle and Master Equation

2.1. Nonlinear Fokker–Planck Kinetics in One Dimension

In the following, we focus on the one-dimensional nonlinear diffusive process that is described by the equation

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{\partial }{\partial x}\left\{Da\left(f\right)b\left(f\right)\frac{\partial }{\partial x}\left[\beta U+ln\frac{a\left(f\right)}{b\left(f\right)}\right]\right\}.\end{array}$$

(2)

where $f=f(t,x)$, $U=U\left(x\right)$ and $D=D\left(x\right)$, while x can indicate either the velocity or the coordinate space variable. Equation (2) represents a general form of the nonlinear Fokker–Planck equation. Upon setting $\gamma \left(f\right)=a\left(f\right)b\left(f\right)$, $\mathsf{\Omega }\left(f\right)=a\left(f\right)b\left(f\right)\frac{d}{df}ln\frac{a\left(f\right)}{b\left(f\right)}$, and $A=D\left(x\right)\beta \frac{\partial U\left(x\right)}{\partial x}$, this equation can also be expressed in the equivalent but more familiar form

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{\partial }{\partial x}\left[A\gamma \left(f\right)+D\mathsf{\Omega }\left(f\right)\frac{\partial f}{\partial x}\right],\end{array}$$

(3)

where the first summand inside the bracket on the right-hand side yields the drift term while the second summand yields the diffusive term.

It is in general difficult or even impossible to explicitly solve Fokker–Planck equations that are nonlinear with respect to the probability density function, or even linear ones for general forms of the external potential function $U\left(x\right)$. In the linear case it is possible to obtain explicit approximations for certain temporal characteristics of the Fokker–Planck equation, such as the escape time from a metastable potential function (by means of the Kramers approximation) or estimates of mean first passage times [49].

Several special cases of the nonlinear Fokker–Planck Equation (3) have long been known. The best known and most frequently studied case in the literature is that of diffusion in porous media, which corresponds to $A=0$, $D=$ constant and $\mathsf{\Omega }\left(f\right)=f^n$ [50]. Equation (3) (expressed in velocity space) has also been used in statistical physics, by setting $A=\beta v$, $D=$ constant and $\gamma \left(f\right)=f$. When $\mathsf{\Omega }\left(f\right)=1\pm f$, Equation (3) describes the kinetics of bosons (positive sign) or fermions (negative sign) [51], whereas when $\mathsf{\Omega }\left(f\right)=f^n$, it describes particle kinetics in non-extensive statistical mechanics [52].

In the case of bosonic and fermionic statistics, Equation (3) was obtained with standard methods used in linear kinetics, i.e., by applying the Kramers–Moyal expansion to a balance equation or alternatively starting from a master equation [51]. The same techniques were also used in [53,54] to derive Equation (3) starting from a master equation, in the special case where the function $\mathsf{\Omega }\left(f\right)$ is the sum of an arbitrary number of powers of the function f, i.e., $\mathsf{\Omega }\left(f\right)={\sum }_i{c}_i{f}^{{\alpha }_i}$.

The most general form of Equation (3) involves arbitrary functions $\gamma \left(f\right)$ and $\mathsf{\Omega }\left(f\right)$. The study of this general form is important for classifying the stationary and thermodynamically stable solutions [1,55,56,57]. On the other hand, the derivation of the nonlinear Fokker–Planck equation from physical principles is crucial for understanding the nature of the $\gamma \left(f\right)$ and $\mathsf{\Omega }\left(f\right)$ functions that determine the solutions.

The nonlinear Fokker–Planck equation in the most general possible from (2) has been obtained by applying the Kramers–Moyal expansion (that takes into account only transitions between near neighbors) to a balance equation based on the (KIP) [1].

Subsequently, Equation (3) was obtained in [58] from a master equation that involves nearest-neighbor transitions, the rates of which are determined by two arbitrary functions $\mathrm{a}\left(f\right)$ and $Y(f,f^{\prime })$. The first function, $\mathrm{a}\left(f\right)$, depends on the particle density function f and is related to the $\gamma \left(f\right)$ function through the simple relation $\gamma \left(f\right)=f\mathrm{a}\left(f\right)$. On the other hand, $Y(f,f^{\prime })$ depends on two distinct particle density functions, f and $f^{\prime }$ and is related to the $\mathsf{\Omega }\left(f\right)$ function by means of a more complicated expression.

2.2. The Master Equation

Let us consider the particle kinetics in a one-dimensional lattice gas comprising N identical particles. Under the hypothesis that only transitions involving the nearest-neighbour lattice sites are allowed during the evolution of the particle system, the probability density function $f_i=f(t,x_i)$ obeys the following master equation

$$\begin{array}{c}\hfill \frac{d{f}_i}{dt}=\pi (t,x_{i-1}\to x_i)-\pi (t,x_i\to x_{i-1})+\pi (t,x_{i+1}\to x_i)-\pi (t,x_i\to x_{i+1}),\end{array}$$

(4)

where $\pi (t,x_i\to x_j)$ is the transition probability from site i to site j.

The transition probability related to the hopping of a particle from site i to site j depends on the nature of the interaction between the lattice and the particle, and this dependence is taken into account through the transition rates $w_{ij}$. The transition probability can of course depend on the particle density of the starting i and arrival j sites.

In Ref. [1], it was postulated that the $f_i$ and $f_j$ populations, of the starting and arrival sites affect the transition probability through the two arbitrary functions $a\left(f_i\right)$, and $b\left(f_j\right)$, so that $\pi (t,x_i\to x_j)$ assumes the following factorized form

$$\pi (t,x_i\to x_j)=w_{ij}a\left(f_i\right)b\left(f_j\right),$$

(5)

and this particular dependence on the particle population of the starting and arrival sites defines the Kinetic Interaction Principle.

2.3. The Kinetic Interaction Principle

Let us now consider the transition probability from site i to site j. The $a\left(f_i\right)$ function must satisfy the obvious condition $a\left(0\right)=0$, because in the case where the starting site is empty no particles can transit toward the arrival site. On the other hand the $b\left(f_j\right)$ function must satisfy the condition $b\left(0\right)=1$, because in the case where the arrival site is empty the transition probability cannot be influenced by the arrival site.

In linear kinetics, $a\left(f_i\right)=f_i$ and $b\left(f_j\right)=1$, is usually posed [59]. In nonlinear kinetics the $b\left(f_j\right)$ function plays an important role, because it can stimulate or inhibit the particle transition to site j and can simulate collective interactions in many body physics. For instance, the expression $b\left(f_j\right)=1-f_j$ can account for the Pauli exclusion principle that governs a fermionic particle system [51,60,61]. Analogously, the expression $b\left(f_j\right)=1+f_j$ accounts for bosonic systems. Other more complicated expressions of the $b\left(f_j\right)$ factor can take into account collective interactions introduced by the Haldane generalized exclusion principle originating from the fractal structure of the single particle Hilbert space.

Another case of nonlinear kinetics corresponding to $a\left(f_i\right)=f_i{\mathrm{e}}^{-u{f}_i}$ where $u>0$ and $b\left(f_j\right)=1$ was investigated in [62], motivated by observations of two different grain growth regimes in sintering studies [63]. In this type of master equation the value of u controls the equilibrium state: a normal diffusive regime is obtained for u below a threshold, while an abnormal diffusion regime that exhibits the Matthew effect of accumulated advantage is obtained at higher u.

If transitions are allowed only among the nearest lattice neighbors, i.e., $i\to j=i\pm 1$, the transition rates $w_{ij}$ are non-zero only for nearest-neighbor transitions. The nearest-neighbor transition rates are $w_i^{\pm }$, where $w_i^{+}$ corresponds to the transition $i\to i+1$ and $w_i^{-}$ corresponds to the transition $i\to i-1$, the KIP assumes the following form

$$\pi (t,x_i\to x_{i\pm 1})=w_i^{\pm }a\left(f_i\right)b\left(f_{i\pm 1}\right),$$

(6)

while the master Equation (4) becomes

$$\begin{array}{c}\hfill \frac{d{f}_i}{dt}=w_{i-1}^{+}a\left(f_{i-1}\right)b\left(f_i\right)-w_i^{-}a\left(f_i\right)b\left(f_{i-1}\right)+w_{i+1}^{-}a\left(f_{i+1}\right)b\left(f_i\right)-w_i^{+}a\left(f_i\right)b\left(f_{i+1}\right).\end{array}$$

(7)

2.4. Incoming and Outgoing Lattice Currents

Let us express the master Equation (7) in the form

$$\begin{array}{c}\hfill \frac{\partial f_i}{\partial t}+\frac{j_i^{+}-j_i^{-}}{\Delta x}=0,\end{array}$$

(8)

where the outgoing, $j_i^{+}$, and incoming, $j_i^{-}$, currents are defined as follows

$$\begin{array}{c}{j}_i^{+}=\Delta x\left[w_i^{+}a\left(f_i\right)b\left(f_{i+1}\right)-w_{i+1}^{-}a\left(f_{i+1}\right)b\left(f_i\right)\right],\hfill \end{array}$$

(9)

$$\begin{array}{c}{j}_i^{-}=\Delta x\left[w_{i-1}^{+}a\left(f_{i-1}\right)b\left(f_i\right)-w_i^{-}a\left(f_i\right)b\left(f_{i-1}\right)\right].\hfill \end{array}$$

(10)

Taking into account the continuity relation between the incoming current at the node i and the outgoing current at node $i-1$, that is,

$$\begin{array}{c}\hfill j_i^{-}=j_{i-1}^{+},\end{array}$$

(11)

the difference $j_i^{+}-j_i^{-}$ between the outgoing and the incoming currents can be expressed exclusively in terms of the outgoing current (or alternatively in terms of the incoming current), according to

$$j_i^{+}-j_i^{-}={\Delta }^{\pm }{j}_i^{\mp },$$

(12)

where the forward difference ${\Delta }^{+}{g}_i$ and the backward difference ${\Delta }^{-}{g}_i$ of the discrete function $g_i$ are defined through

$${\Delta }^{+}{g}_i=g_{i+1}-g_i,$$

(13a)

$${\Delta }^{-}{g}_i=g_i-g_{i-1}.$$

(13b)

2.5. Continuity form of the Master Equation

Based on the definition of the incoming and outgoing currents given above, the master Equation (8) assumes the following compact form

$$\frac{\partial f_i}{\partial t}+\frac{{\Delta }^{\pm }{j}_i^{\mp }}{\Delta x}=0.$$

(14)

The two forms of the master equation, as given by Equation (14), are asymmetric with respect to the currents. By employing the symmetry ${\Delta }^{-}{j}_i^{+}={\Delta }^{+}{j}_i^{-}$, the master equation can be expressed in a symmetric form that involves both the outgoing $j_i^{+}$ and incoming $j_i^{-}$ currents, as well as both the forward ${\Delta }^{+}$ and backward ${\Delta }^{-}$ differences, i.e.,

$$\frac{\partial f_i}{\partial t}+\frac{1}{2}\left(\frac{{\Delta }^{+}{j}_i^{-}}{\Delta x}+\frac{{\Delta }^{-}{j}_i^{+}}{\Delta x}\right)=0.$$

(15)

To write the latter equation in an even more symmetric form we introduce the discrete Fokker–Planck current, $j_i$, according to the average

$$j_i=\frac{j_i^{+}+j_i^{-}}{2}.$$

(16)

When this definition is employed, the incoming $j_i^{-}$ and the outgoing $j_i^{+}$ currents can be expressed in terms of the $j_i$ current

$$j_i^{\pm }=j_i\pm \frac{\Delta x}{2}\frac{j_i^{+}-j_i^{-}}{\Delta x},$$

(17)

and after taking into account the master Equation (8), we obtain

$$j_i^{\pm }=j_i\mp \frac{\Delta x}{2}\frac{\partial f_i}{\partial t}.$$

(18)

Finally, the following first- and second-order symmetric finite difference operators are introduced

$$\begin{array}{c}\hfill \frac{\Delta }{\Delta x}=\frac{1}{2}\left(\frac{{\Delta }^{+}}{\Delta x}+\frac{{\Delta }^{-}}{\Delta x}\right),\end{array}$$

(19)

$$\begin{array}{c}\hfill \frac{{\Delta }^2}{{(\Delta x)}^2}=\frac{1}{\Delta x}\left(\frac{{\Delta }^{+}}{\Delta x}-\frac{{\Delta }^{-}}{\Delta x}\right),\end{array}$$

(20)

and after taking into account Equation (18), the master Equation (15) assumes the form

$$\begin{array}{c}\hfill \frac{\partial f_i}{\partial t}+\frac{\Delta j_i}{\Delta x}+{(\Delta x)}^2\frac{1}{4}\frac{\partial }{\partial t}\frac{{\Delta }^2{f}_i}{{(\Delta x)}^2}=0.\end{array}$$

(21)

The latter expression of the master equation is particulary suitable for obtaining the continuous limit as $\Delta x\to 0$. First, we introduce the particle distribution function $f(t,x)$, according to

$$\begin{array}{c}\hfill f(t,x)=\underset{\Delta x\to 0}{lim}{f}_i,\end{array}$$

(22)

and the Fokker–Planck current $j(t,x)$, through

$$\begin{array}{c}\hfill j(t,x)=\underset{\Delta x\to 0}{lim}{j}_i.\end{array}$$

(23)

Subsequently, after taking into accounts the limits ${\Delta }^{+}/\Delta x\to \partial /\partial x$, ${\Delta }^{-}/\Delta x\to \partial /\partial x$, $\Delta /\Delta x\to \partial /\partial x$, ${\Delta }^2/{(\Delta x)}^2\to {\partial }^2/{(\partial x)}^2$, both Equations (15) and (21) reduce to the continuity equation

$$\begin{array}{c}\hfill \frac{\partial f(t,x)}{\partial t}+\frac{\partial j(t,x)}{\partial x}=0.\end{array}$$

(24)

The above continuity equation is the nonlinear Fokker–Planck equation obtained in the continuous limit, starting from the master Equation (7). To write the nonlinear Fokker–Planck equation explicitly, as last step we have to calculate the Fokker–Planck current $j(t,x)$, starting from its definition (16), and this will be the task of the next section.

3. The Nonlinear Fokker–Planck Current

3.1. Lattice Expression of the Fokker–Planck Current

To calculate the continuous limit of the discrete Fokker–Planck current defined in Equation (23), we express $j_i$, defined through Equations (9), (10) and (16), in terms of the functions $a(·)$ and $b(·)$ thus obtaining

$$\begin{array}{cc}\hfill j_i=& \frac{\Delta x}{2}[w_i^{+}a\left(f_i\right)b\left(f_{i+1}\right)+w_{i-1}^{+}a\left(f_{i-1}\right)b\left(f_i\right)\hfill \\ & \hfill -w_i^{-}a\left(f_i\right)b\left(f_{i-1}\right)-w_{i+1}^{-}a\left(f_{i+1}\right)b\left(f_i\right)]& .\hfill \end{array}$$

(25)

To express the current $j_i$ in a form that can be more efficiently used to compute the continuous limit, we first note that the transition rates $w_{i\mp 1}^{\pm }$ can be expanded in terms of the nearest-neighbour rates $w_i^{\pm }$ as follows

$$\begin{array}{c}\hfill w_{i\mp 1}^{\pm }=w_i^{\pm }\mp \frac{{\Delta }^{\mp }{w}_i^{\pm }}{\Delta x}\Delta x.\end{array}$$

(26)

Furthermore, we assume that the functions $a(·)$ and $b(·)$ are differentiable functions of $f_i$. Then, the $a\left(f_{i+1}\right)$ function can be approximated using the nearest-neighbour approximation of $f_{i+1}$ in terms of $f_i$ and the leading-order Taylor expansion of $a\left(f_{i+1}\right)$ around $f_i$ which lead to

$$\begin{array}{cc}\hfill a\left(f_{i+1}\right) & =a\left(f_i+\frac{f_{i+1}-f_i}{\Delta x}\Delta x\right)\hfill \\ & =a\left(f_i+\frac{{\Delta }^{+}{f}_i}{\Delta x}\Delta x\right)\hfill \\ & \approx a\left(f_i\right)+\frac{da\left(f_i\right)}{d{f}_i}\frac{{\Delta }^{+}{f}_i}{\Delta x}\Delta x.\hfill \end{array}$$

(27)

Following the same procedure, one obtains the approximate expressions for the functions $a(·)$, $b(·)$, at the nearest-neighbour sites of a given site i

$$\begin{array}{c}a\left(f_{i\pm 1}\right)\approx a\left(f_i\right)\pm \frac{da\left(f_i\right)}{d{f}_i}\frac{{\Delta }^{\pm }{f}_i}{\Delta x}\Delta x,\hfill \end{array}$$

(28)

$$\begin{array}{c}b\left(f_{i\pm 1}\right)\approx b\left(f_i\right)\pm \frac{db\left(f_i\right)}{d{f}_i}\frac{{\Delta }^{\pm }{f}_i}{\Delta x}\Delta x.\hfill \end{array}$$

(29)

After substitution of the above expressions for $a\left(f_{i\pm 1}\right)$ and $b\left(f_{i\pm 1}\right)$ in the equation Equation (25) of the discrete Fokker–Planck current $j_i$, one obtains

$$\begin{array}{cc}\hfill j_i\hspace{1em}& \approx \left[(w_i^{+}-w_i^{-})\Delta x-\frac{{\Delta }^{-}{w}_i^{+}+{\Delta }^{+}{w}_i^{-}}{\Delta x}\frac{{(\Delta x)}^2}{2}\right]a\left(f_i\right)b\left(f_i\right)\\ & +a\left(f_i\right)\frac{db\left(f_i\right)}{d{f}_i}\left[w_i^{+}\frac{{\Delta }^{+}{f}_i}{\Delta x}+w_i^{-}\frac{{\Delta }^{-}{f}_i}{\Delta x}\right]\frac{{(\Delta x)}^2}{2}\hfill \\ & -b\left(f_i\right)\frac{da\left(f_i\right)}{d{f}_i}\left[w_i^{+}\frac{{\Delta }^{-}{f}_i}{\Delta x}+w_i^{-}\frac{{\Delta }^{+}{f}_i}{\Delta x}\right]\frac{{(\Delta x)}^2}{2}.\hfill \end{array}$$

(30)

3.2. Continuum Limit Expression of the Fokker–Planck Current

The above expression for the discrete current can be simplified in the continuum limit, $\Delta x\to 0$, by taking into account that there is no difference between the first-order forward, backward, and symmetric finite-differences, i.e.,

$$\begin{array}{c}\hfill \frac{{\Delta }^{\pm }{g}_i}{\Delta x}\approx \frac{\Delta g_i}{\Delta x}.\end{array}$$

(31)

Hence, the discrete Fokker–Planck current $j_i$ assumes the form

$$\begin{array}{cc}\hfill j_i\hspace{1em}& \approx \left[(w_i^{+}-w_i^{-})\Delta x-\frac{\Delta \left(w_i^{+}+w_i^{-}\right)}{\Delta x}\frac{{(\Delta x)}^2}{2}\right]a\left(f_i\right)b\left(f_i\right)\hfill \\ & +(w_i^{+}+w_i^{-})\frac{{(\Delta x)}^2}{2}\left[a\left(f_i\right)\frac{db\left(f_i\right)}{d{f}_i}-b\left(f_i\right)\frac{da\left(f_i\right)}{d{f}_i}\right]\frac{\Delta f_i}{\Delta x}.\hfill \end{array}$$

(32)

The diffusion coefficient, defined by means of

$$\begin{array}{c}\hfill D_i=\frac{w_i^{+}+w_i^{-}}{2}{(\Delta x)}^2,\end{array}$$

(33)

and the drift coefficient, defined by means of

$$\begin{array}{c}\hfill J_i=(w_i^{-}-w_i^{+})\Delta x,\end{array}$$

(34)

are usually introduced into Fokker-Plank kinetics. Following these replacements, the discrete Fokker–Planck current given by Equation (32) assumes the equivalent form

$$\begin{array}{ccc}\hfill j_i=& -& \left(J_i+\frac{\Delta D_i}{\Delta x}\right)a\left(f_i\right)b\left(f_i\right)-D_i\left[b\left(f_i\right)\frac{da\left(f_i\right)}{d{f}_i}-a\left(f_i\right)\frac{db\left(f_i\right)}{d{f}_i}\right]\frac{\Delta f_i}{\Delta x}.\hfill \end{array}$$

(35)

Furthermore, the introduction of the kinetic energy function $U_i$ according to

$$\begin{array}{c}\hfill \beta \frac{\Delta U_i}{\Delta x}=\frac{1}{D_i}\left(J_i+\frac{\Delta D_i}{\Delta x}\right),\end{array}$$

(36)

allows expressing the current in the following form

$$j_i=-D_{i}a\left(f_i\right)b\left(f_i\right)\left[\beta \frac{\Delta U_i}{\Delta x}+\frac{d}{d{f}_i}ln\frac{a\left(f_i\right)}{b\left(f_i\right)}\frac{\Delta f_i}{\Delta x}\right],$$

(37)

or equivalently as

$$j_i=-D_{i}a\left(f_i\right)b\left(f_i\right)\frac{\Delta }{\Delta x}\left[\beta U_i+ln\frac{a\left(f_i\right)}{b\left(f_i\right)}\right].$$

(38)

The discrete current $j_i$, as given by the latter equation, is expressed in a form that permits straightforward evaluation of the $\Delta x\to 0$ continuous limit. Thus, the nonlinear Fokker–Planck current in the continuum limit is given by the following equation

$$j(t,x)=-D(t,x)a\left[f(t,x)\right]b\left[f(t,x)\right]\frac{\partial }{\partial x}\left\{\beta U(t,x)+ln\frac{a\left[f\right(t,x\left)\right]}{b\left[f\right(t,x\left)\right]}\right\}.$$

(39)

4. Fokker–Planck Equation and Discretization Schemes

If we insert the continuum limit of the Fokker–Planck $j(t,x)$ current, given by Equation (39), in the continuity Equation (24) we obtain the nonlinear Fokker–Planck equation

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{\partial }{\partial x}\left\{Da\left(f\right)b\left(f\right)\frac{\partial }{\partial x}\left[\beta U+ln\frac{a\left(f\right)}{b\left(f\right)}\right]\right\}.\end{array}$$

(40)

The above equation recovers the original Fokker–Planck equation given by Equations (2) and (3). Thus, the derivations in the preceding section show how the nonlinear Fokker–Planck equation can be obtained from the discrete Fokker–Planck current of Equation (23), which is derived from the simple, KIP-based master Equation (7).

In the following, we investigate the impact of the discretization of the time derivative in Fokker–Planck equation, as well as the impact of the space derivative discretization in the linear and nonlinear kinetic regimes. We also present an extension of the one-dimensional formalism to two spatial dimensions.

We will express the Fokker–Planck Equation (40) as

$$\begin{array}{cc}\hfill \frac{\partial f(t,x)}{\partial t}=& F(t,x),\hfill \end{array}$$

(41a)

$$\hspace{1em}F(t,x)=\frac{\partial G(t,x)}{\partial x},$$

(41b)

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}G(t,x)=Da\left(f\right)b\left(f\right)\frac{\partial }{\partial x}\left[\beta U+ln\frac{a\left(f\right)}{b\left(f\right)}\right].$$

(41c)

There are various numerical schemes for the discretization of partial differential equations such as the above, including finite-difference schemes, methods based on finite elements, finite-volume methods, spectral approaches, multigrid methods, etc. [64]. Each approach leads to a different set of equations that depends on the mathematical assumptions made. Different discretization schemes for the one-dimensional Fokker–Planck equation are presented in [65].

It is most straightforward to compare the master Equation (7) with discretized Fokker–Planck equations that are obtained by means of finite-difference methods. Discretized versions of the functions $f(t,x)$ are denoted by means of $f_i^n$, where the index i refers to the lattice site and the index n to the time instant.

4.1. Temporal Discretization

Three of the most common finite-difference methods lead to the following discretization schemes of (41a)

$$\begin{array}{cc}\frac{f_i^{n+1}-f_i^n}{\Delta t}=F_i^n,\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{Forward} \mathrm{Euler}\end{array}$$

(42a)

$$\begin{array}{cc}\frac{f_i^{n+1}-f_i^n}{\Delta t}=F_i^{n+1},\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{Backward} \mathrm{Euler}\end{array}$$

(42b)

$$\begin{array}{cc}\frac{f_i^{n+1}-f_i^n}{\Delta t}=\frac{1}{2}\left(F_i^n+F_i^{n+1}\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\text{Crank–Nicolson}\end{array}$$

(42c)

where $\Delta t$ is the time step.

On the other hand, the discretization of the master Equation (7) leads to

$$\frac{f_i^{n+1}-f_i^n}{\Delta t}={\tilde{F}}_i^n,$$

(43)

$${\tilde{F}}_i^n=w_{i-1}^{+n}a\left(f_{i-1}^n\right)b\left(f_i^n\right)-w_i^{-n}a\left(f_i^n\right)b\left(f_{i-1}^n\right)+w_{i+1}^{-n}a\left(f_{i+1}^n\right)b\left(f_i^n\right)-w_i^{+n}a\left(f_i^n\right)b\left(f_{i+1}^n\right)$$

(44)

If we compare the master Equations (43) with (42) we notice that the right-hand side of the former involves the function $F_i^n$ only at the current time step. In contrast, the backward Euler and the Crank-Nicolson schemes involve the values of the function, $F_i^{n+1}$, at the next time step as well.

In the forward Euler discretization the time difference is determined by the function $F_i^n$, which is in general different from ${\tilde{F}}_i^n.$ The function $F_i^n$ is obtained by discretizing the spatial partial derivatives in $F(t,x)$. Using the forward spatial derivative defined in (13), the function ${\tilde{F}}_i^n$ can be shown to depend on the values of $a(·),b(·),U(·)$ at the sites labeled by $i,i+1,i+2$. In contrast, the time difference of the discretized master Equation (43) depends on the lattices sites i and $i\pm 1$.

4.2. Linear Kinetics Regime

The linear kinetics regime of (2) is obtained if $a\left(f\right)=f$ and $b\left(f\right)=1$. After recalling that $\gamma \left(f\right)=a\left(f\right)b\left(f\right)$, $\mathsf{\Omega }\left(f\right)=a\left(f\right)b\left(f\right)\frac{\partial }{\partial f}ln\frac{a\left(f\right)}{b\left(f\right)}$, it follows that $\gamma \left(f\right)=f$ and $\mathsf{\Omega }\left(f\right)=1$. Hence, the linear Fokker–Planck equation derived from (2) assumes the form:

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{dA\left(x\right)}{dx}f+\left[A\left(x\right)+\frac{dD\left(x\right)}{dx}\right]\frac{\partial f}{\partial x}+D\left(x\right)\frac{{\partial }^{2}f}{\partial x^2},\end{array}$$

(45)

with $A=D\left(x\right)\beta \frac{\partial U\left(x\right)}{\partial x}$.

To discretize the above equation in space, we introduce at the place of the first and second order partial derivatives the following symmetric finite differences:

$$\frac{\partial f}{\partial x}\approx \frac{f_{i+1}-f_{i-1}}{2\Delta x},$$

(46a)

$$\frac{{\partial }^{2}f}{\partial x^2}\approx \frac{f_{i+1}-2f_i+f_{i-1}}{{(\Delta x)}^2},$$

(46b)

where $\Delta x$ is the uniform discretization step in space.

Then, the linear Fokker–Planck Equation (45) transforms into the linear master equation

$$\begin{array}{c}\hfill \frac{d{f}_i}{dt}=w_{i-1}^{+}{f}_{i-1}+w_{i+1}^{-}{f}_{i+1}-w_i{f}_i,\end{array}$$

(47a)

with the nearest-neighbor transition rates $w_{i+1}^{-}$ and $w_{i-1}^{+}$ defined by

$$w_{i\pm 1}^{\mp }=\frac{D_i}{{(\Delta x)}^2}\pm \frac{D_{i+1}-D_{i-1}}{4{(\Delta x)}^2}\pm \frac{A_i}{2\Delta x},$$

(47b)

$$w_i=\frac{2D_i}{{(\Delta x)}^2}-\frac{A_{i+1}-A_{i-1}}{2\Delta x}.$$

(47c)

Similarly, the nearest-neighbor master Equation (7) is given in the linear regime by the following equation

$$\frac{d{f}_i}{dt}=w_{i-1}^{+}{f}_{i-1}+w_{i+1}^{-}{f}_{i+1}-\left(w_i^{+}+w_i^{-}\right)f_i.$$

(48)

After observing that $w_i\approx w_i^{+}+w_i^{-}$, the linear master Equation (48) derived by means of the KIP is essentially equivalent to the master Equation (47a) obtained from the discretization of the linear Fokker–Planck Equation (45).

We can thus conclude that in the case of linear kinetics the KIP-based master equation leads in the continuum limit to a Fokker–Planck equation. Discretization of the Fokker–Planck equation using the centered-difference scheme yields the original master equation.

In the case of linear kinetics we can use either the forward or the backward Euler method, or the Cranck-Nicolson method, since the coefficients of the Fokker–Planck equation do not depend on time. In addition, the linear dependence on the probability density function ensures that even for the more demanding Crank-Nicolson scheme, it suffices to solve a linear system.

4.3. Nonlinear Kinetics

In the case of nonlinear kinetics, we consider the Fokker–Planck equation as given in (3). Expansion of the terms that involve the spatial derivatives lead to the following expression

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{dA\left(x\right)}{dx}\gamma \left(f\right)+\left[A\left(x\right)\frac{d\gamma \left(f\right)}{df}+\frac{dD\left(x\right)}{dx}\mathsf{\Omega }\left(f\right)\right]\frac{\partial f}{\partial x}+D\left(x\right)\frac{d\mathsf{\Omega }\left(f\right)}{df}{\left(\frac{\partial f}{\partial x}\right)}^2+D\left(x\right)\mathsf{\Omega }\left(f\right)\frac{{\partial }^{2}f}{\partial x^2}.\end{array}$$

(49)

Using the centered finite difference approximation of the spatial derivatives for the functions $f(t,x)$, $D\left(x\right)$ and $A\left(x\right)$ we obtain the following master equation

$$\begin{array}{c}\hfill \frac{d{f}_i}{dt}={\omega }^{+}\left(f_i\right)f_{i-1}+{\theta }_i\left(f_i\right){\left(f_{i-1}\right)}^2+{\omega }^{-}\left(f_i\right)f_{i+1}+{\theta }_i\left(f_i\right){\left(f_{i+1}\right)}^2-\omega \left(f_i\right)f_i-2{\theta }_i\left(f_i\right)f_{i+1}{f}_{i-1},\end{array}$$

(50a)

where the density-dependent coefficients ${\omega }^{\mp }\left(f_i\right)$, $\omega \left(f_i\right)$, and ${\theta }_i\left(f_i\right)$ are given by

$${\omega }^{\mp }\left(f_i\right)=\frac{D_i}{{(\Delta x)}^2}\mathsf{\Omega }\left(f_i\right)\pm \frac{D_{i+1}-D_{i-1}}{4{(\Delta x)}^2}\mathsf{\Omega }\left(f_i\right)\pm \frac{A_i}{2\Delta x}\frac{d\gamma \left(f_i\right)}{d{f}_i},$$

(50b)

$$\omega \left(f_i\right)=\frac{2D_i}{{(\Delta x)}^2}\mathsf{\Omega }\left(f_i\right)-\frac{A_{i+1}-A_{i-1}}{2\Delta x}\frac{\gamma \left(f_i\right)}{f_i},$$

(50c)

$${\theta }_i\left(f_i\right)=\frac{D_i}{4{(\Delta x)}^2}\frac{d\mathsf{\Omega }\left(f_i\right)}{d{f}_i}.$$

(50d)

The master Equation (7) that is based on the KIP with nearest-neighbor transitions represents different kinetics than the master Equation (50) that is obtained from the discretization of the nonlinear Fokker–Planck Equation (3). The kinetics described by the KIP-based master equation is physically motivated at the microscopic level. On the other hand, it is not straightforward to interpret in physical terms the kinetics described by (50). More precisely, the first two terms on the right-hand side of (50a) can be viewed as representing transitions into the site i from the site $i-1$. The population of the departure site enters linearly the transition rate in the first term, but the population of the arrival site i can be a highly nonlinear and complicated function. In the second term, the dependence on the departure site is also nonlinear in addition to that of the arrival site. Analogous remarks hold for the third and fourth terms which represent transitions from the site $i+1$ towards the site i. The last two terms represent exiting transitions from the sites i. However, there is no obvious physical mechanism that explains why the transition rate for particles exiting the site i should be proportional to the product of the densities at both sites $i-1$ and $i+1$ (as dictated by the term ${\theta }_i\left(f_i\right)f_{i+1}{f}_{i-1}$). Clearly, the master Equation (50) does not have an obvious microscopic interpretation, even though in the continuum limit it yields the nonlinear Fokker–Planck Equation (3).

Based on the above observations, in order to integrate numerically the nonlinear Fokker–Planck equation it is recommended to start with the spatial discretization provided by the master Equation (7). For the time step, the forward or backward Euler methods are simpler alternatives than the Crank–Nicolson method which leads to a nonlinear system of algebraic equations with respective computational cost. More importantly, in order to study the statistical properties of a physical system governed by lattice nonlinear kinetics, the master Equation (7) is clearly the starting point, independently of the method used to integrate it (which is not the focus of the current study).

How should we understand the difference between the master Equation (7) and that produced by the discretization of the nonlinear Fokker–Planck Equation (3), which after all is obtained from the initial master Equation (7) in the continuum limit? The key is that in process of taking the continuum limit information is lost. The master Equation (7) contains the density functions $f_i$, $f_{i-1}$, $f_{i+1}$ of three sites $i,i\pm 1$ through the set of functions $I_M=\{a\left(f_i\right),b\left(f_i\right),a\left(f_{i-1}\right),b\left(f_{i-1}\right),a\left(f_{i+1}\right),b\left(f_{i+1}\right)\}$. The nonlinear Fokker–Planck Equation (3) on the other hand, contains the function $f\left(x\right)$ through its dependence on $a\left(f\right)$, $b\left(f\right)$ and the partial derivatives $\partial f/\partial x$, ${\partial }^{2}f/\partial x^2$. Clearly the set $I_M$ that is used in the master Equation (7) contains more information than the set $I_{FP}$ = {$a\left(f\right)$, $b\left(f\right)$, $\partial f/\partial x$ and ${\partial }^{2}f/\partial x^2$}.

Due to this loss of information in the transition from the discrete (lattice-based) to the continuum model, discretization of the latter generates the master Equation (50), which is radically different from the master Equation (7), although both master equations yield the same Fokker–Planck equation in the continuum limit. It is important to stress that more than one master equations, including (7) imposed by the KIP, yield the same Fokker–Planck equation in the continuum limit. In contrast, applying the standard discretization rules to the Fokker–Planck equation, the master equations derived are not physically motivated and they differ from the master equation generated by the KIP. Only in the case of the linear Fokker–Planck equation, the master equation obtained by discretization coincides with the KIP-based master equation.

5. Drift Current and Diffusion in Nonlinear Kinetics

5.1. Nonlinear Drift and Fick Currents

Let us now write the Fokker–Plank Equation (40) in the following form

$$\frac{\partial f}{\partial t}+\frac{\partial (j_{\mathrm{drift}}+j_{\mathrm{Fick}})}{\partial x}=0,$$

(51)

where the nonlinear drift current $j_{\mathrm{drift}}=j_{\mathrm{drift}}(t,x)$ and the nonlinear Fick current $j_{\mathrm{Fick}}=j_{\mathrm{Fick}}(t,x)$ are defined according to

$$\begin{array}{c}\hfill j_{\mathrm{drift}}=-\beta D\frac{\partial U}{\partial x}a\left(f\right)b\left(f\right),\end{array}$$

(52)

$$\begin{array}{c}\hfill j_{\mathrm{Fick}}=-D\varphi \left(f\right)\frac{\partial f}{\partial x},\end{array}$$

(53)

with

$$\begin{array}{c}\hfill \varphi \left(f\right)=a\left(f\right)b\left(f\right)\frac{d}{df}ln\frac{a\left(f\right)}{b\left(f\right)}.\end{array}$$

(54)

The term $D\varphi \left(f\right)$ represents the nonlinear diffusion coefficient. If $\varphi \left(f\right)=1$, then the ordinary Fick current $j_{\mathrm{Fick}}=-D\frac{\partial f}{\partial x}$ is obtained. The ordinary Fick current appears in (i) linear kinetics for $a\left(f\right)=f$ and $b\left(f\right)=1$, and (ii) in the classical models of boson or fermion kinetics, for $a\left(f\right)=f$ and $b\left(f\right)=1\pm f$ respectively.

This important result, pertaining to the ordinary diffusion that underlies the kinetics of classical bosons and fermions, naturally leads to the question whether there exist other nonlinear kinetics which are governed by standard Fickian diffusion.

5.2. Nonlinear Kinetics with Fickian Diffusion

In the following, we focus on the Fokker–Planck Equation (40) in the case of Fickian diffusion i.e., if the $\varphi \left(f\right)$ function defined in Equation (54) is subject to the $\varphi \left(f\right)=1$ condition. In nonlinear kinetics, the generalized logarithm $\mathsf{\Lambda }\left(f\right)$ is introduced through

$$\begin{array}{c}\hfill \mathsf{\Lambda }\left(f\right)=ln\frac{a\left(f\right)}{b\left(f\right)}.\end{array}$$

(55)

In terms of the generalized logarithm, the stationary and stable solution, $f_s$, of Equation (40) assumes the form

$$\begin{array}{c}\hfill \mathsf{\Lambda }\left(f_s\right)=-\beta [U\left(x\right)-\mu ],\end{array}$$

(56)

where $\mu$ is an arbitrary constant.

It is easy to express the $a\left(f\right)$ and $b\left(f\right)$ functions in terms of the $\mathsf{\Lambda }\left(f\right)$ function. First, we write Equation (55) in the form

$$\begin{array}{c}\hfill \frac{a\left(f\right)}{b\left(f\right)}=exp\left[\mathsf{\Lambda }\left(f\right)\right].\end{array}$$

(57)

Then, based on Equation (54) and taking into account the $\varphi \left(f\right)=1$ condition, we obtain

$$\begin{array}{c}\hfill a\left(f\right)b\left(f\right)=\frac{1}{d\mathsf{\Lambda }\left(f\right)/df}.\end{array}$$

(58)

Finally, the solutions of the two equations above for $a\left(f\right)$ and $b\left(f\right)$ lead to

$$\begin{array}{c}a\left(f\right)={\left[\frac{d\mathsf{\Lambda }\left(f\right)}{df}\right]}^{-1/2}exp\left[\frac{1}{2}\mathsf{\Lambda }\left(f\right)\right],\hfill \end{array}$$

(59)

$$\begin{array}{c}b\left(f\right)={\left[\frac{d\mathsf{\Lambda }\left(f\right)}{df}\right]}^{-1/2}exp\left[-\frac{1}{2}\mathsf{\Lambda }\left(f\right)\right].\hfill \end{array}$$

(60)

The nonlinear Fokker–Planck kinetics described by Equation (40), when $a\left(f\right)$ and $b\left(f\right)$ are given by Equations (59) and (60), simplifies to

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{\partial }{\partial x}\left[D\beta \frac{\partial U}{\partial x}\gamma \left(f\right)+D\frac{\partial f}{\partial x}\right],\end{array}$$

(61)

where the function $\gamma \left(f\right)$ represents the inverse of the derivative of the generalized logarithm, given by

$$\begin{array}{c}\hfill \frac{1}{\gamma \left(f\right)}=\frac{d\mathsf{\Lambda }\left(f\right)}{df}.\end{array}$$

(62)

The Equation (61) describes a process that undergoes standard Fickian diffusion dominated by a nonlinear drift. The latter imposes the non-standard form given by Equation (56) to the stationary distribution function.

5.3. The Case of Kappa Kinetics

As a working example, let us briefly consider the $\kappa$-kinetics [1], where the generalized logarithm is given by

$$\begin{array}{c}\hfill \mathsf{\Lambda }\left(f\right)={ln}_{\kappa }\left(f\right),\end{array}$$

(63)

${ln}_{\kappa }\left(f\right)$ being the $\kappa$-logarithm defined by means of the equation

$$\begin{array}{c}\hfill {ln}_{\kappa }\left(f\right)=\frac{f^{\kappa }-f^{-\kappa }}{2\kappa }.\end{array}$$

(64)

The inverse of the logarithm function, i.e., the $\kappa$-exponential, is given by

$$\begin{array}{c}\hfill {exp}_{\kappa }\left(x\right)={\left(\sqrt{1+{\kappa }^2{x}^2}+\kappa x\right)}^{1/\kappa }.\end{array}$$

(65)

In this case, the $a\left(f\right)$ and $b\left(f\right)$ functions assume the expressions

$$\begin{array}{c}a\left(f\right)=\sqrt{\frac{2f}{f^{\kappa }+f^{-\kappa }}}exp\left[\frac{f^{\kappa }-f^{-\kappa }}{4\kappa }\right],\hfill \end{array}$$

(66)

$$\begin{array}{c}b\left(f\right)=\sqrt{\frac{2f}{f^{\kappa }+f^{-\kappa }}}exp\left[-\frac{f^{\kappa }-f^{-\kappa }}{4\kappa }\right],\hfill \end{array}$$

(67)

while the respective Fokker–Planck equation becomes

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{\partial }{\partial x}\left(D\beta \frac{\partial U}{\partial x}\frac{2f}{f^{\kappa }+f^{-\kappa }}+D\frac{\partial f}{\partial x}\right).\end{array}$$

(68)

It is remarkable that the latter equation is different from the one proposed in [1]. Both equations are correct and describe two different nonlinear kinetics that admit the same stationary solution

$$\begin{array}{c}\hfill f_s={exp}_{\kappa }\left[-\beta (U-\mu )\right].\end{array}$$

(69)

The main difference between the two kinetics is that the present Fokker–Planck equation describes a process which undergoes ordinary Fickian diffusion, second term on the right hand-side of (68), in the presence of nonlinear drift, first term on the right hand-side of (68), while the Fokker–Planck equation of Ref. [1] describes a process that undergoes non-Fickian diffusion in the presence of linear drift.

6. Conclusions

Let us summarize the main results of the present study. Starting from the Kinetic Interaction Principle introduced in [1], we obtained the evolution equation of a general statistical system defined on a one-dimensional lattice. The KIP-based master equation in the continuum limit yielded the already known nonlinear-Fokker–Planck Equation (2), which is widely used in the literature to describe anomalous diffusion in condensed matter physics and nonconventional statistical physics. Our derivation provides a better, physical understanding of the underlying many-body lattice dynamics at the microscopic level.

A second result is related to the possibility of directly using the KIP-based master equation as the numerical discretization scheme for the solution of the Fokker–Planck partial differential Equation (2). This approach resolves the ambiguity that results from different numerical discretization schemes of the Fokker–Planck equation. The proposed discretization scheme based on the master equation is physically motivated and follows from the KIP that describes microscopic nonlinear dynamics. On the contrary, discretization schemes based on finite differences follow from mathematically considerations. In addition, the master equations obtained from the nonlinear Fokker–Planck equation by applying different discretization schemes cannot be obtained starting from the KIP. The master equation proposed in this paper and obtained by means of the KIP is thus uniquely defined by microscopic interactions.

Finally, after noting that the Fokker–Planck current expresses the sum of two distinct contributions, i.e., the nonlinear drift current and the generalized Fick current as shown in Equations (51)–(53), we have demonstrated that some important nonlinear kinetics formulations proposed in the literature can be successfully described by ordinary Fickian diffusion. This can be accomplished by introducing a nonlinear drift term that is combined with Fickian diffusion. This combination leads to the same stationary solution of the Fokker–Planck equation as an equation that involves non-Fickian diffusion and linear drift.

It is straightforward to extend the KIP-based approach for the derivation of physically inspired master equations to lattices in two and three dimensions. For example, in the case of a two dimensional square lattice whose sites are labeled by the integer indices $i,j$, the master Equation (7) with nearest-neighbor transitions becomes

$$\begin{array}{cc}\hfill \frac{d{f}_{i,j}}{dt}=& w_{i-1,j}^{+}a\left(f_{i-1,j}\right)b\left(f_{i,j}\right)-w_{i,j}^{-}a\left(f_{i,j}\right)b\left(f_{i-1,j}\right)+w_{i+1,j}^{-}a\left(f_{i+1,j}\right)b\left(f_{i,j}\right)-w_{i,j}^{+}a\left(f_{i,j}\right)b\left(f_{i+1,j}\right)\hfill \\ \hfill & \hfill +r_{i,j-1}^{+}a\left(f_{i,j-1}\right)b\left(f_{i,j}\right)-r_{i,j}^{-}a\left(f_{i,j}\right)b\left(f_{i,j-1}\right)+r_{i,j+1}^{-}a\left(f_{i,j+1}\right)b\left(f_{i,j}\right)-r_{i,j}^{+}a\left(f_{i,j}\right)b\left(f_{i,j+1}\right).\end{array}$$