Entropy Production in a Fractal System with Diffusive Dynamics

Abstract

We study the entropy production in a fractal system composed of two subsystems, each of which is subjected to an external force. This is achieved by using the H-theorem on the nonlinear Fokker–Planck equations (NFEs) characterizing the diffusing dynamics of each subsystem. In particular, we write a general NFE in terms of Hausdorff derivatives to take into account the metric of each system. We have also investigated some solutions from the analytical and numerical point of view. We demonstrate that each subsystem affects the total entropy and how the diffusive process is anomalous when the fractal nature of the system is considered.

1. Introduction

Perhaps one of the greatest achievements of modern physics is the statistical mechanics formalism, first developed by Boltzmann [1] and later expanded by Gibbs [2], to describe the dynamics of microscopic particles and their connection with observable thermodynamics. Entropy is the crucial ingredient in this theory, which must verify the H-Theorem, establishing that nonequilibrium systems will reach equilibrium after long-time evolution, essentially posing a way of investigating the rule of additivity for systems with different entropies [3,4,5,6]. From a microscopic perspective, the Langevin equation is typically employed when studying brownian motion, which results in the time-dependent position of the particle when friction and noise are considered. From the phenomenological point of view, nonequilibrium statistical mechanics, including entropy calculation, is typically calculated with the Fokker–Planck equation which, in its linear form, results in Gaussian distributions, whereas the NFE is generally used to describe anomalous behavior often seen in long-range interaction [7], memory effects [8], porous media [9,10], and many others (see, for example, Refs. [11,12,13] and references therein), all of which display non-Gaussian distributions, typical of non-Markovian characteristics. A remarkable case in which non-Gaussian behavior may be present is in the case of interacting particles/systems. For example, a set of particles may react with another set, which may result in combination, or conversion of one chemical species into another. In this case, as the process occurs, the interaction dynamically affects the physical parameters of the system as a whole, including the entropy production and the diffusion coefficient, which depend on temperature and how the particles interact with each other [14]. Thus, such cases may be seen as a dynamically coupled system [14].

An anomalous process, such as anomalous diffusion, is a burgeoning field of research, as it is often observed across several research fields, from separation to biological media diffusion. From a theoretical point of view, fractional calculus is commonly used to describe anomalous diffusion processes. However, fractional calculus is nonlocal, making it troublesome in the numerical simulation of long-term and large-scale problems [15]. Also, the Mittag–Leffler decay [16] and the Lévy stable statistics [17] resulting from fractional calculus do not describe well some stretched relaxation and stretched Gaussian statistics [15]. An alternative approach, a local differential operator (conformable derivative [18]), is the Hausdorff derivative [19], which relates the fractal nature of space (or time) to the appearance of anomalous behavior such as observed in anomalous diffusion [15]. Fractal derivatives have been applied in several contexts, from anomalous diffusion to viscoelasticity and water transport [15]. The use of Hausdorff derivatives in generalized Fokker–Planck equations is also scarce. For example, recently [20], the solutions for the generalized Fokker–Planck equation with conformable and integro-differential operators were analyzed, indicating that the mean square displacement (MSD) for the studied case presented some difference when compared to the Caputo derivative [21], but presents the same time dependence as the scaled Brownian motion [20]. The Fokker–Planck equation of fractal curves was also explored in reference [22], where anomalous diffusion is observed. Hausdorff derivatives have also been used to describe anomalous transport in porous media by assuming a non-Euclidian fractal metric, displaying better agreement with experimental data for the heavy tail distribution [23].

This article studies the entropy production in a fractal system composed of two subsystems. Hence, this article is devoted to study a system composed of two different set of diffusing particles, each set a subsystem of the whole system. As each subsystem relaxes, the change in one subsystem affects the other, since the time dependent diffusion coefficient of one subsystem has to change as the distribution of the other subsystem evolve, as previously discussed for interacting particles. In fact, it is well known that morphology directly affects diffusion [24,25], so the change in concentration of interacting particles, commonly reported as molecular crowding [26], may be viewed as a morphology change. Aiming to make our model as general as possible, we use a set of NFEs to describe the evolution of each subsystem and analyze different dependences of the time-dependent diffusion coefficient of one subsystem on the distribution of the other. Also, we solve our equations with the Hausdorff metric so each subsystem may be viewed as fractal in nature. Furthermore, we use this system to calculate the entropy production, which is also generalized by assuming different forms for the entropy of each subsystem. The thermodynamic connection is achieved by using the H-Theorem on the NFE equations, characterizing the diffusing dynamics of each subsystem, each subjected to an external force. In this instance, the diffusion coefficient is temperature-dependent, and the entropy is also temperature-dependent, so the behavior of one system affects the other. Consequently, each subsystem serves as a thermal bath for the other, and the H-Theorem guarantees a connection between the subsystems. This feature also helps us to analyze entropy and the zero-law of the thermodynamics in terms of a relaxation process governed by nonlinear Fokker–Planck equations or the mixing between different regimes of diffusion. In fact, the nonlinear or mixing between different terms connected to the relaxation process directly influences the functional entropy form and, consequently, on the properties such as the additivity when different subsystems are added to compose the system. In this sense, the approach considered here gives a suitable form for the entropy and nonlinear Fokker–Planck equations in connection with these processes in a thermostatistic context for usual or anomalous relaxation processes. In this manner, the results will show that the dynamic of each one has a direct influence on the other as a thermal bath since the coupling appears in the diffusive term, which is related to how the system will spread. We express a general NFE regarding Hausdorff derivatives, considering each system’s fractal metric. We show that the fractal order significantly affects entropy production and leads to an anomalous diffusion behavior for each design.

2. Nonlinear Fokker–Planck Equations and Hausdorff Derivative

We begin by setting up the NFEs that explain the behavior of each part of a composite system:

(1)$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}{\rho }_1(x_1,t)& =& -\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\mathcal{J}}_1(x_1,t)\hfill \end{array}$

(2)$\begin{array}{ccc}\hfill {\mathcal{J}}_1(x_1,t)& =& -\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)+\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left[F_1\left(x_1\right){\rho }_1(x_1,t)\right]\hfill \end{array}$

which implies

(3)$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_1(x_1,t)=\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left[F_1\left(x_1\right){\rho }_1(x_1,t)\right]\right\}\end{array}$

and

(4)$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}{\rho }_2(x_2,t)& =& -\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\mathcal{J}}_2(x_2,t)\hfill \end{array}$

(5)$\begin{array}{ccc}\hfill {\mathcal{J}}_2(x_2,t)& =& -\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)+\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left[F_2\left(x_2\right){\rho }_2(x_2,t)\right]\hfill \end{array}$

which implies

(6)$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_2(x_2,t)=\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left[F_2\left(x_2\right){\rho }_2(x_2,t)\right]\right\}.\end{array}$

In Equations (1)–(6), ${\rho }_{1\left(2\right)}$ is the particle distribution in system $1\left(2\right)$. An external force $F_i\left(x_i\right)$ is applied to each subsystem, with i being 1 or 2. This force is connected to the potential energy ${\varphi }_i$ as $F_i=-{\partial }_{{\xi }_{\mathcal{I},i}\left(x_i\right)}{\varphi }_i\left(x_i\right)$, $x_1$ and $x_2$ are defined in the range $(-\infty ,\infty )$ where the diffusion proceeds, $\mathcal{D}$ is the diffusion coefficient, and the spatial operator is the Hausdorff derivative [15,27], as shall be defined below. Furthermore, $P_{1\left(2\right)}({\rho }_{1\left(2\right)},t)$ is a functional depending on the distribution of particles here used to generalize the problem. In fact, this research will make use of $P_1({\rho }_1,t)$ and $P_2({\rho }_2,t)$ to illustrate a certain phenomenon, as previously seen in porous media [28], anomalous diffusion [29], overdamped systems [30], and the Boltzmann equation with a correlation term [31]. Equations (3) and (6) also extend the equations used in Refs. [32,33,34,35,36,37] to analyze the H-theorem and the entropy production enable us to consider different contexts. One of them is the relaxation to an equilibrium, a system composed of subsystems that are governed by Equations (3) and (6), which may be connected to the zero law of the thermodynamics in generalized thermostatistics contexts [38,39,40]. The spatial differential operator in Equations (1)–(5), the Hausdorff derivative [15,27], is defined as follows:

(7)$\begin{array}{c}\hfill \frac{\partial }{\partial {\xi }_{\mathcal{I},i}\left(x_i\right)}{\rho }_i(x_i,t)=\underset{x_i^{\prime }\to x_i}{\lim }\frac{{\rho }_i(x_i,t)-{\rho }_i(x_i^{\prime },t)}{{\xi }_{\mathcal{I},i}\left(x_i\right)-{\xi }_{\mathcal{I},i}\left(x_i^{\prime }\right)}=\frac{1}{{\xi }_i\left(x_i\right)}\frac{\partial }{\partial x_i}{\rho }_i(x_i,t),\end{array}$

where ${\xi }_{\mathcal{I},i}\left(x_i\right)={\int }^{x_i}d\overline{x}{\xi }_i\left(\overline{x}\right)$ which, as previously discussed, may be connected to the fractal aspects of each system [19].

2.1. H-Theorem

We begin by applying the H-Theorem, taking into account $P_1({\rho }_1,t)$ formally equal to $P_2({\rho }_2,t)$. We will then explore the consequences of making a different selection and how it impacts the entropy of the combined system. Subsequently, we will calculate the Helmholtz free energy and its rate of change, as outlined in Refs. [35,41,42]. The free energy is expressed as $F=U-TS$, with the internal energy, U, given as:

(8)$\begin{array}{c}\hfill U={\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\Xi (x_1,x_2){\rho }_1(x_1,t){\rho }_2(x_2,t),\end{array}$

with $\Xi (x_1,x_2)={\varphi }_1\left(x_1\right)+{\varphi }_2\left(x_2\right)$, and the entropy S calculated as:

(9)$\begin{array}{c}\hfill \mathcal{S}=k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)s({\rho }_1,{\rho }_2)\end{array}$

where $s({\rho }_1,{\rho }_2)$ represents a generalized entropy form. For example, if $s({\rho }_1,{\rho }_2)$=${\rho }_1\ln {\rho }_1$+${\rho }_2\ln {\rho }_2$, we recover the classical Boltzmann entropy. The total free energy of the system composed of two subsystems is expressed by Equations (8) and (9).

(10)$\begin{array}{c}\hfill F={\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)[\Xi (x_1,x_2){\rho }_1(x_1,t){\rho }_2(x_2,t)-kTs({\rho }_1,{\rho }_2)].\end{array}$

Before examining the H-Theorem with Equation (10), we will assume that $P_1({\rho }_1,t)$ and $P_2({\rho }_2,t)$ have the same structure and that the entropy is a function of the product of the probability densities of each subsystem, i.e., $s({\rho }_1,{\rho }_2)=s\left({\rho }_1{\rho }_2\right)$. This allows us to demonstrate that

(11)$\begin{array}{c}\hfill \frac{d}{dt}F={\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\Xi (x_1,x_2)-kT\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right]\frac{\partial }{\partial t}\left[{\rho }_1(x_1,t){\rho }_2(x_2,t)\right],\end{array}$

in which ${\rho }_{12}={\rho }_1{\rho }_2$. By performing some calculations and using Equations (2) and (5), we have that

(12)$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\Xi (x_1,x_2){\rho }_2-kT{\rho }_2\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left[F_1\left(x_1\right){\rho }_1(x_1,t)\right]\right\}\hfill \\ & +& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\Xi (x_1,x_2){\rho }_1-kT{\rho }_1\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left[F_2\left(x_2\right){\rho }_2(x_2,t)\right]\right\}.\hfill \end{array}$

We now assume the following conditions: ${\rho }_i(x_i\to \pm \infty ,t)\to 0$ and ${\partial }_{{\xi }_{\mathcal{I},i}\left(x_i\right)}{\rho }_i(x_i\to \pm \infty ,t)$→ 0. Thus, Equation (12) becomes, after integration by parts:

(13)$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& -{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_2-kT{\rho }_2^2\frac{\partial {\rho }_1}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)+\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1(x_1,t)\right\}\hfill \\ & -& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1-kT{\rho }_1^2\frac{\partial {\rho }_2}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)+\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_2(x_2,t)\right\}\hfill \end{array}$

(for more details, see the Appendix A). Equation (13) is a general result and can be connected to different relaxation processes depending on the choice of $P_1({\rho }_1,t)$ and $P_2({\rho }_2,t)$ present in Equations (3) and (6). It can also be used to evidence the interaction between the subsystems in connection with the relaxation process each subsystem exhibits. This feature will be evident below with the analysis of the condition required for the H-Theorem and the previous assumption for the entropy, i.e., $s({\rho }_1,{\rho }_2)=s\left({\rho }_1{\rho }_2\right)$. Also, this result maintains the additivity in Penrose sense [43], i.e., $S\left({\rho }_{12}\right)=S\left({\rho }_1{\rho }_2\right)$ required for a system composed of independent subsystems when the standard entropy is employed. Thus, the conditions required by Equation (13) to verify the H-theorem will define a suitable entropy for the relaxation process described in terms of Equations (3) and (6) for a choice of $P_1({\rho }_1,t)$ and $P_2({\rho }_2,t)$, and consequently, the properties of this entropy in connection with a thermostatistical context.

To proceed with our analysis, we consider that

(14)$\begin{array}{c}\hfill P_i({\rho }_i,t)={\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\mathcal{D}}_{j,\gamma }\left(t\right){\rho }_i^{\overline{\gamma }}(x_i,t),\end{array}$

with $j\ne i$, $p\left(\overline{\gamma }\right)$ is a distribution, and ${\mathcal{D}}_{j,\overline{\gamma }}\left(t\right)={\int }_{-\infty }^{\infty }d{x}_j{\xi }_j\left(x_j\right){\rho }_j^{\overline{\gamma }}(x_j,t)$. Note that the distribution $p\left(\overline{\gamma }\right)$ is connected with the nonlinear term present in Equations (3) and (6), which leads us a diffusive term with different diffusion regimes. In connection with the porous media equation and the Tsallis framework [11], $\overline{\gamma }$ may be related to the nonlinearity present in Equations (3) and (6) after substituting Equation (14) and also with the extension of the entropy to accommodate a thermostatistic context. The range of $\gamma$ is connected with the choice performed for the distribution $p\left(\overline{\gamma }\right)$, which defines how will be the behavior of the diffusive term. One possibility will be analyzed later by considering two different regimes of diffusion, with $p\left(\overline{\gamma }\right)=\delta (\overline{\gamma }-1)/2+\delta (\overline{\gamma }-\nu )/2$, with $\max (1,\nu )\le \gamma$. The choice performed for ${\mathcal{D}}_{j,\overline{\gamma }}\left(t\right)$ also implies that each subsystem influences the other, i.e., Equations (2) and (5) are coupled by the diffusive term. Thus, ${\mathcal{D}}_{j,\overline{\gamma }}\left(t\right)$ introduces interactions between the subsystems during the thermalization process, where each subsystem works as an additional thermal bath to the other. By substituting Equation (14) into Equation (13), we have that

(15)$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& -{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\frac{1}{{\rho }_1\left(x_1\right)}\hfill \\ & \times & \left\{{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_2{\rho }_1-kT{\rho }_2^2{\rho }_1\frac{\partial {\rho }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right]\right.\hfill \\ & \times & \left.{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1{\rho }_2+\mathcal{D}\frac{\partial {\rho }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_2\frac{\partial }{\partial {\rho }_{12}}\left({\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\rho }_2^{\overline{\gamma }}{\rho }_1^{\overline{\gamma }}\right)\right]\right\}\hfill \\ & -& {\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\frac{1}{{\rho }_2\left(x_2\right)}\hfill \\ & \times & \left\{{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\left[\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1{\rho }_2-kT{\rho }_1^2{\rho }_2\frac{\partial {\rho }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right]\right.\hfill \\ & \times & {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\left.\left[\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1{\rho }_2+\mathcal{D}\frac{\partial {\rho }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1\frac{\partial }{\partial {\rho }_{12}}\left({\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\rho }_2^{\overline{\gamma }}{\rho }_1^{\overline{\gamma }}\right)\right]\right\}.\hfill \end{array}$

We verify that

(16)$\begin{array}{ccc}\hfill \frac{d}{dt}F\le 0\mathrm{if}-kT{\rho }_j^2{\rho }_i\frac{{\partial }^2}{\partial {\rho }_{ij}^2}s\left({\rho }_{ij}\right)& =& \mathcal{D}{\rho }_j\frac{\partial }{\partial {\rho }_{ij}}\left({\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\rho }_i^{\overline{\gamma }}{\rho }_j^{\overline{\gamma }}\right),\hfill \end{array}$

(17)$\begin{array}{ccc}& =& \mathcal{D}{\rho }_j\frac{\partial }{\partial {\rho }_{ij}}\left({\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\rho }_{ij}^{\overline{\gamma }}\right),\hfill \end{array}$

for $i=1,2$ and $j=1,2$ with $i\ne j$, $\mathcal{D}=kT$, and ${\rho }_{ij}={\rho }_i{\rho }_j$, which implies

(18)$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& -{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\frac{1}{{\rho }_1\left(x_1\right)}\hfill \\ & \times & {\left\{{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_2{\rho }_1-kT{\rho }_2^2{\rho }_1\frac{\partial {\rho }_1}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right]\right\}}^2\hfill \\ & -& {\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\frac{1}{{\rho }_2\left(x_2\right)}\hfill \\ & \times & {\left\{{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\left[\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1{\rho }_2-kT{\rho }_1^2{\rho }_2\frac{\partial {\rho }_2}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right]\right\}}^2.\hfill \end{array}$

By solving Equation (16) under the conditions defined in Refs. [34,35,41,42], we obtain

(19)$\begin{array}{c}\hfill s\left({\rho }_{12}\right)={\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right)\frac{1}{\overline{\gamma }-1}\left({\rho }_{12}-{\rho }_{12}^{\overline{\gamma }}\right).\end{array}$

The entropy for the composite system is given by

(20)$\begin{array}{c}\hfill \mathcal{S}=k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right)\frac{1}{\overline{\gamma }-1}\left({\rho }_{12}-{\rho }_{12}^{\overline{\gamma }}\right),\end{array}$

which can also be rewritten as

(21)$\begin{array}{c}\hfill \mathcal{S}=k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right)\frac{1}{\overline{\gamma }-1}\left({\rho }_1{\rho }_2-{\left({\rho }_1{\rho }_2\right)}^{\overline{\gamma }}\right)\end{array}$

and, consequently, as

(22)$\begin{array}{c}\hfill \mathcal{S}=k{\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right)\frac{1}{\overline{\gamma }-1}\left(1-{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\left[{\rho }_1\left(x_1\right){\rho }_2\left(x_2\right)\right]}^{\overline{\gamma }}\right).\end{array}$

Equation (22) may have several particular cases, such as the Tsallis and Kaniadakis entropies, depending on the choice of $p\left(\overline{\gamma }\right)$. It is also worth mentioning that the Boltzmann–Gibbs entropy is recovered from Equation (22), which is connected with the linear Fokker–Planck equation as discussed in Ref. [44]. In addition, different from the standard situation, to satisfy the H-theorem and preserve the form $s=s\left({\rho }_1{\rho }_2\right)$, the nonlinear Fokker–Planck equations imply that the systems interact with each other. Equation (22) can be connected to the unusual additivity present in the Tsallis formalism as follows:

(23)$\begin{array}{c}\hfill S={\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right)(S_1\left(\overline{\gamma }\right)+S_2\left(\overline{\gamma }\right)+[(1-\overline{\gamma })/k]S_1\left(\overline{\gamma }\right)S_2\left(\overline{\gamma }\right)),\end{array}$

where

(24)$\begin{array}{c}\hfill {\mathcal{S}}_{1\left(2\right)}\left(\overline{\gamma }\right)=\frac{k}{\overline{\gamma }-1}\left(1-{\int }_{-\infty }^{\infty }d{x}_{1\left(2\right)}{\xi }_{1\left(2\right)}\left(x_{1\left(2\right)}\right){\left[{\rho }_{1\left(2\right)}\left(x_{1\left(2\right)}\right)\right]}^{\overline{\gamma }}\right).\end{array}$

This result, connected to each system’s relaxation, shows how each part’s entropy is added to compose the total entropy. Equation (23) has been applied in several situations such as black hole [45], inanimated and living matter [46], and interacting particles [47].

The NFE that emerges from the previous analysis, i.e., from Equation (14), which is related to the nonlinear term and verifies the H-theorem, for each subsystem can be written as follows:

(25)$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_1(x_1,t)=\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left\{{\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\overline{\mathcal{D}}}_{1,\overline{\gamma }}\left(t\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1^{\overline{\gamma }}(x_1,t)-F_1\left(x_1\right){\rho }_1(x_1,t)\right\}\end{array}$

and

(26)$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_2(x_2,t)=\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left\{{\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\overline{\mathcal{D}}}_{2,\overline{\gamma }}\left(t\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_2^{\overline{\gamma }}(x_2,t)-F_2\left(x_2\right){\rho }_2(x_2,t)\right\},\end{array}$

with ${\overline{\mathcal{D}}}_{i,\overline{\gamma }}\left(t\right)={\mathcal{D}}_{i,\gamma }\left(t\right)\mathcal{D}$, which demonstrates the influence of one of the subsystems on the other. In addition, the components that constitute the diffusive element can also be associated with anomalous diffusion processes with distinct diffusion regimes. This feature was investigated, for example, in Ref. [48] by considering the Fokker–Planck equation with a different form for the diffusive term or Langevin equations with additive noises in Ref. [49]. It is also worth mentioning that the H-theorem shows us a suitable way for the dynamic processes in connection with the system’s entropy.

We will now contemplate a general situation in which the diffusion terms have a distinct nonlinear relationship with the distributions. This implies that the systems have distinct dynamic characteristics regulated by the nonlinear connection with the distribution found in the diffusive term. Using the preceding equations and having in mind Equation (10), we may write

(27)$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\Xi (x_1,x_2)\frac{\partial }{\partial t}\left[{\rho }_1(x_1,t){\rho }_2(x_2,t)\right]\right.\hfill \\ & -& \left.kT\left[\frac{\partial }{\partial {\rho }_1}s({\rho }_1,{\rho }_2)\frac{\partial }{\partial t}{\rho }_1(x_1,t)+\frac{\partial }{\partial {\rho }_2}s({\rho }_1,{\rho }_2)\frac{\partial }{\partial t}{\rho }_2(x_2,t)\right]\right\},\hfill \end{array}$

which implies

(28)$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\Xi (x_1,x_2){\rho }_2-kT\frac{\partial }{\partial {\rho }_1}s({\rho }_1,{\rho }_2)\right\}\hfill \\ & \times & \frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)+\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1(x_1,t)\right\}\hfill \\ & +& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\Xi (x_1,x_2){\rho }_1-kT\frac{\partial }{\partial {\rho }_2}s({\rho }_1,{\rho }_2)\right\}\hfill \\ & \times & \frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)+\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_2(x_2,t)\right\}.\hfill \end{array}$

After some calculations, it is possible to show that

(29)$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& -{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\frac{1}{{\rho }_1\left(x_1\right)}{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\hfill \\ & \times & \left\{\frac{\partial \varphi \left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_2\left(x_2\right){\rho }_1\left(x_1\right)-kT\left[{\rho }_1\frac{{\partial }^2}{\partial {\rho }_1^2}s({\rho }_1,{\rho }_2)\right]\frac{\partial {\rho }_1}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\right\}\hfill \\ & \times & \left\{\left[\mathcal{D}\frac{\partial }{\partial {\rho }_1}{P}_1({\rho }_1,t)\right]\frac{\partial {\rho }_1}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}+\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1(x_1,t)\right\}\hfill \\ & -& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\hfill \\ & \times & \frac{1}{{\rho }_2\left(x_2\right)}\left\{\frac{\partial \varphi \left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1\left(x_1\right){\rho }_2\left(x_2\right)-kT\left[{\rho }_2\frac{{\partial }^2}{\partial {\rho }_2^2}s({\rho }_1,{\rho }_2)\right]\frac{\partial {\rho }_2}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\right\}\hfill \\ & \times & \frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left\{\left[\mathcal{D}\frac{\partial }{\partial {\rho }_2}{P}_2({\rho }_2,t)\right]\frac{\partial {\rho }_2}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}+\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_2(x_2,t)\right\}.\hfill \end{array}$

Now, we assume, for example, the case

(30)$\begin{array}{c}\hfill P_1({\rho }_1,t)={\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\mathcal{D}}_{2,\nu }\left(t\right){\rho }_1^{\overline{\gamma }}(x_1,t)\end{array}$

and

(31)$\begin{array}{c}\hfill P_2({\rho }_2,t)={\int }_0^{\nu }d\overline{\nu }p\left(\overline{\nu }\right){\mathcal{D}}_{1,\overline{\gamma }}\left(t\right){\rho }_2^{\overline{\nu }}(x_2,t),\end{array}$

with

(32)$\begin{array}{c}\hfill {\mathcal{D}}_{2,\nu }\left(t\right)={\int }_0^{\nu }d\overline{\nu }p\left(\overline{\nu }\right)\frac{1}{\overline{\nu }-1}{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\rho }_2^{\overline{\nu }}(x_2,t)\end{array}$

and

(33)$\begin{array}{c}\hfill {\mathcal{D}}_{1,\gamma }\left(t\right)={\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right)\frac{1}{\overline{\gamma }-1}{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\rho }_1^{\overline{\gamma }}(x_1,t).\end{array}$

This implies that Equations (2) and (5) have different forms, and thus, the two subsystems have distinct relaxation processes. In Ref. [14], a particular case was studied by looking at the interaction between the two subsystems. Each choice has its implications for the total entropy of the composite system. By examining Equations (32) and (33), we can see that the entropy must satisfy the following equations:

(34)$\begin{array}{c}\hfill -{\rho }_1\frac{{\partial }^2}{\partial {\rho }_1^2}s({\rho }_1,{\rho }_2)={\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\int }_0^{\nu }d\overline{\nu }p\left(\overline{\nu }\right)\frac{\overline{\gamma }}{\overline{\nu }-1}{\rho }_2^{\overline{\nu }}{\rho }_1^{\overline{\gamma }-1}\end{array}$

and

(35)$\begin{array}{c}\hfill -{\rho }_2\frac{{\partial }^2}{\partial {\rho }_2^2}s({\rho }_1,{\rho }_2)={\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\int }_0^{\nu }d\overline{\nu }p\left(\overline{\nu }\right)\frac{\overline{\nu }}{\overline{\gamma }-1}{\rho }_2^{\overline{\nu }-1}{\rho }_1^{\overline{\gamma }}\end{array}$

to verify

(36)$\begin{array}{c}\hfill \frac{d}{dt}F\le 0,\end{array}$

and, consequently, to satisfy the H-Theorem. A solution for the previous system of equations is

(37)$\begin{array}{c}\hfill s({\rho }_1,{\rho }_2)={\int }_0^{\gamma }d\overline{\gamma }\frac{p\left(\overline{\gamma }\right)}{\overline{\gamma }-1}{\int }_0^{\nu }d\overline{\nu }\frac{p\left(\overline{\nu }\right)}{\overline{\gamma }-1}\left({\rho }_1{\rho }_2-{\rho }_2^{\overline{\nu }}{\rho }_1^{\overline{\gamma }}\right).\end{array}$

This result allows us to write the total entropy of this system as follows:

(38)$\begin{array}{c}\hfill \mathcal{S}=k{\int }_0^{\gamma }d\overline{\gamma }\frac{p\left(\overline{\gamma }\right)}{\overline{\gamma }-1}{\int }_0^{\nu }d\overline{\nu }\frac{p\left(\overline{\nu }\right)}{\overline{\nu }-1}\left[1-{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\rho }_2^{\overline{\nu }}(x_2,t){\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\rho }_1^{\overline{\gamma }}(x_1,t)\right].\end{array}$

This result for the entropy is distinct from the one given by Equation (22), which was derived from a different selection of NFEs. It is the consequence of combining subsystems with distinct relaxation processes, each of which has its entropy. Additionally, Equation (38) is linked to the combination of Tsallis entropies with different $q-$ indices [50,51,52]. The solution can be obtained in this context by using $q-$ exponential functions.

Let us now consider a particular example of the results mentioned above when different behaviors for $P_1({\rho }_1,t)$ and $P_2({\rho }_2,t)$ are chosen. We consider the case obtained from Equations (30) and (31) for $p\left(\overline{\gamma }\right)=\delta (\overline{\gamma }-\gamma )$ and $p\left(\overline{\nu }\right)=\delta (\overline{\nu }-\nu )$ in absence of external force, yielding

(39)$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_1=\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left[{\mathcal{D}}_{2,\nu }\left(t\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1^{\gamma }(x_1,t)\right]\end{array}$

and

(40)$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_2=\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left[{\mathcal{D}}_{1,\gamma }\left(t\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_2^{\nu }(x_2,t)\right].\end{array}$

Hence, it is simple to verify that

(41)$\begin{array}{c}\hfill {\rho }_1(x_1,t)={\exp }_{\gamma }\left[-{\beta }_1\left(t\right){\xi }_{\mathcal{I},1}^2\left(x_1\right)\right]/{\mathcal{Z}}_1\left(t\right)\end{array}$

and

(42)$\begin{array}{c}\hfill {\rho }_2(x_1,t)={\exp }_{\nu }\left[-{\beta }_2\left(t\right){\xi }_{\mathcal{I},2}^2\left(x_2\right)\right]/{\mathcal{Z}}_2\left(t\right),\end{array}$

are solutions for Equations (39) and (40). The previous solutions for the NFEs can be also obtained from the maximum principle of entropy when the entropy

(43)$\begin{array}{c}\hfill S_{1\left(2\right)}=\frac{k}{\gamma \left(\nu \right)-1}\left(1-\int d{x}_{1\left(2\right)}{\xi }_{1\left(2\right)}\left(x_{1\left(2\right)}\right){\left[{\rho }_{1\left(2\right)}(x_{1\left(2\right)},t)\right]}^{\gamma \left(\nu \right)}\right)\end{array}$

with the constraints

(44)$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }d{x}_{1\left(2\right)}{\xi }_{1\left(2\right)}\left(x_{1\left(2\right)}\right){\rho }_{1\left(2\right)}^{\gamma \left(\nu \right)}(x_{1\left(2\right)},t){\left[{\xi }_{\mathcal{I},1\left(2\right)}\left(x_{1\left(2\right)}\right)\right]}^2={\sigma }_{{\xi }_{\mathcal{I}}\left(x\right),1\left(2\right)}^2\left(t\right)\end{array}$

and

(45)$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }d{x}_{1\left(2\right)}{\xi }_{1\left(2\right)}\left(x_{1\left(2\right)}\right){\rho }_{1\left(2\right)}(x_{1\left(2\right)},t)=1.\end{array}$

A particular situation was worked out in Ref. [11] for the porous media equation. From Equations (41) and (42), we can verify how the distribution is affected by the space metric and the parameter $\gamma \left(\nu \right)$. Figure 1 shows ${\rho }_{1\left(2\right)}{\mathcal{Z}}_{1\left(2\right)}$ vs. ${\beta }_{1\left(2\right)}x$ for ${\xi }_{\mathcal{I},1\left(2\right)}\left(x_{1\left(2\right)}\right)={\left|x_{1\left(2\right)}\right|}^{{\alpha }_{1\left(2\right)}}$. If ${\alpha }_{1\left(2\right)}=\gamma \left(\nu \right)=1$, as expected, we recover the usual behavior given by a Gaussian distribution. On the other hand, one can clearly see from figure Figure 1 that if the parameter coming from the Hausdorff derivative (${\alpha }_{1\left(2\right)}$) or the parameter $\gamma$ (which describes how the distribution of one kind of particle affects the time-dependent diffusion coefficient of the other, or, in different words, how molecular crowding affects diffusion) the behavior becomes non-Gaussian. Thus, depending on the choice of parameters $\nu$, $\gamma$, and ${\alpha }_{1\left(2\right)}$, the distribution may present different behaviors that are asymptotically characterized by short- or long-tailed distributions. In the latter case, it is possible to connect the results with the Lévy distributions. In addition, a particular scenario of Equations (39) and (40) has been worked out in Ref. [53] in connection with the fractal dimensions. Similar situations are found in the context of the Tsallis statistics, which are described by power-law distributions. However, small changes in the parameters ${\alpha }_{1\left(2\right)}$ and $\gamma \left(\nu \right)$ strongly affect the distribution of particles.

With ${\beta }_1\left(t\right)$, ${\beta }_2\left(t\right)$, ${\mathcal{Z}}_1\left(t\right)$, and ${\mathcal{Z}}_2\left(t\right)$ obtained from the following set of equations:

(46)$\begin{array}{ccc}\hfill \frac{1}{2{\beta }_1}\frac{d}{dt}{\beta }_1& =& -\frac{2\gamma }{\nu -1}\frac{{\mathcal{I}}_{\nu }}{{\mathcal{N}}_{\nu }^{\nu }{\mathcal{N}}_{\gamma }^{\gamma -1}}{\beta }_2^{(\nu -1)/2}{\beta }_1^{(\gamma +1)/2}\hfill \end{array}$

and

(47)$\begin{array}{ccc}\hfill \frac{1}{2{\beta }_2}\frac{d}{dt}{\beta }_2& =& -\frac{2\nu }{\gamma -1}\frac{{\mathcal{I}}_{\gamma }}{{\mathcal{N}}_{\gamma }^{\gamma }{\mathcal{N}}_{\nu }^{\nu -1}}{\beta }_2^{(\nu +1)/2}{\beta }_1^{(\gamma -1)/2}\hfill \end{array}$

with

(48)$\begin{array}{c}\hfill {\mathcal{I}}_{\kappa }=\left\{\begin{array}{cc}\frac{\mathrm{\Gamma }\left(\frac{1}{2}\right)\mathrm{\Gamma }\left(1+\frac{\kappa }{\kappa -1}\right)}{\sqrt{\kappa -1}\mathrm{\Gamma }\left(\frac{3}{2}+\frac{\kappa }{\kappa -1}\right)}& 1\le \kappa <2\\ \frac{\mathrm{\Gamma }\left(\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{\kappa }{1-\kappa }-\frac{1}{2}\right)}{\sqrt{1-\kappa }\mathrm{\Gamma }\left(\frac{\kappa }{1-\kappa }\right)}& 0\le \kappa \le 1\end{array}\right.,{\mathcal{N}}_{\kappa }=\left\{\begin{array}{cc}\frac{\mathrm{\Gamma }\left(\frac{1}{2}\right)\mathrm{\Gamma }\left(1+\frac{1}{\kappa -1}\right)}{\sqrt{\kappa -1}\mathrm{\Gamma }\left(\frac{3}{2}+\frac{1}{\kappa -1}\right)}& 1\le \kappa <2\\ \frac{\mathrm{\Gamma }\left(\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{1}{1-\kappa }-\frac{1}{2}\right)}{\sqrt{1-\kappa }\mathrm{\Gamma }\left(\frac{1}{1-\kappa }\right)}& 0\le \kappa \le 1\end{array}\right.\end{array}$

where $\kappa =\gamma$ or $\nu$ and ${\mathcal{Z}}_{1\left(2\right)}\left(t\right)\sqrt{{\beta }_{1\left(2\right)}\left(t\right)}={\mathcal{Z}}_{1\left(2\right)}\left(0\right)\sqrt{{\beta }_{1\left(2\right)}\left(0\right)}=constant$. Note that depending on the choice of the parameters $\gamma$ and $\nu$, Equations (41) and (42) may present a compact or a long-tailed behavior. In the last case, it is possible to connect the solutions with the Lévy distributions as performed in Refs. [54,55]. The solutions for Equations (46) and (47) can be found, and they are given by

(49)$\begin{array}{c}\hfill {\beta }_1\left(t\right)={\left[\frac{\nu +\gamma }{\nu -1}\frac{\gamma {\mathcal{I}}_{\nu }}{{\mathcal{N}}_{\nu }{\mathcal{N}}_{\gamma }^{\gamma -1}}{\mathcal{C}}^{\frac{\nu -1}{2}}t\right]}^{-\frac{2}{\nu +\gamma }}\end{array}$

and ${\beta }_2\left(t\right)=\mathcal{C}{\beta }_1\left(t\right)$ with

(50)$\begin{array}{c}\hfill \mathcal{C}=\frac{\gamma (\gamma -1)}{\nu (\nu -1)}\left(\frac{{\mathcal{I}}_{\nu }{\mathcal{N}}_{\gamma }}{{\mathcal{I}}_{\gamma }{\mathcal{N}}_{\nu }}\right).\end{array}$

By using the previous equations, it is possible to obtain the mean square displacement. Let us now consider the case of ${\xi }_{1\left(2\right)}\left(x\right)={\left|x\right|}^{{\alpha }_{1\left(2\right)}}$, which represents a fractal metric as proposed by Chen [19], and the definition that follows:

(51)$\begin{array}{c}\hfill {\sigma }_{x,1\left(2\right)}^2\left(t\right)=⟨{\left(x_{1\left(2\right)}-{⟨x_{1\left(2\right)}⟩}_{1\left(2\right)}\right)}^2⟩.\end{array}$

After performing some calculations, it is possible to show that

(52)$\begin{array}{c}\hfill {\sigma }_{x,1\left(2\right)}^2\left(t\right)=\frac{{\zeta }_{\gamma \left(\nu \right)}}{{\mathcal{N}}_{\gamma \left(\nu \right)}}\frac{1}{{\left[{\beta }_{1\left(2\right)}\left(t\right)\right]}^{\frac{1}{{\alpha }_{1\left(2\right)}}}},\end{array}$

where ${\zeta }_{\kappa }={\int }_0^{\infty }du{u}^{1/{\alpha }_{1\left(2\right)}-1/2}{\exp }_{\kappa }[-u]$, which implies in ${\sigma }_{x,1\left(2\right)}^2\left(t\right)\propto t^{2/\left[{\alpha }_{1\left(2\right)}(\gamma +\nu )\right]}$. Figure 2 illustrates the behavior of the mean square displacement, where the regions for the sub- and superdiffusion are shown. The black line corresponds to the usual diffusion. The choice of ${\alpha }_1$ and ${\alpha }_2$ connected to the fractal derivative directly influences the diffusion process. Another point concerning the behavior showed in Figure 2 is the influence of the dynamic of each system on the other since the results for the mean square displacement depend on the parameters $\gamma$ and $\nu$ connected to the nonlinearity present in the diffusive term. These feature connected to Figure 2 shows that the anomalous behavior is directly connected with the fractal metric of space [19] and has also been worked out in Refs. [56,57,58,59] by considering diffusion on fractal objects.

It is possible to demonstrate other scenarios by performing numerical computations, that is, by numerically solving Equations (2) and (5). Figure 3 and Figure 4 illustrate the case for which $P_{1\left(2\right)}({\rho }_{1\left(2\right)},t)$ is given by Equation (14) with $p\left(\overline{\gamma }\right)=\delta (\overline{\gamma }-1)/2+\delta (\nu -\overline{\gamma })/2$, $F_1(x_1,t)=-k_1(\theta \left(x_1\right)-\theta (-x_1)){\left|x_1\right|}^{{\alpha }_1}$ ($\theta \left(x\right)$ is the Heaviside function), and $F_2(x_1,t)=0$. Notice that the initial moments, with centered distributions, represent $t={10}^{-3}$, while the more spread distributions represent $t=1$. Remarkably, a simple change in the fractal value ${\alpha }_{1\left(2\right)}$ results in a large change in the observed diffusion regime of the particles. Numerical computations were performed for values $\nu$ greater than and less than one. The system was set in the range of −5000 to 5000, with a step size of $dx=0.02$ and a time step of $dt=0.000001$ to generate the results shown in the figures. The choices of $dx$ and $dt$ meet the requirement of $Ddt/\left(d{x}^2\right)<1/2$ for the stability of the solutions when the initial condition evolves in time to meet the boundary conditions [60,61]. Figure 3 show the dynamical behavior of the distributions ${\rho }_{1\left(2\right)}(x,t)$ for $\nu =1.2$, and two values of the fractal exponent ${\alpha }_{1\left(2\right)}$, which produce a non-Gaussian result. Note that the behavior presented in Figure 3 results from combining two different diffusive terms, one linear and the other nonlinear, besides the unusual metric for the space. Figure 4 demonstrates the anomalous nature of diffusive behavior when the fractal exponent changes from 0.9 to 1.2 for $\nu =0.95$ and $\nu =1.2$. It is also interesting to mention that the results presented in this figure show that the system under the influence of the external force has a different behavior for the mean square displacement, that is, the influence of the external force limits diffusion. On the other hand, the system without external force can spread freely. The influence is also present on the time dependence of the diffusion coefficient, which depends on the integral of the distribution in a nonlinear power law of the other distribution.

2.2. Entropy Production

We can examine the entropy production associated with Equation (21) by looking at the dynamics of ${\rho }_1(x_1,t)$ and ${\rho }_2(x_2,t)$ given by Equations (25) and (26). Differentiating Equation (21) with respect to time gives us the result:

(53)$\begin{array}{ccc}\hfill \frac{d}{dt}\mathcal{S}\left(t\right)& =& k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right]\frac{\partial }{\partial t}\left[{\rho }_1(x_1,t){\rho }_2(x_2,t)\right]\hfill \\ & =& -k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\rho }_2\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\mathcal{J}}_1(x_1,t)\hfill \\ & -& k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\rho }_1{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\mathcal{J}}_2(x_2,t)\hfill \end{array}$

Consequently, integration by parts is performed with the conditions that ${\mathcal{J}}_1(x_1\to \pm \infty ,t)$ and ${\mathcal{J}}_2(x_2\to \pm \infty ,t)$ both approaching zero:

(54)$\begin{array}{ccc}\hfill \frac{d}{dt}\mathcal{S}\left(t\right)& =& k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[{\rho }_2^2\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1\right]{\mathcal{J}}_1(x_1,t)\hfill \\ & +& k{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[{\rho }_1^2\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_2\right]{\mathcal{J}}_2(x_2,t).\hfill \end{array}$

It is feasible to make Equation (53) simpler by utilizing the equations from the H-Theorem:

(55)$\begin{array}{c}\hfill -kT{\rho }_1{\rho }_2^2\frac{\partial {\rho }_1}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)=\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)\end{array}$

and

(56)$\begin{array}{c}\hfill -kT{\rho }_2{\rho }_1^2\frac{\partial {\rho }_2}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)=\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t),\end{array}$

to obtain

(57)$\begin{array}{ccc}\hfill \frac{d}{dt}\mathcal{S}\left(t\right)& =& -\frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)F_1\left(x_1\right){\mathcal{J}}_1(x_1,t)-\frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)F_2\left(x_2\right){\mathcal{J}}_2(x_1,t)\hfill \\ & +& \frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\frac{{\mathcal{J}}_1^2(x_1,t)}{{\rho }_1(x_1,t)}+\frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\frac{{\mathcal{J}}_2^2(x_2,t)}{{\rho }_2(x_2,t)},\hfill \end{array}$

where ${\mathcal{J}}_1(x_1,t)$ and ${\mathcal{J}}_2(x_2,t)$ are defined by Equations (1) and (5), with $P_1({\rho }_1,t)$ and $P_2({\rho }_2,t)$ given by Equations (30) and (31). Equation (53) can also be written as follows:

(58)$\begin{array}{c}\hfill \frac{d}{dt}\mathcal{S}=\mathrm{\Pi }-\mathrm{\Phi }.\end{array}$

The entropy that is exchanged between the two subsystems ${\rho }_1$ and ${\rho }_2$ and their environment is called the flux of entropy. This can be expressed as follows:

(59)$\begin{array}{c}\hfill \mathrm{\Phi }=\frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)F_1\left(x_1\right){\mathcal{J}}_1(x_1,t)+\frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_1\right)F_2\left(x_2\right){\mathcal{J}}_2(x_1,t),\end{array}$

and the entropy-production term:

(60)$\begin{array}{c}\hfill \mathrm{\Pi }=\frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\frac{{\mathcal{J}}_1^2(x_1,t)}{{\rho }_1(x_1,t)}+\frac{1}{T}{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\frac{{\mathcal{J}}_2^2(x_2,t)}{{\rho }_2(x_2,t)}.\end{array}$

Since T and ${\rho }_i(x_i,t)$ are both positive, the desired result is obtained: $\mathrm{\Pi }\ge 0$. This result for the entropy production, given by Equation (57) and, thus, Equation (58) can also be confirmed for any entropy condition.

It is also important to mention that the entropy production in this framework for the external force considered here has the same behavior for the different choices of ${\xi }_{\mathcal{I},1\left(2\right)}\left(x_{1\left(2\right)}\right)$. This feature is directly connected with the definition of the integral used to obtain the entropy. Figure 5 and Figure 6 illustrate the behavior of entropy $\left(\mathcal{S}\right)$ and entropy production $\left(\dot{\mathcal{S}}\right)$. Figure 5 presents $\left(\mathcal{S}\right)$ and $\left(\dot{\mathcal{S}}\right)$ for different values of $\nu$. Figure 6 illustrates the behavior $\left(\mathcal{S}\right)$ and $\left(\dot{\mathcal{S}}\right)$ for the values ${\alpha }_{1\left(2\right)}=1.2$ and $\nu =1.2$ taking into account different external forces. Note that in each case, we consider a diffusion process with two different diffusive terms, one linear and the other nonlinear. Thus, the entropy evaluated for these cases results from the combination of different diffusive regimes. For $\nu <1$, the system is essentially governed by long-tailed distributions, and for $\nu >1$, short-tailed distributions govern the system. These features have a direct influence on entropy and entropy production, as shown in these figures. The entropy of the system is a growing function of time that eventually reaches a plateau for a long time, as expected. Interestingly, the external force, as shown in Figure 6, modifies the system’s entropy. In fact, entropy and entropy production depend on the external force acting on the systems, that is, how it can confine the system or not, during the particle spreading process.

3. Conclusions

We have investigated entropy production in a fractal system composed of two subsystems, each subject to an external force. This is achieved by using the H-theorem on the nonlinear Fokker–Planck equations (NFEs), characterizing the diffusing dynamics of each subsystem. To consider the metric of space in which the systems are embedded and, hence, the fractal nature that leads to anomalous diffusion, we expressed the general NFE in terms of Hausdorff derivatives. We investigated some solutions from an analytical and numerical point of view. From our results, it is clear that the diffusive regime is directly related to the system metric, meaning that the distributions characterize anomalous diffusion, which may represent usual, sub-, or super-diffusive processes. It is also interesting to note that each system has an influence on the spread of the other through the diffusive term and the external forces applied to the systems. In this manner, the results have shown that the dynamic of each one has a direct influence on the other as a thermal bath, since the coupling appears in the diffusive term, which is related to how the system will spread. Our results, as they have been obtained, represent a very general approach to describe particle dynamics and thermodynamic connection in systems composed of interacting particles, such as diffusion in crowded biological media [26] and many others. In particular, our results help to demonstrate how entropy production occurs in such systems which, in turn, may be fundamental to understanding several aspects of interacting particle systems and their connection with measurable quantities. Furthermore, the generalization proposed in this work allows one or both subsystems to be fractal in nature, that is, the space morphology is also present within the model. On the other hand, our results indicate that the space metric does not affect the entropy production of the system. Instead, it is directly affected by other aspects of the subsystems related to the intrinsic nature of diffusion. We anticipate that our findings will contribute to a better understanding of the connection between complex systems, nonlinear sciences, and the metric of space in which the process takes place.

Author Contributions

Conceptualization, R.S.Z., E.K.L., L.R.d.S. and M.K.L.; methodology, R.S.Z., E.K.L., L.R.d.S. and M.K.L.; validation, R.S.Z., E.K.L., L.R.d.S. and M.K.L.; formal analysis, R.S.Z., E.K.L., L.R.d.S. and M.K.L.; investigation, R.S.Z., E.K.L., L.R.d.S. and M.K.L.; writing—original draft preparation, R.S.Z., E.K.L., L.R.d.S. and M.K.L.; writing—review and editing, R.S.Z., E.K.L., L.R.d.S. and M.K.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in this article.

Acknowledgments

The authors thank the CNPq for the partial financial support.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Figure 1 shows the behavior of particle distribution obtained from [Forumla omitted. See PDF.] versus [Forumla omitted. See PDF.] obtained from Equations (41) and (42) for different values of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. The solid black line corresponds to the standard diffusion, and the other lines correspond to the generalized cases. To make things simpler, we assume [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 2. Figure 2 shows the behavior of the mean square displacement obtained for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] from Equations (3) and (26). The solid black line corresponds to the standard diffusion. The green dashed line corresponds to the superdiffusion. The red dashed-dotted line is associated with subdiffusion. To make things simpler, we assume [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 3. Behavior obtained for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] from Equations (3) and (6) in Figure 3a and Figure 3b, where [Forumla omitted. See PDF.] is given by Equation (14) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] ([Forumla omitted. See PDF.] is the Heviside function), and [Forumla omitted. See PDF.]. The black dashed lines correspond to [Forumla omitted. See PDF.] for different values of t. The red dashed-dotted lines correspond to [Forumla omitted. See PDF.] for different values of t. We consider, for simplicity, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 4. Figure 4a and Figure 4b show the behavior of the mean square displacement obtained with the distributions [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] from Equations (3) and (6), where [Forumla omitted. See PDF.] is given by Equation (14) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] ([Forumla omitted. See PDF.] is the Heaviside function) and [Forumla omitted. See PDF.]. The dashed-dotted lines correspond to the mean square displacement obtained from [Forumla omitted. See PDF.] for different values of [Forumla omitted. See PDF.]. dashed lines correspond to the mean square displacement obtained from [Forumla omitted. See PDF.] for different values of [Forumla omitted. See PDF.]. We consider, for simplicity, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 5. Figure 5 shows the behavior obtained for Equation (22) with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] from Equations (3) and (6), where [Forumla omitted. See PDF.] is given by Equation (14) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] ([Forumla omitted. See PDF.] is the Heviside function), and [Forumla omitted. See PDF.]. The red dashed-dotted line corresponds to the case [Forumla omitted. See PDF.]. The green dashed line corresponds to the case [Forumla omitted. See PDF.]. The inset corresponds to the behavior of Equation (57) for these values of [Forumla omitted. See PDF.]. We consider, for simplicity, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 6. Figure 6 shows the behavior obtained for Equation (22) with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] from Equations (3) and (6), where [Forumla omitted. See PDF.] is given by Equation (14) with [Forumla omitted. See PDF.]. The black dashed-dotted line corresponds to the case [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The green dashed line corresponds to the case [Forumla omitted. See PDF.] ([Forumla omitted. See PDF.] is the Heviside function), and [Forumla omitted. See PDF.]. The blue dashed-dotted-dotted line corresponds to the case [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. For simplicity, we consider [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Appendix A

Let us consider in detail the calculations of some equations. We start with the results obtained for Equation (13). For this, we consider the following Equation (A1), (A1)$$\begin{array}{c}\hfill \frac{d}{dt}F={\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\Xi (x_1,x_2)-kT\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right]\frac{\partial }{\partial t}\left[{\rho }_1(x_1,t){\rho }_2(x_2,t)\right],\end{array}$$ in which ${\rho }_{12}={\rho }_1{\rho }_2$, as defined before. After performing the time derivative, Equation (A1) can be written as follows (A2)$$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\Xi (x_1,x_2)-kT\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right]\hfill \\ & \times & \left[{\rho }_2(x_2,t)\frac{\partial }{\partial t}{\rho }_1(x_1,t)+{\rho }_1(x_1,t)\frac{\partial }{\partial t}{\rho }_2(x_2,t)\right].\hfill \end{array}$$

Now, we need to consider the nonlinear Fokker–Planck equations, i.e., Equations (3) and (6), to simplify the differential operators in time. By substituting Equations (3) and (6) in Equation (A2), we have that (A3)$$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\Xi (x_1,x_2){\rho }_2-kT{\rho }_2\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left[F_1\left(x_1\right){\rho }_1(x_1,t)\right]\right\}\hfill \\ & +& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\Xi (x_1,x_2){\rho }_1-kT{\rho }_1\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left[F_2\left(x_2\right){\rho }_2(x_2,t)\right]\right\}.\hfill \end{array}$$

Before performing integration by parts, we can use the properties of the fractal derivative and write the previous equation as follows:(A4)$$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& {\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\int }_{-\infty }^{\infty }d{x}_1\left\{\Xi (x_1,x_2){\rho }_2-kT{\rho }_2\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \frac{\partial }{\partial x_1}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left[F_1\left(x_1\right){\rho }_1(x_1,t)\right]\right\}\hfill \\ & +& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2\left\{\Xi (x_1,x_2){\rho }_1-kT{\rho }_1\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \frac{\partial }{\partial x_2}\left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left[F_2\left(x_2\right){\rho }_2(x_2,t)\right]\right\}.\hfill \end{array}$$

We perform some integration by parts and assume the following conditions: ${\rho }_i(x\to \pm \infty ,t)\to 0$ and ${\partial }_{{\xi }_{\mathcal{I},i}\left(x_i\right)}{\rho }_i(x_i,t){|}_{x_i\to \pm \infty }\to 0$. Thus, Equation (A4) yields:(A5)$$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& -{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right){\int }_{-\infty }^{\infty }d{x}_1\frac{\partial }{\partial x_1}\left\{\Xi (x_1,x_2){\rho }_2-kT{\rho }_2\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\left[F_1\left(x_1\right){\rho }_1(x_1,t)\right]\right\}\hfill \\ & -& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2\frac{\partial }{\partial x_2}\left\{\Xi (x_1,x_2){\rho }_1-kT{\rho }_1\frac{\partial }{\partial {\rho }_{12}}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)-\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\left[F_2\left(x_2\right){\rho }_2(x_2,t)\right]\right\}.\hfill \end{array}$$ and, consequently, (A6)$$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& -{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_2-kT{\rho }_2^2\frac{\partial {\rho }_1}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{P}_1({\rho }_1,t)+\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_1(x_1,t)\right\}\hfill \\ & -& {\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right){\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left\{\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1-kT{\rho }_1^2\frac{\partial {\rho }_2}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right\}\hfill \\ & \times & \left\{\mathcal{D}\frac{\partial }{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{P}_2({\rho }_2,t)+\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_2(x_2,t)\right\}.\hfill \end{array}$$

Note that to obtain Equation (A7), i.e., (A7)$$\begin{array}{ccc}\hfill \frac{d}{dt}F& =& -{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\frac{1}{{\rho }_1}\hfill \\ & \times & {\left\{{\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\left[\frac{\partial {\varphi }_1\left(x_1\right)}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}{\rho }_2{\rho }_1-kT{\rho }_2^2{\rho }_1\frac{\partial {\rho }_1}{\partial {\xi }_{\mathcal{I},1}\left(x_1\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right]\right\}}^2\hfill \\ & -& {\int }_{-\infty }^{\infty }d{x}_2{\xi }_2\left(x_2\right)\frac{1}{{\rho }_2}\hfill \\ & \times & {\left\{{\int }_{-\infty }^{\infty }d{x}_1{\xi }_1\left(x_1\right)\left[\frac{\partial {\varphi }_2\left(x_2\right)}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}{\rho }_1{\rho }_2-kT{\rho }_1^2{\rho }_2\frac{\partial {\rho }_2}{\partial {\xi }_{\mathcal{I},2}\left(x_2\right)}\frac{{\partial }^2}{\partial {\rho }_{12}^2}s\left({\rho }_{12}\right)\right]\right\}}^2.\hfill \end{array}$$ by taking Equations (30) and (31) into account, we have to consider (A8)$$\begin{array}{ccc}\hfill -kT{\rho }_j^2{\rho }_i\frac{{\partial }^2}{\partial {\rho }_{ij}^2}s\left({\rho }_{ij}\right)& =& \mathcal{D}{\rho }_j\frac{\partial }{\partial {\rho }_{ij}}\left({\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\rho }_i^{\overline{\gamma }}{\rho }_j^{\overline{\gamma }}\right),\hfill \end{array}$$ (A9)$$\begin{array}{ccc}& =& \mathcal{D}{\rho }_j\frac{\partial }{\partial {\rho }_{ij}}\left({\int }_0^{\gamma }d\overline{\gamma }p\left(\overline{\gamma }\right){\rho }_{ij}^{\overline{\gamma }}\right),\hfill \end{array}$$ where $i=1,2$ and $j=1,2$ with $i\ne j$, $\mathcal{D}=kT$, ${\rho }_{ij}={\rho }_i{\rho }_j$, and ${\int }_{-\infty }^{\infty }d{x}_i{\xi }_i\left(x_i\right){\rho }_i(x_i,t)=1$, which implies (A10)$$\begin{array}{c}\hfill \frac{d}{dt}F\le 0.\end{array}$$

Similar calculations were performed to obtain Equations (29)–(38).

References

1. Boltzmann, L. Further Studies on the Thermal Equilibrium of Gas Molecules. The Kinetic Theory of Gases; World Scientific: Singapore, 2003; pp. 262-349.

2. Gibbs, J.W. Elementary Principles in Statistical Mechanics: Developed with Especial Reference to the Rational Foundation of Thermodynamics; Cambridge Library Collection-Mathematics, Cambridge University Press: Cambridge, UK, 2010.

3. Lenzi, E.; Mendes, R.; Da Silva, L. Statistical mechanics based on Renyi entropy. Physica A; 2000; 280, pp. 337-345. [DOI: https://dx.doi.org/10.1016/S0378-4371(00)00007-8]

4. Rajagopal, A.; Abe, S. Implications of form invariance to the structure of nonextensive entropies. Phys. Rev. Lett.; 1999; 83, 1711. [DOI: https://dx.doi.org/10.1103/PhysRevLett.83.1711]

5. Kaniadakis, G.; Scarfone, A.; Sparavigna, A.; Wada, T. Composition law of κ-entropy for statistically independent systems. Phys. Rev. E; 2017; 95, 052112. [DOI: https://dx.doi.org/10.1103/PhysRevE.95.052112] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28618546]

6. Tsallis, C. Nonadditive entropy: The concept and its use. Eur. Phys. J. A; 2009; 40, 257. [DOI: https://dx.doi.org/10.1140/epja/i2009-10799-0]

7. Chavanis, P. Relaxation of a test particle in systems with long-range interactions: Diffusion coefficient and dynamical friction. Eur. Phys. J. B-Condens. Matter Complex Syst.; 2006; 52, pp. 61-82. [DOI: https://dx.doi.org/10.1140/epjb/e2006-00268-8]

8. Bologna, M.; Tsallis, C.; Grigolini, P. Anomalous diffusion associated with nonlinear fractional derivative Fokker-Planck-like equation: Exact time-dependent solutions. Phys. Rev. E; 2000; 62, pp. 2213-2218. [DOI: https://dx.doi.org/10.1103/PhysRevE.62.2213]

9. Lukyanov, A.V.; Sushchikh, M.; Baines, M.J.; Theofanous, T. Superfast nonlinear diffusion: Capillary transport in particulate porous media. Phys. Rev. Lett.; 2012; 109, 214501. [DOI: https://dx.doi.org/10.1103/PhysRevLett.109.214501] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23215595]

10. Aronson, D.G. The porous medium equation. Nonlinear Diffusion Problems: Lectures Given at the 2nd 1985 Session of the Centro Internazionale Matermatico Estivo (CIME) Held at Montecatini Terme, Italy 10–18 June 1985; Springer: Berlin/Heidelberg, Germany, 2006; pp. 1-46.

11. Tsallis, C.; Bukman, D.J. Anomalous diffusion in the presence of external forces: Exact time-dependent solutions and their thermostatistical basis. Phys. Rev. E; 1996; 54, R2197. [DOI: https://dx.doi.org/10.1103/PhysRevE.54.R2197] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9965440]

12. Drazer, G.; Wio, H.S.; Tsallis, C. Anomalous diffusion with absorption: Exact time-dependent solutions. Phys. Rev. E; 2000; 61, pp. 1417-1422. [DOI: https://dx.doi.org/10.1103/PhysRevE.61.1417]

13. Frank, T.D. Nonlinear Fokker-Planck Equations: Fundamentals and Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2005.

14. Marin, D.; Ribeiro, M.; Ribeiro, H.; Lenzi, E. A nonlinear Fokker–Planck equation approach for interacting systems: Anomalous diffusion and Tsallis statistics. Phys. Lett. A; 2018; 382, pp. 1903-1907. [DOI: https://dx.doi.org/10.1016/j.physleta.2018.05.009]

15. Liang, Y.; Chen, W.; Cai, W. Hausdorff Calculus: Applications to Fractal Systems; Fractional Calculus in Applied Sciences and Engineering, De Gruyter: Berlin, Germany, 2019.

16. Liang, Y.J.; Chen, W. Bridge fatigue life prediction using Mittag–Leffler distribution. Fatigue Fract. Eng. Mater. Struct.; 2014; 37, pp. 255-264. [DOI: https://dx.doi.org/10.1111/ffe.12110]

17. Liang, Y.; Chen, W. A survey on computing Lévy stable distributions and a new MATLAB toolbox. Signal Process.; 2013; 93, pp. 242-251. [DOI: https://dx.doi.org/10.1016/j.sigpro.2012.07.035]

18. Abdeljawad, T. On conformable fractional calculus. J. Comput. Appl. Math.; 2015; 279, pp. 57-66. [DOI: https://dx.doi.org/10.1016/j.cam.2014.10.016]

19. Chen, W. Time–space fabric underlying anomalous diffusion. Chaos Solitons Fractals; 2006; 28, pp. 923-929. [DOI: https://dx.doi.org/10.1016/j.chaos.2005.08.199]

20. Wei, Q.; Wang, W.; Zhou, H.; Metzler, R.; Chechkin, A. Time-fractional Caputo derivative versus other integrodifferential operators in generalized Fokker-Planck and generalized Langevin equations. Phys. Rev. E; 2023; 108, 024125. [DOI: https://dx.doi.org/10.1103/PhysRevE.108.024125]

21. Caputo, M. Linear Models of Dissipation whose Q is almost Frequency Independent—II. Geophys. J. Int.; 1967; 13, pp. 529-539. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1967.tb02303.x]

22. Satin, S.E.; Parvate, A.; Gangal, A. Fokker–Planck equation on fractal curves. Chaos Solitons Fractals; 2013; 52, pp. 30-35. [DOI: https://dx.doi.org/10.1016/j.chaos.2013.03.013]

23. Liang, Y.; Su, N.; Chen, W. A time-space Hausdorff derivative model for anomalous transport in porous media. Fract. Calc. Appl. Anal.; 2019; 22, pp. 1517-1536. [DOI: https://dx.doi.org/10.1515/fca-2019-0079]

24. Rosseto, M.; Ribeiro de Almeida, R.; Lenzi, E.; Evangelista, L.; Zola, R. Probing modulated liquid crystal media with dielectric spectroscopy. J. Mol. Liq.; 2023; 390, 122943. [DOI: https://dx.doi.org/10.1016/j.molliq.2023.122943]

25. Naderi, S.; Pouget, E.; Ballesta, P.; van der Schoot, P.; Lettinga, M.P.; Grelet, E. Fractional Hoppinglike Motion in Columnar Mesophases of Semiflexible Rodlike Particles. Phys. Rev. Lett.; 2013; 111, 037801. [DOI: https://dx.doi.org/10.1103/PhysRevLett.111.037801]

26. Banks, D.S.; Fradin, C. Anomalous Diffusion of Proteins Due to Molecular Crowding. Biophys. J.; 2005; 89, pp. 2960-2971. [DOI: https://dx.doi.org/10.1529/biophysj.104.051078]

27. Tarasov, V.E. Anisotropic fractal media by vector calculus in non-integer dimensional space. J. Math. Phys.; 2014; 55, 083510. [DOI: https://dx.doi.org/10.1063/1.4892155]

28. Vázquez, J.L. The Porous Medium Equation: Mathematical Theory; Oxford University Press: Oxford, UK, 2007.

29. Casas, G.A.; Nobre, F.D. Nonlinear Fokker-Planck equations in super-diffusive and sub-diffusive regimes. J. Math. Phys.; 2019; 60, 053301. [DOI: https://dx.doi.org/10.1063/1.5001189]

30. Plastino, A.R.; Wedemann, R.S.; Tsallis, C. Nonlinear fokker-planck equation for an overdamped system with drag depending on direction. Symmetry; 2021; 13, 1621. [DOI: https://dx.doi.org/10.3390/sym13091621]

31. Deppman, A.; Khalili Golmankhaneh, A.; Megías, E.; Pasechnik, R. From the Boltzmann equation with non-local correlations to a standard non-linear Fokker-Planck equation. Phys. Lett. B; 2023; 839, 137752. [DOI: https://dx.doi.org/10.1016/j.physletb.2023.137752]

32. Evangelista, L.R.; Lenzi, E.K. Nonlinear Fokker–Planck Equations, H-Theorem and Generalized Entropy of a Composed System. Entropy; 2023; 25, 1357. [DOI: https://dx.doi.org/10.3390/e25091357]

33. Schwämmle, V.; Nobre, F.D.; Curado, E.M.F. Consequences of the H theorem from nonlinear Fokker-Planck equations. Phys. Rev. E; 2007; 76, 041123. [DOI: https://dx.doi.org/10.1103/PhysRevE.76.041123]

34. dos Santos, M.; Lenzi, M.; Lenzi, E. Nonlinear Fokker–Planck equations, H–theorem, and entropies. Chin. J. Phys.; 2017; 55, pp. 1294-1299. [DOI: https://dx.doi.org/10.1016/j.cjph.2017.07.003]

35. Dos Santos, M.; Lenzi, E. Entropic nonadditivity, H theorem, and nonlinear Klein-Kramers equations. Phys. Rev. E; 2017; 96, 052109. [DOI: https://dx.doi.org/10.1103/PhysRevE.96.052109] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29347710]

36. Chavanis, P.H. Nonlinear mean field Fokker-Planck equations. Application to the chemotaxis of biological populations. Eur. Phys. J. B; 2008; 62, pp. 179-208. [DOI: https://dx.doi.org/10.1140/epjb/e2008-00142-9]

37. Schwämmle, V.; Curado, E.M.; Nobre, F.D. A general nonlinear Fokker-Planck equation and its associated entropy. Eur. Phys. J. B; 2007; 58, pp. 159-165. [DOI: https://dx.doi.org/10.1140/epjb/e2007-00217-1]

38. Martınez, S.; Pennini, F.; Plastino, A. Thermodynamics’ zeroth law in a nonextensive scenario. Phys. A Stat. Mech. Its Appl.; 2001; 295, pp. 416-424. [DOI: https://dx.doi.org/10.1016/S0378-4371(01)00121-2]

39. Abe, S. Temperature of nonextensive systems: Tsallis entropy as Clausius entropy. Phys. A Stat. Mech. Its Appl.; 2006; 368, pp. 430-434. [DOI: https://dx.doi.org/10.1016/j.physa.2006.04.001]

40. Biró, T.; Ván, P. Zeroth law compatibility of nonadditive thermodynamics. Phys. Rev. E; 2011; 83, 061147. [DOI: https://dx.doi.org/10.1103/PhysRevE.83.061147] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21797342]

41. Casas, G.; Nobre, F.; Curado, E. Entropy production and nonlinear Fokker-Planck equations. Phys. Rev. E; 2012; 86, 061136. [DOI: https://dx.doi.org/10.1103/PhysRevE.86.061136]

42. Plastino, A.; Wedemann, R.; Nobre, F. H-theorems for systems of coupled nonlinear Fokker-Planck equations. Europhys. Lett.; 2022; 139, 11002. [DOI: https://dx.doi.org/10.1209/0295-5075/ac7c30]

43. Penrose, O. Foundations of Statistical Mechanics: A Deductive Treatment; Courier Corporation: Singapore, 2005.

44. Lebowitz, J.L.; Bergmann, P.G. Irreversible Gibbsian ensembles. Ann. Phys.; 1957; 1, pp. 1-23. [DOI: https://dx.doi.org/10.1016/0003-4916(57)90002-7]

45. Tsallis, C.; Cirto, L.J. Black hole thermodynamical entropy. Eur. Phys. J. C; 2013; 73, pp. 1-7. [DOI: https://dx.doi.org/10.1140/epjc/s10052-013-2487-6]

46. Tsallis, C. An introduction to nonadditive entropies and a thermostatistical approach to inanimate and living matter. Contemp. Phys.; 2014; 55, pp. 179-197. [DOI: https://dx.doi.org/10.1080/00107514.2014.900977]

47. Plastino, A.; Curado, E.; Nobre, F.; Tsallis, C. From the nonlinear fokker-planck equation to the vlasov description and back: Confined interacting particles with drag. Phys. Rev. E; 2018; 97, 022120. [DOI: https://dx.doi.org/10.1103/PhysRevE.97.022120]

48. Lenzi, E.; Mendes, R.; Tsallis, C. Crossover in diffusion equation: Anomalous and normal behaviors. Phys. Rev. E; 2003; 67, 031104. [DOI: https://dx.doi.org/10.1103/PhysRevE.67.031104] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12689052]

49. Tateishi, A.; Lenzi, E.; da Silva, L.; Ribeiro, H.; Picoli, S., Jr.; Mendes, R. Different diffusive regimes, generalized Langevin and diffusion equations. Phys. Rev. E; 2012; 85, 011147. [DOI: https://dx.doi.org/10.1103/PhysRevE.85.011147] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22400552]

50. Li, W.; Wang, Q.A.; Nivanen, L.; Le Méhauté, A. On different q-systems in nonextensive thermostatistics. Eur. Phys. J. B-Condens. Matter Complex Syst.; 2005; 48, pp. 95-100. [DOI: https://dx.doi.org/10.1140/epjb/e2005-00368-y]

51. Wang, Q.A.; Nivanen, L.; Méhauté, A.L. A composition of different q nonextensive systems with the normalized expectation based on escort probability. arXiv; 2006; arXiv: cond-mat/0601255

52. Nivanen, L.; Pezeril, M.; Wang, Q.; Le Méhauté, A. Applying incomplete statistics to nonextensive systems with different q indices. Chaos Solitons Fractals; 2005; 24, pp. 1337-1342. [DOI: https://dx.doi.org/10.1016/j.chaos.2004.09.064]

53. Malacarne, L.; Mendes, R.; Pedron, I.; Lenzi, E. Nonlinear equation for anomalous diffusion: Unified power-law and stretched exponential exact solution. Phys. Rev. E; 2001; 63, 030101. [DOI: https://dx.doi.org/10.1103/PhysRevE.63.030101]

54. Prato, D.; Tsallis, C. Nonextensive foundation of Lévy distributions. Phys. Rev. E; 1999; 60, 2398. [DOI: https://dx.doi.org/10.1103/PhysRevE.60.2398]

55. Tsallis, C.; Levy, S.V.; Souza, A.M.; Maynard, R. Statistical-mechanical foundation of the ubiquity of Lévy distributions in nature. Phys. Rev. Lett.; 1995; 75, 3589. [DOI: https://dx.doi.org/10.1103/PhysRevLett.75.3589]

56. O’Shaughnessy, B.; Procaccia, I. Analytical solutions for diffusion on fractal objects. Phys. Rev. Lett.; 1985; 54, 455. [DOI: https://dx.doi.org/10.1103/PhysRevLett.54.455]

57. O’Shaughnessy, B.; Procaccia, I. Diffusion on fractals. Phys. Rev. A; 1985; 32, 3073. [DOI: https://dx.doi.org/10.1103/PhysRevA.32.3073]

58. Zhao, D.; Yu, G.; Li, W. Diffusion on fractal objects modeling and its physics-informed neural network solution. Fractals; 2021; 29, 2150071. [DOI: https://dx.doi.org/10.1142/S0218348X21500717]

59. Pascal, J.; Pascal, H. On nonlinear diffusion in fractal structures. Phys. A Stat. Mech. Its Appl.; 1994; 208, pp. 351-358. [DOI: https://dx.doi.org/10.1016/0378-4371(94)00052-2]

60. Crank, J. The Mathematics of Diffusion; Oxford University Press: Oxford, UK, 1979.

61. Burden, R.L.; Faires, J.D. Numerical Analysis; 9th ed. PWS Publishing Co.: Boston, MA, USA, 2011.