Anomalous Diffusion and Non-Markovian Reaction of Particles near an Adsorbing Colloidal Particle

Abstract

We investigate the diffusion phenomenon of particles in the vicinity of a spherical colloidal particle where particles may be adsorbed/desorbed and react on the surface of the colloidal particle. The mathematical model comprises a generalized diffusion equation to govern bulk dynamics and kinetic equations which can describe non-Debye relaxations and is used for the colloid’s surface. For the reaction processes, we also consider the presence of convolution kernels, which offer the flexibility of describing a single process or process with intermediate reactions before forming the final species. Our analysis focuses on analytical and numerical calculations to obtain the particles’ behavior on the colloidal particle’s surface and to determine how it affects the diffusion of particles around it. The solutions obtained show various behaviors that can be connected to anomalous diffusion phenomena and may be used to describe the ever-richer science of colloidal particles better.

1. Introduction

Colloidal dispersions are of fundamental importance for human beings. They are integral to many everyday materials, including paints, inks, food products, and biological fluids. Moreover, these heterogeneous systems straddle the boundary between classical and quantum science. They exhibit a complex interplay of phenomena that, once understood, can influence the production of state-of-the-art materials [1]. One of the most ubiquitous phenomena involving colloidal particles (often spherical) is the adsorption of particles [2] present in solution (free ions, for example) [1], which dictates how colloids interact and how stable the system is [1,3,4]. Thus, it is essential to understand how the presence of a colloidal particle affects the diffusing species around it and how adsorption and reaction occur near the spherical particle.

In fact, beyond colloidal sciences, diffusion is a fundamental process that permeates nature from physics to biology [5,6]. For instance, in active transport [7,8], molecular crowding [9,10,11], cellular membranes [12], kinetics on surfaces [13,14,15], subdiffusion in thin membranes [16,17,18], electrical response [19,20], and diffusion on fractals [21,22]. Near a colloidal particle, these scenarios are complex, and combine processes that occur in bulk and on the surface. As a result, we may have a normal or anomalous diffusion present in these systems, characterized by a linear dependence on the mean square displacement, i.e., $⟨{(r-⟨r⟩)}^2⟩\sim t$, typical of a Markovian process. The second case, i.e., anomalous diffusion, has a nonlinear time dependence for the mean square displacement, e.g., $⟨{(r-⟨r⟩)}^2⟩\sim t^{\alpha }$ ($\alpha <1$ and $\alpha >1$ correspond to sub- and superdiffusion [20,23], respectively) typical of non-Markovian processes.

Different approaches from the experimental and theoretical points of view have been successfully used to analyze and interpret the complexity of these contexts. Some representative examples are diffusion in cells [24,25], dialysis [26,27], microscopy techniques [28,29], and electrochemical methods [30,31] on the experimental side. From the theoretical point of view, generalized Langevin equations [32], Fokker–Planck equations [33], master equations [34,35], and nonlinear [36,37] and/or fractional Fokker–Planck equations [20,38] have been used to tackle the relevant problems raised in these scenarios. The advances in each field, experimental or theoretical, are challenging and bring the possibility of exploring different aspects of these systems.

Here, we investigate the dynamics of two species of neutral particles diffusing in the vicinity of a colloidal particle of radius R, which is a system governed by a generalized diffusion equation, where particles can undergo adsorption–desorption or promote a reaction process on the surface of a spherical colloidal particle with the formation of different species (see Figure 1). We model the size of the colloidal particle much larger than the diffusing particles; i.e., particles 1 and 2 have dimensions much smaller than the colloidal particles with a spherical surface where the adsorption–desorption process will occur. In particular, colloidal particles usually range from 1 nanometer to 1 $\mathsf{\mu }$m in diameter. They are small enough to remain suspended in a medium, such as a liquid or gas, and can exhibit properties like Brownian motion. These dimensions allow colloids to form systems, including sols, gels, and emulsions. Considering the adsorption–desorption process governed by kinetic equations, stochastic dynamics can explain a wide variety of phenomena [39,40]. For the reaction processes, we also consider the presence of convolution kernels, which offer the flexibility of describing a single process or process with intermediate reactions before forming the final species. Our analysis employs analytical and numerical methods to obtain the behavior of the particles on the surface and in the bulk. These developments are performed in Section 2 and Section 3, wherein we consider two species of particles considering the dynamic coupling provided by the boundary conditions. The discussion of the results and some concluding remarks are presented in Section 4 and Section 5.

2. Diffusion, Reaction, and Surface Dynamics

Let us start our analysis by assuming that the particles of species 1 and 2 around the colloidal particle are governed by the following diffusion equation:

(1)$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial t}}{\rho }_{1\left(2\right)}(r,t)={\mathcal{D}}_{1\left(2\right)}{\mathfrak{F}}_t\left\{{\nabla }^2{\rho }_{1\left(2\right)}(r,t)\right\},\end{array}$

where ${\rho }_{1\left(2\right)}(r,t)$ represents the density of particles related to each species. In Equation (1), the fractional operator ${\mathfrak{F}}_t\{\cdots \}$ is defined as follows:

(2)$\begin{array}{c}\hfill {\mathfrak{F}}_t\left\{{\rho }_{1\left(2\right)}(x,t)\right\}={\displaystyle \frac{\partial }{\partial t}}{\int }_0^{t}d{t}^{\prime }{\mathfrak{K}}_{1\left(2\right)}(t-t^{\prime })\rho (x,t^{\prime }),\end{array}$

where the kernel ${\mathfrak{K}}_{1\left(2\right)}\left(t\right)$ is connected with the memory effects present in the bulk. Before proceeding, we remark that the statement of the problem is such that according to the kernel, different time-fractional derivatives can be obtained. Some examples are

(3)${\mathfrak{K}}_{1\left(2\right)}\left(t\right)={\displaystyle \frac{t^{-\gamma }}{\Gamma (1-\gamma )}},0<\gamma <1,$

which is connected with the Riemann–Liouville differential operator [20];

${\mathfrak{K}}_{1\left(2\right)}\left(t\right)={\displaystyle \frac{\mathcal{N}\left(\gamma \right)}{1-\gamma }}{e}^{-\gamma t/(1-\gamma )},$

where $\mathcal{N}\left(\gamma \right)$ is the normalization factor, which is connected to the Caputo–Fabrizio differential operator [41]; and

${\mathfrak{K}}_{1\left(2\right)}\left(t\right)={\displaystyle \frac{{\mathcal{N}}^{\prime }\left(\gamma \right)}{1-\gamma }}{E}_{\gamma }\left[-\gamma t^{\gamma }/(1-\gamma )\right],$

where $E_{\gamma }\left(x\right)$ is the Mittag–Leffler function and ${\mathcal{N}}^{\prime }\left(\gamma \right)$ is the normalization factor, which is the Antangana–Baleanu differential operator [42,43]. It is also possible to consider different kernels, and consequently, to obtain other differential operators [44,45]. Additional contexts where general fractional differential operators have been used can be found in Refs. [46,47,48,49]. Equation (1) may be obtained from different approaches, one of them is the continuous random walk [20,38] by considering suitable waiting time $\omega \left(t\right)$ and jumping $\lambda \left(r\right)$ distributions. In this approach, we have a probability density function $\psi (r;t)$, from which it is possible to obtain the waiting time distribution and the jumping distributions as follows:

(4)$\begin{array}{c}\hfill \omega \left(t\right)={\int }_{\mathcal{V}}dr{r}^2\psi (r;t)\mathrm{a}\mathrm{n}\mathrm{d}\lambda \left(r\right)={\int }_0^{\infty }dt\psi (r;t),\end{array}$

where $\mathcal{V}$ is the region accessible to the system to diffuse into. By using $\psi (r,t)$, the distribution related to the diffusion process $\rho (r,t)$ is found by combining the equations

(5)$\begin{array}{c}\hfill \eta (r,t)={\int }_{\mathcal{V}}d{r}^{\prime }{r\prime }^2{\int }_0^{t}d{t}^{\prime }\eta \left(r^{\prime },t^{\prime }\right)\psi \left(r,r^{\prime };t-t^{\prime }\right)d{t}^{\prime }+\delta \left(t\right)\phi \left(r\right)\end{array}$

and

(6)$\begin{array}{c}\hfill \rho (r,t)={\int }_0^t\eta \left(r,t-t^{\prime }\right)\Phi \left(t^{\prime }\right)d{t}^{\prime },\end{array}$

with $\Phi \left(t\right)=1-{\int }_0^t\omega \left(\tau \right)d\tau$. In Equation (5), $\eta (r,t)$ can be connected to the probability density of just arriving at r at time t with the event of having just arrived at $r^{\prime }$ at time $t^{\prime }$, $\eta (r^{\prime },t^{\prime })$. The second term in Equation (5) is the initial condition, for simplicity, chosen to be $\phi \left(r^{\prime }\right)$. By using the Laplace transform ($\overline{\rho }(x,s)$=$\mathcal{L}\left\{\rho \right(x,t);s\}$) and an integral transform connected to the symmetry properties of the $\psi (r,t)$ with relation to the spatial variable, denoted by $\tilde{\rho }(k,t)={\mathfrak{F}}_{\psi }\{\rho (r,t);k\}$, it is possible to obtain

(7)$\begin{array}{c}\hfill \tilde{\overline{\rho }}(k,s)={\displaystyle \frac{1-\overline{\omega }\left(s\right)}{s}}{\displaystyle \frac{\tilde{\phi }\left(k\right)}{1-\tilde{\overline{\psi }}(k,s)}}.\end{array}$

$\omega \left(t\right)$ is given by [50]:

(8)$\begin{array}{c}\hfill \omega \left(t\right)={\displaystyle \frac{1}{1+1/\mathfrak{K}\left(t\right)}}\end{array}$

and asymptotically $\widehat{\lambda }\left(k\right)\sim 1-\mathcal{D}{\tilde{\lambda }}_{\psi }\left(k\right)$, with a short-tailed behavior. By assuming that $\tilde{\overline{\psi }}(k,s)=\tilde{\lambda }\left(k\right)\overline{\omega }\left(s\right)$, we obtain

(9)$\begin{array}{c}\hfill \tilde{\overline{\rho }}(k,s)={\displaystyle \frac{1/\left[s\mathfrak{K}\right(s\left)\right]}{1/\mathfrak{K}\left(s\right)+{\widehat{\lambda }}_{\psi }\left(k\right)}}.\end{array}$

By performing the inverse integral transforms, we obtain that

(10)$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial t}}\rho (r,t)=\mathcal{D}{\displaystyle \frac{\partial }{\partial t}}{\int }_0^{t}d{t}^{\prime }\mathfrak{K}(t-t^{\prime }){\nabla }_{\psi }^2\rho (r^{\prime },t),\end{array}$

with ${\nabla }_{\psi }^2(\cdots )\equiv {\mathfrak{F}}_{\psi }^{-1}\{{\tilde{\lambda }}_{\psi }\left(k\right);r\}$, which is essentially Equation (1). In particular, the symmetry of our problem is suitable for the Bessel functions with spherical symmetries, which implies that the integral transform can be connected to the Weber integral transform [51] and ${\lambda }_{\psi }\left(k\right)=-k^2$; we have that ${\nabla }_{\psi }^2(\cdots )\equiv r^{-2}{\partial }_r\left(r^2{\partial }_r(\cdots )\right)$.

Equation (1) is also subjected to the boundary conditions

(11)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{D}}_1{\left.\oint d\mathcal{A}·{\mathfrak{F}}_t\left\{\nabla {\rho }_1(r,t)\right\}\right|}_{r=\mathcal{R}}& ={\displaystyle \frac{d}{dt}}{\rho }_{\mathrm{surf},1}\left(t\right)+{\int }_0^t{k}_{s,11}(t-t^{\prime }){\rho }_1(\mathcal{R},t^{\prime })d{t}^{\prime }\hfill \\ & -{\int }_0^t{k}_{s,12}(t-t^{\prime }){\rho }_2(\mathcal{R},t^{\prime })d{t}^{\prime }\hfill \end{array}\end{array}$

and

(12)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{D}}_2{\left.\oint d\mathcal{A}·{\mathfrak{F}}_t\left\{\nabla {\rho }_2(r,t^{\prime })\right\}\right|}_{r=\mathcal{R}}& ={\displaystyle \frac{d}{dt}}{\rho }_{\mathrm{surf},2}\left(t\right)+{\int }_0^t{k}_{s,22}(t-t^{\prime }){\rho }_2(\mathcal{R},t^{\prime })d{t}^{\prime }\hfill \\ & -{\int }_0^t{k}_{s,21}(t-t^{\prime }){\rho }_1(\mathcal{R},t^{\prime })d{t}^{\prime },\hfill \end{array}\end{array}$

which connect the flux through the surface of the spherical colloidal particle (in $r=R)$ with the changes due to the processes undergone by particles in the vicinity. This set of equations implies that species 1 and 2 can be sprayed on the surface and, through a reaction process, promotes the production of species 2 and 1. Furthermore, these equations can be used as a toy model to analyze a biological cell that, by a membrane sorption process, detects the presence of different species, for example, species 1, and produces another species, for example, species 2, which is released to the bulk and can react with species 1 in the bulk, detected by the sorption process, or stimulates another process or the production of other species. In addition to these boundary conditions, we have to take into account that

(13)$\begin{array}{c}\hfill {\mathcal{D}}_{1\left(2\right)}{\left.\oint d\mathcal{A}·{\mathfrak{F}}_t\left\{\nabla {\rho }_{1\left(2\right)}(r,t)\right\}\right|}_{r=\infty }=0.\end{array}$

Note that in Equations (11) and (12), we have terms related to the adsorption–desorption processes, ${\rho }_{\mathrm{surf},1\left(2\right)}\left(t\right)$, and terms related to reaction processes on the surface with consumption or production of species. These reaction processes are governed by the kernels $k_{s,11\left(22\right)}\left(t\right)$ and $k_{s,12\left(21\right)}\left(t\right)$, which define if the reaction is irreversible or reversible. For the adsorption–desorption processes on the surface, we consider

(14)$\begin{array}{c}\hfill {\rho }_{\mathrm{surf},1\left(2\right)}\left(t\right)={\rho }_{\mathrm{surf},1\left(2\right)}^{\left(0\right)}\left(t\right)+{\int }_0^{t}d{t}^{\prime }{\kappa }_{1\left(2\right)}(t-t^{\prime }){\rho }_{1\left(2\right)}(\mathcal{R},t^{\prime }),\end{array}$

where ${\rho }_{\mathrm{surf},1\left(2\right)}^{\left(0\right)}\left(t\right)$ can be related to an arbitrary initial condition, and the kernel ${\kappa }_{1\left(2\right)}\left(t\right)$ is related to the kinetic aspects occurring on the surface. Typical expressions for the kernel are ${\kappa }_{1\left(2\right)}\left(t\right)=constant$ and ${\kappa }_{1\left(2\right)}\left(t\right)=\left(k/{\tau }_{1\left(2\right)}\right)exp\left(-t/{\tau }_{1\left(2\right)}\right)$. The first expression represents a reaction process of the particles with the surface when the surface absorbs the particles. The second can be related to an adsorption–desorption process with a characteristic time ${\tau }_{1\left(2\right)}$. In particular, this case can be directly related to a kinetic process of first order, i.e.,

(15)$\begin{array}{c}\hfill {\displaystyle \frac{d}{dt}}{\rho }_{\mathrm{surf},1\left(2\right)}\left(t\right)=k_{1\left(2\right)}{\rho }_{1\left(2\right)}(\mathcal{R},t)-{\displaystyle \frac{1}{{\tau }_{1\left(2\right)}}}{\rho }_{\mathrm{surf},1\left(2\right)}\left(t\right),\end{array}$

where $k_{1\left(2\right)}{\tau }_{1\left(2\right)}$ may be related to the distance from the colloidal particle that defines the region of the interaction between the surface and the bulk, where the particles in front of the surface can be adsorbed. Again, in view of the generality of the approach, other expressions for $\kappa \left(t\right)$ may also be taken into account, giving rise to different kinetic processes, such as the ones analyzed in Refs. [52,53,54,55]. From the previous equations, we have

(16)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}\left\{{\rho }_{\mathrm{surf},1\left(2\right)}\left(t\right)+4\pi {\int }_{\mathcal{R}}^{\infty }dr{r}^2{\rho }_{1\left(2\right)}(r,t)\right\}& =-{\int }_0^t{k}_{s,11\left(22\right)}(t-t^{\prime }){\rho }_{1\left(2\right)}(\mathcal{R},t^{\prime })d{t}^{\prime }\hfill \\ & +{\int }_0^t{k}_{s,12\left(21\right)}(t-t^{\prime }){\rho }_{2\left(1\right)}(\mathcal{R},t^{\prime })d{t}^{\prime }.\hfill \end{array}\end{array}$

The first term on the left-hand side is connected to the particles on the surface, and the other term is connected to the particles in the bulk. The terms on the right-hand side are related to the reaction process on the surface, where one species promotes the reaction process by forming the other species.

3. Time-Dependent Solutions

Let us focus on the time-dependent solutions of Equation (1) when the previous boundary conditions are invoked. We start our analysis by considering that the initial conditions of the system are given by ${\rho }_{1\left(2\right)}(r,0)={\phi }_{1\left(2\right)}\left(r\right)$ and ${\rho }_{\mathrm{surf},1\left(2\right)}\left(0\right)$= 0, which account for a case where there are no adsorbed particles on the colloid at $t=0$, so all particles are initially in the bulk. To solve the diffusion equation related to each species, which are coupled by the boundary conditions, we may use the Green’s function approach and the Laplace transform, i.e.,

(17)$\begin{array}{c}\hfill \begin{array}{cc}& \mathcal{L}\{\rho (r,t);s\}={\int }_0^{\infty }dt{e}^{-st}\rho (r,t)=\widehat{\rho }(r,s)\mathrm{and}\hfill \\ & {\mathcal{L}}^{-1}\{\widehat{\rho }(r,s);t\}={\displaystyle \frac{1}{2\pi i}}{\int }_{c-i\infty }^{c+i\infty }dt{e}^{-st}\widehat{\rho }(r,s)=\rho (r,t).\hfill \end{array}\end{array}$

By using the Green’s function approach, the solutions for these equations can be found, and they are given by

(18)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_1(r,t)& =4\pi {\int }_{\mathcal{R}}^{\infty }{\mathcal{G}}_1(r,r^{\prime };t){\phi }_1\left(r^{\prime }\right){r\prime }^{2}d{r}^{\prime }\hfill \\ & +{\int }_0^{t}d{t}^{\prime }{\mathcal{G}}_1(r,\mathcal{R};t-t^{\prime }){\int }_0^{t^{\prime }}d{t}^{″}{k}_{12}(t^{\prime }-t^{″}){\rho }_2(\mathcal{R},t^{″})\hfill \end{array}\end{array}$

and

(19)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_2(r,t)& =4\pi {\int }_{\mathcal{R}}^{\infty }{\mathcal{G}}_2(r,r^{\prime };t){\phi }_1\left(r^{\prime }\right){r\prime }^{2}d{r}^{\prime }\hfill \\ & +{\int }_0^{t}d{t}^{\prime }{\mathcal{G}}_2(r,\mathcal{R};t-t^{\prime }){\int }_0^{t^{\prime }}d{t}^{″}{k}_{21}(t^{\prime }-t^{″}){\rho }_1(\mathcal{R},t^{″}).\hfill \end{array}\end{array}$

Note that in Equations (18) and (19) the first term promotes the spread of the initial condition and the other terms represent the influence on the diffusive process of the reaction processes occurring on the surface. The Green’s function is obtained by solving the equation

(20)$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial t}}{\mathcal{G}}_{1\left(2\right)}(r,r^{\prime },t)-{\mathcal{D}}_{1\left(2\right)}{\mathfrak{F}}_t\left\{{\nabla }^2{\mathcal{G}}_{1\left(2\right)}(r,r^{\prime },t)\right\}={\displaystyle \frac{1}{4\pi r^2}}\delta (r-r^{\prime })\delta \left(t\right),\end{array}$

with the conditions ${\mathcal{G}}_{1\left(2\right)}(r,r^{\prime };t)=0$ for $t\le 0$,

(21)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{D}}_{1\left(2\right)}\oint d\mathcal{A}·{\left.{\mathfrak{F}}_t\left\{{\nabla }_{1\left(2\right)}{\mathcal{G}}_{1\left(2\right)}(r,r^{\prime },t)\right\}\right|}_{r=\mathcal{R}}& ={\displaystyle \frac{\partial }{\partial t}}{\int }_0^{t}d{t}^{\prime }{\kappa }_{1\left(2\right)}(t-t^{\prime }){\mathcal{G}}_{1\left(2\right)}(\mathcal{R},r^{\prime },t^{\prime })\hfill \\ & +{\int }_0^t{k}_{s,11\left(22\right)}(t-t^{\prime }){\mathcal{G}}_{1\left(2\right)}(R,r^{\prime },t^{\prime })d{t}^{\prime },\hfill \end{array}\end{array}$

and

(22)$\begin{array}{ccc}\hfill {\left.{\mathcal{D}}_{1\left(2\right)}\oint d\mathcal{A}·{\mathfrak{F}}_t\left\{{\nabla }_{1\left(2\right)}{\mathcal{G}}_{1\left(2\right)}(r,r^{\prime },t)\right\}\right|}_{r=\infty }& =& 0.\hfill \end{array}$

The boundary condition given by Equation (21) contains an integrodifferential operator, which depends on the expression of ${\kappa }_{1\left(2\right)}\left(t\right)$, related to the sorption–desorption processes on the surface. For a kinetic process governed by Equation (14), the kernel is given by an exponential function and may be related to the operators analyzed in Ref. [43]. Another expression, such as, in particular, accounting for a power-law behavior, may be related to the Riemann–Liouville fractional differential operator.

The solution for Equation (20) in the Laplace domain is

(23)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\mathcal{G}}}_{1\left(2\right)}(r,r^{\prime };s)& ={\displaystyle \frac{1}{8\pi r{r}^{\prime }\sqrt{s{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}}}\left(e^{-\sqrt{{\displaystyle \frac{s}{{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}}}|r-r^{\prime }|}-e^{-\sqrt{{\displaystyle \frac{s}{{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}}}|r+r^{\prime }-2\mathcal{R}|}\right)\hfill \\ & +{\displaystyle \frac{1}{r{r}^{\prime }}}{\displaystyle \frac{{\mathcal{R}}^2{e}^{-\sqrt{\frac{s}{{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}}|r+r^{\prime }-2\mathcal{R}|}}{4\pi {\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)\mathcal{R}\left(\mathcal{R}\sqrt{s/{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}+1\right)+{\widehat{\eta }}_{1\left(2\right)}\left(s\right)}},\hfill \end{array}\end{array}$

where ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)={\mathcal{D}}_{1\left(2\right)}s{\widehat{\mathfrak{K}}}_{1\left(2\right)}\left(s\right)$ and ${\widehat{\eta }}_{1\left(2\right)}\left(s\right)=s{\widehat{\kappa }}_{1\left(2\right)}\left(s\right)+{\widehat{k}}_{s,11\left(22\right)}\left(s\right)$, and the last term represents the influence of the adsorption–desorption processes on the spreading of the system.

To perform the inverse Laplace transform of Equation (23), we use some expressions for ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)$ and ${\widehat{\eta }}_{1\left(2\right)}\left(s\right)$ to illustrate the behavior of the Green’s function, as shown by Figure 2, Figure 3 and Figure 4. The cases that we consider are:

1. ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)={\mathcal{D}}_{1\left(2\right)}=\mathrm{constant}$ for ${\widehat{\eta }}_{1\left(2\right)}\left(s\right)$, and ${\widehat{\eta }}_{1\left(2\right)}\left(s\right)={\eta }_{1\left(2\right)}=\mathrm{constant}$ (Figure 2);

2. ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)={\mathcal{D}}_{1\left(2\right)}{s}^{1-\gamma }$, with ${\widehat{\eta }}_{1\left(2\right)}\left(s\right)$ (Figure 3);

3. ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)=s{\mathcal{D}}_{1\left(2\right)}/(s+1/\tau )$, with ${\widehat{\eta }}_{1\left(2\right)}\left(s\right)$ (Figure 4).

The previous choices for ${\mathcal{D}}_{1\left(2\right)}\left(s\right)$ are related to different forms of the differential operators present in Equation (1) used to perform the inverse Laplace transform. The first case, i.e., ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)={\mathcal{D}}_{1\left(2\right)}=constant$, corresponds to use of the kernel ${\widehat{\mathfrak{K}}}_{1\left(2\right)}\left(s\right)=1/s$, which leads us to the standard form of the diffusion equation. The second case, where ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)={\mathcal{D}}_{1\left(2\right)}{s}^{1-\gamma }$, implies the kernel ${\widehat{\mathfrak{K}}}_{1\left(2\right)}\left(s\right)=t^{\gamma -1}/\Gamma \left(\gamma \right)$, which corresponds a Riemann–Liouville fractional operator. This choice for the kernel leads us to obtain the fractional diffusion equations discussed in Ref. [38]. The last case is given by ${\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)=s{\mathcal{D}}_{1\left(2\right)}/(s+1/\tau )$. It corresponds to considering the kernel ${\widehat{\mathfrak{K}}}_{1\left(2\right)}\left(s\right)=e^{-t/\tau }$, which implies that the fractional derivative considered is the Caputo–Fabrizio one in this case. Other choices for the kernel may be connected to different fractional operators with singular or nonsingular kernels.

Notice that Figure 2 and Figure 3 show how the propagator ${\mathcal{G}}_{1\left(2\right)}(r,r^{\prime };s)$ behaves from the surface of the colloid for different forms of ${\mathcal{D}}_{1\left(2\right)}\left(s\right)$ when $t=1.0$, and different values of $\gamma$ and $\alpha$; while Figure 4 shows how ${\mathcal{G}}_{1\left(2\right)}(r,r^{\prime };s)$ behaves near the colloid as time evolves. The Green’s function, after performing the inverse Laplace transform, can be written for the first case as

(24)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{G}}_{1\left(2\right)}(r,r^{\prime };t)& ={\displaystyle \frac{1}{8\pi r{r}^{\prime }\sqrt{\pi {\mathcal{D}}_{1\left(2\right)}t}}}\left\{e^{-{\displaystyle \frac{{(r-r^{\prime })}^2}{4{\mathcal{D}}_{1\left(2\right)}t}}}-e^{-{\displaystyle \frac{{(r+r^{\prime }2\mathcal{R})}^2}{4{\mathcal{D}}_{1\left(2\right)}t}}}\right\}\hfill \\ & +{\displaystyle \frac{1}{r{r}^{\prime }}}{\mathcal{G}}_{1\left(2\right)}^{\left(1\right)}(r+r^{\prime }-2\mathcal{R};t)+{\displaystyle \frac{1}{r{r}^{\prime }}}\sum _{n=1}^{\infty }{(-1)}^n{\int }_0^{t}d{t}_n{\mathcal{I}}_{1\left(2\right)}^{\left(1\right)}(t-t_n)\hfill \\ & \times {\int }_0^{t_n}d{t}_{n-1}{\mathcal{I}}_{1\left(2\right)}^{\left(1\right)}(t_n-t_{n-1})\cdots {\int }_0^{t_2}d{t}_1{\mathcal{I}}_{1\left(2\right)}^{\left(1\right)}(t_2-t_1){\mathcal{G}}_{1\left(2\right)}^{\left(1\right)}(r+r^{\prime }-2\mathcal{R};t_1),\hfill \end{array}\end{array}$

with

(25)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{G}}_{1\left(2\right)}^{\left(1\right)}(r,t)& =r{\int }_0^{t}d{t}^{\prime }{\mathcal{I}}_{1\left(2\right)}^{\left(0\right)}(t-t^{\prime }){\displaystyle \frac{e^{-\frac{r^2}{4{\mathcal{D}}_{1\left(2\right)}{t}^{\prime }}}}{\sqrt{4\pi {\mathcal{D}}_{1\left(2\right)}{t}^{\prime 3}}}},\hfill \\ \hfill {\mathcal{I}}_{1\left(2\right)}^{\left(0\right)}\left(t\right)& ={\displaystyle \frac{1}{4\pi \sqrt{{\mathcal{D}}_{1\left(2\right)}}}}\left\{{\displaystyle \frac{1}{\sqrt{\pi }t}}-\left[e^{{\mathcal{D}}_{1\left(2\right)}t/{\mathcal{R}}^2}\sqrt{{\displaystyle \frac{{\mathcal{D}}_{1\left(2\right)}}{{\mathcal{R}}^2}}}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\sqrt{{\displaystyle \frac{{\mathcal{D}}_{1\left(2\right)}t}{{\mathcal{R}}^2}}}\right)\right]\right\},\hfill \end{array}\end{array}$

and ${\mathcal{I}}_{1\left(2\right)}^{\left(1\right)}=\left(1/{\mathcal{R}}^2\right){\int }_0^{t}d{t}^{\prime }{\eta }_{1\left(2\right)}\left(t^{\prime }\right){\mathcal{I}}_{1\left(2\right)}^{\left(0\right)}(t-t^{\prime })$, where $\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(x\right)$ is the complementary error function. For the particular case ${\widehat{\eta }}_{1\left(2\right)}\left(s\right)=constant$, we have that

(26)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{G}}_{1\left(2\right)}(r,r^{\prime };t)& ={\displaystyle \frac{1}{8\pi r{r}^{\prime }\sqrt{\pi {\mathcal{D}}_{1\left(2\right)}t}}}\left\{e^{-{\displaystyle \frac{{(r-r^{\prime })}^2}{4{\mathcal{D}}_{1\left(2\right)}t}}}-e^{-{\displaystyle \frac{{(r+r^{\prime }2\mathcal{R})}^2}{4{\mathcal{D}}_{1\left(2\right)}t}}}\right\}\hfill \\ & +{\displaystyle \frac{1}{r{r}^{\prime }}}\left(r+r^{\prime }-2\mathcal{R}\right){\int }_0^{t}d{t}^{\prime }\Xi \left(t^{\prime }\right){\displaystyle \frac{e^{-\frac{{(r+r^{\prime }-2\mathcal{R})}^2}{4{\mathcal{D}}_{1\left(2\right)}t}}}{\sqrt{4\pi }{\mathcal{D}}_{1\left(2\right)}{\left(t-t^{\prime }\right)}^{3/2}}},\hfill \end{array}\end{array}$

where

(27)$\begin{array}{c}\hfill \Xi \left(t\right)=\left[1/\sqrt{\pi t^{\prime }}-\overline{\eta }{e}^{{\overline{\eta }}^{2}t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\overline{\eta }\sqrt{t^{\prime }}\right)\right]/\left(4\pi \sqrt{{\mathcal{D}}_{1\left(2\right)}}\right)\end{array}$

with $\overline{\eta }=\left(\sqrt{{\mathcal{D}}_{1\left(2\right)}}/\mathcal{R}\right)$$[1$ + $\eta$/$\left(4\pi {\mathcal{D}}_{1\left(2\right)}\mathcal{R}\right)]$.

In the second case, the Green’s function is

(28)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{G}}_{1\left(2\right)}(r,r^{\prime };t)& ={\displaystyle \frac{1}{8\pi r{r}^{\prime }}}\left\{{\mathfrak{G}}_{1\left(2\right)}^{\left(0\right)}(r-r^{\prime };t)-{\mathfrak{G}}_{1\left(2\right)}^{\left(0\right)}(r+r^{\prime }-2\mathcal{R};t)\right\}+{\displaystyle \frac{1}{r{r}^{\prime }}}{\mathfrak{G}}_{1\left(2\right)}^{\left(1\right)}(r+r^{\prime }-2\mathcal{R};t)\hfill \\ & +{\displaystyle \frac{1}{r{r}^{\prime }}}\sum _{n=1}^{\infty }{(-1)}^n{\int }_0^{t}d{t}_n{{\rm Y}}_{1\left(2\right)}^{\left(0\right)}(t-t_n)\hfill \\ & \times {\int }_0^{t_n}d{t}_{n-1}{{\rm Y}}_{1\left(2\right)}^{\left(0\right)}(t_n-t_{n-1})\cdots {\int }_0^{t_2}d{t}_1{{\rm Y}}_{1\left(2\right)}^{\left(0\right)}(t_2-t_1){\mathfrak{G}}_{1\left(2\right)}^{\left(1\right)}(r+r^{\prime }-2\mathcal{R};t_1),\hfill \end{array}\end{array}$

with

(29)$\begin{array}{c}\hfill {\mathfrak{G}}_{1\left(2\right)}^{\left(0\right)}(r;t)={\displaystyle \frac{1}{\sqrt{{\mathcal{D}}_{1\left(2\right)}{t}^{\gamma }}}}{\mathrm{H}}_{11}^{10}\left[{\displaystyle \frac{\left|r\right|}{\sqrt{{\mathcal{D}}_{1\left(2\right)}{t}^{\gamma }}}}{|}_{\left(0,1\right)}^{\left(1-{\displaystyle \frac{\gamma }{2}},{\displaystyle \frac{\gamma }{2}}\right)}\right],\end{array}$

(30)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{G}}_{1\left(2\right)}^{\left(1\right)}(r;t)& ={\displaystyle \frac{1}{4\pi \sqrt{{\mathcal{D}}_{1\left(2\right)}}}}\hfill \\ & \times {\int }_0^t{\displaystyle \frac{d{t}^{\prime }}{t^{\prime \gamma }}}{(t-t^{\prime })}^{{\displaystyle \frac{\gamma }{2}}-1}{\mathrm{E}}_{{\displaystyle \frac{\gamma }{2}},{\displaystyle \frac{\gamma }{2}}}\left(-\sqrt{{\displaystyle \frac{{\mathcal{D}}_{1\left(2\right)}}{{\mathcal{R}}^2}}}{(t-t^{\prime })}^{{\displaystyle \frac{\gamma }{2}}}\right){\mathrm{H}}_{11}^{10}\left[{\displaystyle \frac{\left|r\right|}{\sqrt{{\mathcal{D}}_{1\left(2\right)}{t}^{\prime \gamma }}}}{|}_{\left(0,1\right)}^{\left(1-\gamma ,{\displaystyle \frac{\gamma }{2}}\right)}\right],\hfill \end{array}\end{array}$

and

(31)$\begin{array}{c}\hfill {{\rm Y}}_{1\left(2\right)}^{\left(0\right)}\left(t\right)={\displaystyle \frac{1}{4\pi \sqrt{{\mathcal{D}}_{1\left(2\right)}}{\mathcal{R}}^2}}{\int }_0^{t}d{t}^{\prime }{\eta }_{1\left(2\right)}\left(t^{\prime }\right){(t-t^{\prime })}^{-{\displaystyle \frac{\gamma }{2}}}{\mathrm{E}}_{{\displaystyle \frac{\gamma }{2}},1-{\displaystyle \frac{\gamma }{2}}}\left(-\sqrt{{\displaystyle \frac{{\mathcal{D}}_{1\left(2\right)}}{{\mathcal{R}}^2}}}{(t-t^{\prime })}^{{\displaystyle \frac{\gamma }{2}}}\right).\end{array}$

Finally, for the third case, we have

(32)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{G}}_{1\left(2\right)}(r,r^{\prime };t)& ={\displaystyle \frac{1}{8\pi r{r}^{\prime }}}\left\{{\mathfrak{G}}_{1\left(2\right)}^{\left(2\right)}(r-r^{\prime };t)-{\mathfrak{G}}_{1\left(2\right)}^{\left(2\right)}(r+r^{\prime }-2\mathcal{R};t)\right\}\hfill \\ & +{\displaystyle \frac{1}{r{r}^{\prime }}}{\mathfrak{G}}_{1\left(2\right)}^{\left(3\right)}(r+r^{\prime }-2\mathcal{R};t)+{\displaystyle \frac{1}{r{r}^{\prime }}}\sum _{n=1}^{\infty }{(-1)}^n{\int }_0^{t}d{t}_n{{\rm Y}}_{1\left(2\right)}^{\left(1\right)}(t-t_n)\hfill \\ & \times {\int }_0^{t_n}d{t}_{n-1}{{\rm Y}}_{1\left(2\right)}^{\left(1\right)}(t_n-t_{n-1})\cdots {\int }_0^{t_2}d{t}_1{{\rm Y}}_{1\left(2\right)}^{\left(1\right)}(t_2-t_1){\mathfrak{G}}_{1\left(2\right)}^{\left(3\right)}(r+r^{\prime }-2\mathcal{R};t_1),\hfill \end{array}\end{array}$

with

(33)$\begin{array}{c}\hfill {\mathfrak{G}}_{1\left(2\right)}^{\left(2\right)}(r;t)={\displaystyle \frac{e^{-\frac{\left|r\right|}{\sqrt{{\mathcal{D}}_{1\left(2\right)}t}}}}{\sqrt{{\mathcal{D}}_{1\left(2\right)}t}}}+{\displaystyle \frac{1}{\tau }}{\int }_0^{t}d{t}^{\prime }{e}^{-{\displaystyle \frac{t^{\prime }}{\tau }}}{\displaystyle \frac{e^{-\frac{\left|r\right|}{\sqrt{{\mathcal{D}}_{1\left(2\right)}{t}^{\prime }}}}}{\sqrt{{\mathcal{D}}_{1\left(2\right)}{t}^{\prime }}}},\end{array}$

(34)$\begin{array}{c}\hfill {\mathfrak{G}}_{1\left(2\right)}^{\left(3\right)}(r;t)={\mathfrak{G}}_{1\left(2\right)}^{\left(4\right)}(r+r^{\prime }-2\mathcal{R};t)+{\displaystyle \frac{1}{\tau }}{\int }_0^{t}d{t}^{\prime }{e}^{-{\displaystyle \frac{t}{\tau }}}{\mathfrak{G}}_{1\left(2\right)}^{\left(4\right)}(r+r^{\prime }-2\mathcal{R};t^{\prime }),\end{array}$

and

(35)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{G}}_{1\left(2\right)}^{\left(4\right)}(r;t)& =r{e}^{-{\displaystyle \frac{t}{\tau }}}{\int }_0^{t}d{t}^{\prime }{\Xi }^{\left(0\right)}(t-t^{\prime }){\displaystyle \frac{e^{-\frac{r^2}{4{\mathcal{D}}_{1\left(2\right)}{t}^{\prime }}}}{\sqrt{4\pi {\mathcal{D}}_{1\left(2\right)}{t}^{\prime 3}}}},\hfill \end{array}\end{array}$

in which

${\Xi }^{\left(0\right)}\left(t\right)={\displaystyle \frac{1}{4\pi \sqrt{{\mathcal{D}}_{1\left(2\right)}}}}\left\{{\displaystyle \frac{1}{\sqrt{\pi }t}}-\left[e^{{\mathcal{D}}_{1\left(2\right)}t/{\mathcal{R}}^2}\sqrt{{\displaystyle \frac{{\mathcal{D}}_{1\left(2\right)}}{{\mathcal{R}}^2}}}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\sqrt{{\displaystyle \frac{{\mathcal{D}}_{1\left(2\right)}t}{{\mathcal{R}}^2}}}\right)\right]\right\}.$

In addition, in Equation (32) the following quantities are present:

(36)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {{\rm Y}}_{1\left(2\right)}^{\left(1\right)}\left(t\right)& ={\displaystyle \frac{1}{{\mathcal{R}}^2}}{\int }_0^{t}d{t}^{\prime }{\eta }_{1\left(2\right)}\left(t^{\prime }\right)e^{-(t-t^{\prime })/\tau }{\Xi }^{\left(0\right)}(t-t^{\prime })\hfill \\ & +{\displaystyle \frac{1}{\tau {\mathcal{R}}^2}}{\int }_0^{t}d{t}^{\prime }{\int }_0^{t^{\prime }}d{t}^{″}{\eta }_{1\left(2\right)}\left(t^{″}\right)e^{-(t^{\prime }-t^{″})/\tau }{\Xi }^{\left(0\right)}(t^{\prime }-t^{″}).\hfill \end{array}\end{array}$

The last case allows us to obtain the stationary behavior for $t\to \infty$ in the case ${\eta }_{1\left(2\right)}\left(t\right)\to {\eta }_{1\left(2\right)}^{\left(st\right)}$ for $t\to \infty$, namely,

(37)$\begin{array}{c}\hfill {\widehat{\mathcal{G}}}_{1\left(2\right)}^{\left(st\right)}(r,r^{\prime })={\displaystyle \frac{1}{8\pi r{r}^{\prime }\sqrt{{\mathcal{D}}_{1\left(2\right)}}}}\left(e^{-|r-r^{\prime }|/\sqrt{{\mathcal{D}}_{1\left(2\right)}}}-e^{-|r+r^{\prime }-2\mathcal{R}|/\sqrt{{\mathcal{D}}_{1\left(2\right)}}}\right).\end{array}$

After some calculations in the Laplace domain, we obtain from these equations the following results:

(38)$\begin{array}{c}\hfill {\widehat{\rho }}_1(\mathcal{R},s)=\widehat{\Lambda }\left(s\right)\left[{\widehat{\mathcal{I}}}_1\left(s\right)+{\widehat{k}}_{s,12}\left(s\right){\widehat{\mathcal{I}}}_2\left(s\right){\widehat{\mathcal{G}}}_1(\mathcal{R},\mathcal{R},s)\right]\end{array}$

and

(39)$\begin{array}{c}\hfill {\widehat{\rho }}_2(\mathcal{R},s)=\widehat{\Lambda }\left(s\right)\left[{\widehat{\mathcal{I}}}_2\left(s\right)+{\widehat{k}}_{s,21}\left(s\right){\widehat{\mathcal{I}}}_1\left(s\right){\widehat{\mathcal{G}}}_2(\mathcal{R},\mathcal{R},s)\right],\end{array}$

where

(40)$\begin{array}{c}\hfill \widehat{\Lambda }\left(s\right)={\displaystyle \frac{1}{1-{\widehat{k}}_{s,21}\left(s\right){\widehat{k}}_{s,12}\left(s\right){\widehat{\mathcal{G}}}_1(\mathcal{R},\mathcal{R},s){\widehat{\mathcal{G}}}_2(\mathcal{R},\mathcal{R},s)}},\end{array}$

with ${\widehat{\mathcal{I}}}_{1\left(2\right)}\left(s\right)=4\pi {\int }_0^{\infty }d{r}^{\prime }{r^{\prime }}^2{\widehat{\mathcal{G}}}_{1\left(2\right)}(\mathcal{R},r^{\prime },s){\phi }_{1\left(2\right)}\left(r^{\prime }\right)$. These equations yield the survival probability

(41)$\begin{array}{c}\hfill S_{1\left(2\right)}\left(t\right)=4\pi {\int }_{\mathcal{R}}^{\infty }dr{r}^2{\rho }_{1\left(2\right)}(r,t),\end{array}$

which is related to the quantity of substance present in the bulk. After some calculations, we obtain, for each species,

(42)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\mathcal{S}}}_{1\left(2\right)}\left(s\right)& ={\displaystyle \frac{1}{s}}4\pi {\int }_{\mathcal{R}}^{\infty }d\tilde{r}\tilde{r}{\phi }_{1\left(2\right)}\left(\tilde{r}\right)\left(\tilde{r}-\mathcal{R}{e}^{-\sqrt{{\displaystyle \frac{s}{{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}}}(\tilde{r}-\mathcal{R})}\right)\hfill \\ & +{\displaystyle \frac{1}{s}}{\displaystyle \frac{4\pi {\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)\mathcal{R}\left(\mathcal{R}\sqrt{s/{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}+1\right)}{4\pi {\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)\mathcal{R}\left(\mathcal{R}\sqrt{s/{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}+1\right)+{\widehat{\eta }}_{1\left(2\right)}\left(s\right)}}\hfill \\ & \times \left({\widehat{k}}_{s,12\left(21\right)}\left(s\right){\rho }_{2\left(1\right)}(\mathcal{R},s)+4\pi \mathcal{R}{\int }_{\mathcal{R}}^{\infty }d\tilde{r}\tilde{r}{\phi }_{1\left(2\right)}\left(\tilde{r}\right)e^{-\sqrt{{\displaystyle \frac{s}{{\widehat{\mathcal{D}}}_{1\left(2\right)}\left(s\right)}}}(\tilde{r}-\mathcal{R})}\right).\hfill \end{array}\end{array}$

To proceed, let us analyze some special cases contained in the previous calculations. The first we mention is obtained by considering ${\kappa }_1\left(t\right)=0$ and ${\kappa }_2\left(t\right)=0$; i.e., the situation that corresponds to the absence of adsorption process on the surface, with $k_{s,11}\left(t\right)=k_{s,21}\left(t\right)=k_1\delta \left(t\right)$ and $k_{s,22}\left(t\right)=k_{s,12}\left(t\right)=k_2\delta \left(t\right)$. This case implies that the surface absorbs particles 1 and 2 to promote the reaction process $1\rightleftarrows 2$, where each substance promotes the formation of the other, then is released to the bulk. This reaction process on the surface is typical of a reversible reaction. In this representative case, the distributions related to each species can be written as follows:

(43)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\rho }}_1(r,s)& =4\pi {\int }_{\mathcal{R}}^{\infty }d{r}^{\prime }{r\prime }^2{\overline{\mathcal{G}}}_1(r,r^{\prime },s){\phi }_1\left(r^{\prime }\right)\hfill \\ & +{\mathcal{G}}_1(r,\mathcal{R};s){\displaystyle \frac{4\pi \mathcal{R}{k}_2\left[\widehat{\Phi }\left(s\right)+k_1\right]}{\widehat{\Phi }\left(s\right)\left[\widehat{\Phi }\left(s\right)+(k_1+k_2)\right]}}{\int }_{\mathcal{R}}^{\infty }d{r}^{\prime }{r}^{\prime }{\phi }_2\left(r^{\prime }\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})}\hfill \\ & +{\mathcal{G}}_1(r,\mathcal{R};s){\displaystyle \frac{4\pi \mathcal{R}{k}_2{k}_1}{\widehat{\Phi }\left(s\right)\left[\widehat{\Phi }\left(s\right)+(k_1+k_2)\right]}}{\int }_{\mathcal{R}}^{\infty }d{r}^{\prime }{r}^{\prime }{\phi }_1\left(r^{\prime }\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})},\hfill \end{array}\end{array}$

and

(44)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\rho }}_2(r,s)& =4\pi {\int }_{\mathcal{R}}^{\infty }d{r}^{\prime }{r\prime }^2{\widehat{\mathcal{G}}}_2(r,r^{\prime },s){\phi }_2\left(r^{\prime }\right)\hfill \\ & +{\mathcal{G}}_2(r,\mathcal{R};s){\displaystyle \frac{4\pi \mathcal{R}{k}_1\left[\widehat{\Phi }\left(s\right)+k_2\right]}{\widehat{\Phi }\left(s\right)\left[\widehat{\Phi }\left(s\right)+(k_1+k_2)\right]}}{\int }_{\mathcal{R}}^{\infty }d{r}^{\prime }{r}^{\prime }{\phi }_1\left(r^{\prime }\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})}\hfill \\ & +{\mathcal{G}}_2(r,\mathcal{R};s){\displaystyle \frac{4\pi \mathcal{R}{k}_1{k}_2}{\widehat{\Phi }\left(s\right)\left[\widehat{\Phi }\left(s\right)+(k_1+k_2)\right]}}{\int }_{\mathcal{R}}^{\infty }d{r}^{\prime }{r}^{\prime }{\phi }_2\left(r^{\prime }\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})},\hfill \end{array}\end{array}$

with $\widehat{\Phi }\left(s\right)=4\pi \widehat{\mathcal{D}}\left(s\right)\mathcal{R}\left(\mathcal{R}\sqrt{s/\widehat{\mathcal{D}}\left(s\right)}+1\right)$ and ${\widehat{\mathcal{D}}}_1\left(s\right)={\widehat{\mathcal{D}}}_2\left(s\right)=\widehat{\mathcal{D}}\left(s\right)$.

The survival probability, related to the quantity of each species of particle present in the bulk, in the case being considered are

(45)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\mathcal{S}}}_1\left(s\right)={\displaystyle \frac{1}{s}}4\pi {\int }_{\mathcal{R}}^{\infty }d\tilde{r}{\tilde{r}}^2{\phi }_1\left(\tilde{r}\right)& +{\displaystyle \frac{1}{s}}{\displaystyle \frac{4\pi \mathcal{R}{k}_2}{\widehat{\Phi }\left(s\right)+k_1+k_2}}{\int }_{\mathcal{R}}^{\infty }d\tilde{r}{\phi }_2\left(\tilde{r}\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})}\hfill \\ & -{\displaystyle \frac{1}{s}}{\displaystyle \frac{4\pi \mathcal{R}{k}_1}{\widehat{\Phi }\left(s\right)+k_1+k_2}}{\int }_{\mathcal{R}}^{\infty }d\tilde{r}{\phi }_1\left(\tilde{r}\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})}\hfill \end{array}\end{array}$

and

(46)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\mathcal{S}}}_2\left(s\right)={\displaystyle \frac{1}{s}}4\pi {\int }_{\mathcal{R}}^{\infty }d\tilde{r}{\tilde{r}}^2{\phi }_1\left(\tilde{r}\right)& -{\displaystyle \frac{1}{s}}{\displaystyle \frac{4\pi \mathcal{R}{k}_2}{\widehat{\Phi }\left(s\right)+k_1+k_2}}{\int }_{\mathcal{R}}^{\infty }d\tilde{r}{\phi }_2\left(\tilde{r}\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})}\hfill \\ & +{\displaystyle \frac{1}{s}}{\displaystyle \frac{4\pi \mathcal{R}{k}_1}{\widehat{\Phi }\left(s\right)+k_1+k_2}}{\int }_{\mathcal{R}}^{\infty }d\tilde{r}{\phi }_1\left(\tilde{r}\right)e^{-\sqrt{{\displaystyle \frac{s}{\widehat{\mathcal{D}}\left(s\right)}}}(\tilde{r}-\mathcal{R})}.\hfill \end{array}\end{array}$

From these expressions, ${\widehat{\mathcal{S}}}_1\left(s\right)+{\widehat{\mathcal{S}}}_2\left(s\right)=1/s$, and consequently, ${\mathcal{S}}_1\left(t\right)+{\mathcal{S}}_2\left(t\right)=1$.

Using these results makes it possible to obtain the first-passage-time distribution related to each species. The first-passage-time distribution is connected to the mean first passage time, which is the mean time taken for the system, in our case, a species, to reach a certain value. It is defined as

(47)$\begin{array}{c}\hfill {\mathcal{F}}_{1\left(2\right)}\left(t\right)=-{\displaystyle \frac{\partial }{\partial t}}\left\{4\pi {\int }_{\mathcal{R}}^{\infty }dr{r}^2{\rho }_{1\left(2\right)}(r,t)\right\}=-{\displaystyle \frac{\partial }{\partial t}}{\mathcal{S}}_{1\left(2\right)}\left(t\right).\end{array}$

The next move is to provide numerical simulations from the stochastic equations’ perspective. To proceed this way, we handle the discrete form of the Langevin equations as follows:

(48)$\begin{array}{c}\hfill \begin{array}{cc}\hfill x_{t+h}& =x_t+\sqrt{2\mathcal{D}h}{\zeta }_x\left(t\right),\hfill \\ \hfill y_{t+h}& =y_t+\sqrt{2\mathcal{D}h}{\zeta }_y\left(t\right),\hfill \\ \hfill z_{t+h}& =z_t+\sqrt{2\mathcal{D}h}{\zeta }_z\left(t\right).\hfill \end{array}\end{array}$

For the whole numerical simulation, we consider the Cartesian coordinates with $r^2\left(t\right)=x^2\left(t\right)+y^2\left(t\right)+z^2\left(t\right)$, and take into account ${\zeta }_i\left(t\right)$ ($i=x,y,$ and z) as white Gaussian noise, with a normalized deviation generated by the Box–Muller method [56,57] (see the Appendix A). Furthermore, its mean value is zero $(⟨{\zeta }_i\left(t\right)⟩=0)$, $⟨{\zeta }_i\left(t\right){\zeta }_j\left(t^{\prime }\right)⟩=0$, for $i\ne j$, and $⟨{\zeta }_i\left(t\right){\zeta }_i\left(t^{\prime }\right)⟩\propto \delta (t-t^{\prime })$. We also consider, without loss of generality, the same diffusion coefficient for the particles of species 1 and 2, i.e., ${\mathcal{D}}_1={\mathcal{D}}_2=\mathcal{D}$. Note that the set of stochastic equations present in Equation (48) corresponds to the standard Brownian motion, i.e., usual diffusion. In this sense, we use Equation (25) for the Green’s function to obtain the solution and perform the comparison between the analytical approach and the results obtained from the Langevin equation.

The absorption process of the particles by the surface is treated stochastically. When a particle is absorbed by the surface, it promotes the formation of a particle of a different species and vice versa; they are released back to the bulk as soon as the reaction occurs. Concerning absorption, we consider the probability $p_{12}h$ for species 1 and $p_{21}h$ for species 2, with $p_{12}>p_{21}$. In the case of total absorption of particles by the surface, in the absence of the reaction process, the absorption probability is $ph$. In both scenarios, absorption can only occur when the particle touches the surface, i.e., when $r\le R$ is verified. Figure 5 illustrates the diffusion process obtained with the set of Langevin equations given by Equation (48), with the surface absorbing the particles, and a reaction process promoting the formation of particles of different species.

Figure 6 and Figure 7 show the numerical simulations and the analytical results. We observe a good agreement between these approaches in both scenarios, that is, in the case of the total absorption process of the particles by the surface and in the case of absorption with the reaction process on the surface, with the release to the bulk of the particles originating from the reaction process. We also consider, in the first scenario, the first-passage-time distribution with good agreement between the results obtained in the numerical and analytical approaches.

4. Discussion

We have investigated diffusion processes exhibiting spherical symmetry, representing a colloidal particle. The surface can either adsorb–desorb, fully absorb particles, or catalyze reactions that promote particle formation when the surface absorbs particles. The processes on the bulk and surface are coupled by the boundary conditions, which implies that the dynamic processes on the surface influence the bulk dynamics and vice versa. This point is evidenced by the boundary conditions, which couple the kinetics on the surface with the bulk dynamics, i.e., the flux of particles nearby. This feature changes the bulk dynamics and may introduce different diffusion processes, as shown in Refs. [53,58]. The particles obtained from the reaction process are released into the bulk. The adsorbed particles are also desorbed in the bulk after some characteristic time, defined by a relaxation process, which can be of the Debye or non-Debye type depending on the kinetic equations considered for the surface effects. The non-Debye relaxation depends on the kernel choice in the convolution integrals, which define the boundary conditions. These kernels may also be used in scenarios with intermediate reaction processes on the surface, where the formation of one species depends on an intermediate reaction during the kinetic process. We also considered a generalized diffusion equation to govern the dynamics of the particles in the bulk. In particular, we considered a fractional differential operator that can be connected to different scenarios, such as the Riemann–Liouville fractional operator (singular kernel), Caputo–Fabrizio fractional operator (nonsingular kernel), Atangana–Baleanu fractional operator (nonsingular kernel), and others. Note that this equation, as discussed in Section 2, may be connected to a random walk, which depends on the characteristics of the media. These characteristics are manifested by the choice of the probability density function associated with the diffusive dynamics of the species present in the system. Thus, formulated in this manner, this approach can describe a rich class of diffusion processes, either Markovian or non-Markovian.

Non-Markovian diffusion is usually characterized, for example, by memory effects, long-range correlations, and intermittent dynamics. From the point of view of the random walk formulation, these scenarios associated with non-Markovian processes lead us to obtain, for example, a long-tailed behavior for the waiting time distributions, which is different from the usual one. In the Laplace domain, we found a general solution for the diffusion equation by taking into account the boundary conditions given by Equations (11) and (12), which is represented by Equations (18) and (19) for species 1 and 2, with the Green function given by Equation (23). We chose singular and nonsingular kernels related to some representative fractional differential operators to determine the inverse Laplace transform. The results for each case showed how the fractional operator influences solutions, and consequently, the bulk dynamics. In particular, the solution exhibited stationary behavior in the case of a nonsingular kernel.

On the other hand, we also carried out some numerical simulations using stochastic equations, i.e., Langevin equations, to simulate the random motion of the particles in the bulk. The absorption process was defined as stochastic, with a given probability of absorption when interacting with the spherical surface. After that, we considered that the absorbed particles can promote the formation of particles of different species by a determined stochastic process. The results obtained from the numerical simulation were compared with analytical results for the survival probability (particles in bulk) and first-passage-time distribution, as shown in Figure 6 and Figure 7, with a good agreement.

5. Conclusions

The system investigated has shown that the dynamics of particles are influenced by surface and bulk effects. The surface effects considered here are connected with the adsorption–desorption processes or reactions with the formation of another species, which is released to the bulk. The influence of the bulk on the particle dynamics is connected to the choice of the kernel present in the fractional operator, which reflects how the media changes the spreading of the system by introducing effects that are present in standard diffusion processes. By combining these effects (surface and bulk), a large class of diffusion processes may be investigated, particularly ones related to anomalous diffusion. This approach can also be used to investigate the effect of systems with different regimes of diffusion on the choice of the kernel, see Refs. [50,59], or connected to other processes such as stochastic resetting [42]. Finally, we hope that it can be useful in discussing surface effects coupled with bulk dynamics in confined or restricted geometries.

Author Contributions

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

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

We thank Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES—Brazil), National Council for Scientific and Technological Development (CNPq—Brazil), and Politecnico di Torino (Italy).

Conflicts of Interest

The authors declare no conflicts 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. This figure illustrates particles (species 1 and 2 (orange and purple)) near or surrounding the spherical colloidal particle in a colloidal suspension. The particles are at a radial position r interacting with a colloidal particle with a spherical surface of radius R.

Figure 2. Trend of the Green’s function for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for different values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] with [Forumla omitted. See PDF.]. We consider, for illustrative purposes, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], in arbitrary units.

Figure 3. Trend of the Green’s function for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], for different values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. We consider [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], in arbitrary units.

Figure 4. The profile of the Green’s function for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], for different values of t. For illustrative purposes, we consider [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], in arbitrary units.

Figure 5. The probability density maps in units of particles per [Forumla omitted. See PDF.] obtained from the numerical simulation of Langevin equations when the surface may adsorb the particles and, by a reaction process, promote the formation of particles of different species. In (a), particles of species 1, and in (b) particles of species 2. The particles are created by a reaction process at [Forumla omitted. See PDF.] and, after that, they are released into the bulk to diffuse. The parameters used, for illustrative purposes, are [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] (initial position), [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 6. Comparison of the adsorption ([Forumla omitted. See PDF.]) and the first-passage-time distribution ([Forumla omitted. See PDF.]) from the analytical and numerical simulations for different values of k. For the analytical approach, we use [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. For the reaction term, we consider [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.]. The curves are drawn for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], in arbitrary units. For the simulations, [Forumla omitted. See PDF.] for all systems, then the absorption probabilities are [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] from the bottom to the top curve.

Figure 7. Comparison of the survival probability for the two species (1 and 2) obtained from the numerical and analytical calculations for different values of [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] and 100. For the analytical approach, we use [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. For the reaction term, we consider [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.]. The absorption constants are related as follows: [Forumla omitted. See PDF.], where for each [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. We use the values [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], in arbitrary units, for illustrative purposes. For the simulations, [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for all systems.

Appendix A. Numerical Approach

In this section, we show in more detail how the numerical simulations were carried out. The programming language used was C++, with dependence on the pseudo-random number generator (RNG) “maximally equidistributed combined Tausworthe generator” [60] tauss88 from the library Boost [61].

The code shown is for a system with two types of particles. We begin with diffusion in x, y, and z:

D=sqrt(2.0*D*h); //for all steps to be simulated.  for(n=0;n<steps;n++){ //n1 and n2 are the numbers of particles of type one and type two.  n1=0;n2=0; //for all particles in the system  for(i=0;i<npar;i++){ //dist is the distance from the origin or r.  //ext exists so if n1 becomes n2 then we don’t check for reaction again.  dist=0;ext=1; //this will compute diffusion in x,y, and z  for(j=0;j<3;j++){ //double **pxyz (pxyz[npar][3]) is the position matrix  //the funcion gauss() creates a gaussian number via box-muller  pxyz[i][j]+=D*gauss(); //x^2+y^2+z^2  dist+=pxyz[i][j]*pxyz[i][j]; }

After we have computed the diffusion and $r^2$ for a single particle, we can check for interaction with the surface.

//if the particle is type 1  if(calc[i]==0){ n1++; //if the particle is inside or on the surface (since h > 0)  if(sqrt(dist)<=R){ // a is how much you need to multiply x,y, and z  //so that the particle is just above the surface  a=((R+0.01)/sqrt(dist)); for(j=0;j<3;j++){ pxyz[i][j]*=a; } //luck is the RNG/max(RNG) therefore 0<=luck<=1  luck=1.0*generator()/(1.0*boost::taus88().max()); if(luck<p12*h){ //particle n1 becomes n2  calc[i]=1; // so we don’t check for absortion again  ext=0; n1--; n2++; } } } //n2 follows the same ideia from n1:  if(calc[i]==1 & ext){ n2++; if(sqrt(dist)<=R){ a=((R+0.01)/sqrt(dist)); for(j=0;j<3;j++){ pxyz[i][j]*=a; } luck=1.0*generator()/(1.0*boost::taus88().max()); if(luck<p21*h){ calc[i]=0; n2--; n1++; } } } }

The dynamics of the system are established; now, we can generate Gaussian noise with two linearly distributed random numbers $u_1$ and $u_2$, where (A1)$$\begin{array}{c}\hfill \begin{array}{cc}& 0<u_1\le 1\hfill \\ & 0\le u_2\le 1\hfill \\ & Z_0=\sqrt{-2.0\times ln\left(u_1\right)}\times cos\left(2.0\pi u_2\right),\hfill \\ & Z_1=\sqrt{-2.0\times ln\left(u_1\right)}\times sin\left(2.0\pi u_2\right).\hfill \end{array}\end{array}$$

Both $Z_0$ and $Z_1$ are Gaussian-distributed numbers; see Figure A1. In the code, only $Z_0$ was used in the function “gauss”:

double gauss(){ double u1,u2,val; u1=1.0*generator()/(1.0*boost::taus88().max()); u2=1.0*generator()/(1.0*boost::taus88().max()); while(u1==0.0){ u1=1.0*generator()/(1.0*boost::taus88().max()); } val=sqrt(-2.0*log(u1))*cos(2.0*PI*u2); return(val); }

Figure A1. Probability density distribution of the Gaussian [Forumla omitted. See PDF.] generated via the Box–Muller method.

References

1. Jones, R. Soft Condensed Matter; Oxford Master Series in Physics; OUP Oxford: Oxford, UK, 2002.

2. Rosenberg, R.T.; Dan, N. Self-Assembly of Colloidosome Shells on Drug-Containing Hydrogels. J. Biomater. Nanobiotechnol.; 2011; 2, pp. 1-7. [DOI: https://dx.doi.org/10.4236/jbnb.2011.21001]

3. Smith, G.N.; Eastoe, J. Controlling colloid charge in nonpolar liquids with surfactants. Phys. Chem. Chem. Phys.; 2013; 15, pp. 424-439. [DOI: https://dx.doi.org/10.1039/C2CP42625K] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23187453]

4. Comiskey, B.; Albert, J.D.; Yoshizawa, H.; Jacobson, J. An electrophoretic ink for all-printed reflective electronic displays. Nature; 1998; 394, pp. 253-255. [DOI: https://dx.doi.org/10.1038/28349]

5. Mereghetti, P.; Kokh, D.; McCammon, J.A.; Wade, R.C. Diffusion and association processes in biological systems: Theory, computation and experiment. BMC Biophys.; 2011; 4, 2. [DOI: https://dx.doi.org/10.1186/2046-1682-4-2]

6. Benelli, R.; Weiss, M. From sub- to superdiffusion: Fractional Brownian motion of membraneless organelles in early C. elegans embryos. New J. Phys.; 2021; 23, 063072. [DOI: https://dx.doi.org/10.1088/1367-2630/ac0853]

7. Brangwynne, C.P.; Koenderink, G.H.; MacKintosh, F.C.; Weitz, D.A. Intracellular transport by active diffusion. Trends Cell Biol.; 2009; 19, pp. 423-427. [DOI: https://dx.doi.org/10.1016/j.tcb.2009.04.004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19699642]

8. Bechinger, C.; Di Leonardo, R.; Löwen, H.; Reichhardt, C.; Volpe, G.; Volpe, G. Active particles in complex and crowded environments. Rev. Mod. Phys.; 2016; 88, 045006. [DOI: https://dx.doi.org/10.1103/RevModPhys.88.045006]

9. Höfling, F.; Franosch, T. Anomalous transport in the crowded world of biological cells. Rep. Prog. Phys.; 2013; 76, 046602. [DOI: https://dx.doi.org/10.1088/0034-4885/76/4/046602]

10. Javanainen, M.; Hammaren, H.; Monticelli, L.; Jeon, J.H.; Miettinen, M.S.; Martinez-Seara, H.; Metzler, R.; Vattulainen, I. Anomalous and normal diffusion of proteins and lipids in crowded lipid membranes. Faraday Discuss.; 2013; 161, pp. 397-417. [DOI: https://dx.doi.org/10.1039/C2FD20085F]

11. Sokolov, I.M. Models of anomalous diffusion in crowded environments. Soft Matter; 2012; 8, pp. 9043-9052. [DOI: https://dx.doi.org/10.1039/c2sm25701g]

12. Elowitz, M.B.; Surette, M.G.; Wolf, P.E.; Stock, J.B.; Leibler, S. Protein mobility in the cytoplasm of Escherichia coli. J. Bacteriol.; 1999; 181, pp. 197-203. [DOI: https://dx.doi.org/10.1128/JB.181.1.197-203.1999] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9864330]

13. Kopelman, R. Fractal reaction kinetics. Science; 1988; 241, pp. 1620-1626. [DOI: https://dx.doi.org/10.1126/science.241.4873.1620] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17820893]

14. Seri-Levy, A.; Avnir, D. Kinetics of diffusion-limited adsorption on fractal surfaces. J. Phys. Chem.; 1993; 97, pp. 10380-10384. [DOI: https://dx.doi.org/10.1021/j100142a019]

15. Giona, M.; Giustiniani, M. Adsorption kinetics on fractal surfaces. J. Phys. Chem.; 1996; 100, pp. 16690-16699. [DOI: https://dx.doi.org/10.1021/jp961518l]

16. Kosztołowicz, T. Random walk model of subdiffusion in a system with a thin membrane. Phys. Rev. E; 2015; 91, 022102. [DOI: https://dx.doi.org/10.1103/PhysRevE.91.022102]

17. Kosztołowicz, T. Subdiffusion in a system consisting of two different media separated by a thin membrane. Int. J. Heat Mass Transf.; 2017; 111, pp. 1322-1333. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.04.058]

18. Kosztołowicz, T.; Dutkiewicz, A. Boundary conditions at a thin membrane for the normal diffusion equation which generate subdiffusion. Phys. Rev. E; 2021; 103, 042131. [DOI: https://dx.doi.org/10.1103/PhysRevE.103.042131]

19. Bisquert, J.; Compte, A. Theory of the electrochemical impedance of anomalous diffusion. J. Electroanal. Chem.; 2001; 499, pp. 112-120. [DOI: https://dx.doi.org/10.1016/S0022-0728(00)00497-6]

20. Evangelista, L.R.; Lenzi, E.K. Fractional Diffusion Equations and Anomalous Diffusion; Cambridge University Press: Cambridge, UK, 2018.

21. Metzler, R.; Glöckle, W.G.; Nonnenmacher, T.F. Fractional model equation for anomalous diffusion. Phys. A; 1994; 211, pp. 13-24. [DOI: https://dx.doi.org/10.1016/0378-4371(94)90064-7]

22. O’Shaughnessy, B.; Procaccia, I. Analytical Solutions for Diffusion on Fractal Objects. Phys. Rev. Lett.; 1985; 54, pp. 455-458. [DOI: https://dx.doi.org/10.1103/PhysRevLett.54.455]

23. Spiechowicz, J.; Łuczka, J. Subdiffusion via dynamical localization induced by thermal equilibrium fluctuations. Sci. Rep.; 2017; 7, 16451. [DOI: https://dx.doi.org/10.1038/s41598-017-16601-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29184075]

24. Bronaugh, R.L.; Stewart, R.F.; Congdon, E.R.; Giles, A.L., Jr. Methods for in vitro percutaneous absorption studies I. Comparison with in vivo results. Toxicol. Appl. Pharmacol.; 1982; 62, pp. 474-480. [DOI: https://dx.doi.org/10.1016/0041-008X(82)90148-X] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/7071862]

25. Abd, E.; Yousef, S.A.; Pastore, M.N.; Telaprolu, K.; Mohammed, Y.H.; Namjoshi, S.; Grice, J.E.; Roberts, M.S. Skin models for the testing of transdermal drugs. Clin. Pharmacol. Adv. Appl.; 2016; 8, pp. 163-176. [DOI: https://dx.doi.org/10.2147/CPAA.S64788]

26. Wang, L.; Cao, T.; Dykstra, J.E.; Porada, S.; Biesheuvel, P.; Elimelech, M. Salt and water transport in reverse osmosis membranes: Beyond the solution-diffusion model. Environ. Sci. Technol.; 2021; 55, pp. 16665-16675. [DOI: https://dx.doi.org/10.1021/acs.est.1c05649]

27. Soltanieh, M.; GILL’, W.N. Review of reverse osmosis membranes and transport models. Chem. Eng. Commun.; 1981; 12, pp. 279-363. [DOI: https://dx.doi.org/10.1080/00986448108910843]

28. Peters, R. Translational diffusion in the plasma membrane of single cells as studied by fluorescence microphotolysis. Cell Biol. Int. Rep.; 1981; 5, pp. 733-760. [DOI: https://dx.doi.org/10.1016/0309-1651(81)90231-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/6269761]

29. Qian, H.; Sheetz, M.P.; Elson, E.L. Single particle tracking. Analysis of diffusion and flow in two-dimensional systems. Biophys. J.; 1991; 60, pp. 910-921. [DOI: https://dx.doi.org/10.1016/S0006-3495(91)82125-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/1742458]

30. Zistler, M.; Wachter, P.; Wasserscheid, P.; Gerhard, D.; Hinsch, A.; Sastrawan, R.; Gores, H. Comparison of electrochemical methods for triiodide diffusion coefficient measurements and observation of non-Stokesian diffusion behaviour in binary mixtures of two ionic liquids. Electrochim. Acta; 2006; 52, pp. 161-169. [DOI: https://dx.doi.org/10.1016/j.electacta.2006.04.050]

31. Wang, P.; Anderko, A. Modeling self-diffusion in mixed-solvent electrolyte solutions. Ind. Eng. Chem. Res.; 2003; 42, pp. 3495-3504. [DOI: https://dx.doi.org/10.1021/ie030050n]

32. Schienbein, M.; Gruler, H. Langevin equation, Fokker-Planck equation and cell migration. Bull. Math. Biol.; 1993; 55, pp. 585-608. [DOI: https://dx.doi.org/10.1016/S0092-8240(05)80241-1]

33. Risken, H.; Risken, H. Fokker-Planck Equation; Springer: Berlin/Heidelberg, Germany, 1996.

34. Kenkre, V.M.; Montroll, E.W.; Shlesinger, M.F. Generalized master equations for continuous-time random walks. J. Stat. Phys.; 1973; 9, pp. 45-50. [DOI: https://dx.doi.org/10.1007/BF01016796]

35. Kenkre, V. The generalized master equation and its applications. Statistical Mechanics and Statistical Methods in Theory and Application; Springer: New York, NY, USA, 1977; pp. 441-461.

36. Pedron, I.T.; Mendes, R.; Malacarne, L.C.; Lenzi, E.K. Nonlinear anomalous diffusion equation and fractal dimension: Exact generalized Gaussian solution. Phys. Rev. E; 2002; 65, 041108. [DOI: https://dx.doi.org/10.1103/PhysRevE.65.041108] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12005807]

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

38. Metzler, R.; Klafter, J. The random walk’s guide to anomalous diffusion: A fractional dynamics approach. Phys. Rep.; 2000; 339, pp. 1-77. [DOI: https://dx.doi.org/10.1016/S0370-1573(00)00070-3]

39. Górska, K.; Horzela, A.; Pogány, T. Non-Debye relaxations: Smeared time evolution, memory effects, and the Laplace exponents. Commun. Nonlinear Sci. Numer. Simul.; 2021; 99, 105837. [DOI: https://dx.doi.org/10.1016/j.cnsns.2021.105837]

40. Górska, K.; Horzela, A.; Penson, K.A. Non-Debye Relaxations: The Ups and Downs of the Stretched Exponential vs. Mittag–Leffler’s Matchings. Fractal Fract.; 2021; 5, 265. [DOI: https://dx.doi.org/10.3390/fractalfract5040265]

41. Caputo, M.; Fabrizio, M. Applications of new time and spatial fractional derivatives with exponential kernels. Prog. Fract. Differ. Appl.; 2016; 2, pp. 1-11. [DOI: https://dx.doi.org/10.18576/pfda/020101]

42. Tateishi, A.A.; Ribeiro, H.V.; Lenzi, E.K. The Role of Fractional Time-Derivative Operators on Anomalous Diffusion. Front. Phys.; 2017; 5, 52. [DOI: https://dx.doi.org/10.3389/fphy.2017.00052]

43. Atangana, A.; Baleanu, D. New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model. arXiv; 2016; arXiv: 1602.03408

44. Baleanu, D.; Fernandez, A. On fractional operators and their classifications. Mathematics; 2019; 7, 830. [DOI: https://dx.doi.org/10.3390/math7090830]

45. Haider, S.S.; Rehman, M.u.; Abdeljawad, T. On Hilfer fractional difference operator. Adv. Differ. Equ.; 2020; 2020, 122. [DOI: https://dx.doi.org/10.1186/s13662-020-02576-2]

46. Kochubei, A.N. Equations with general fractional time derivatives–Cauchy problem. Fractional Differential Equations; De Gruyter: Berlin, Germany, 2019; pp. 223-234. [DOI: https://dx.doi.org/10.1515/9783110571660-011]

47. Luchko, Y.; Yamamoto, M. The general fractional derivative and related fractional differential equations. Mathematics; 2020; 8, 2115. [DOI: https://dx.doi.org/10.3390/math8122115]

48. Luchko, Y. General fractional integrals and derivatives with the Sonine kernels. Mathematics; 2021; 9, 594. [DOI: https://dx.doi.org/10.3390/math9060594]

49. Al-Refai, M.; Luchko, Y. General Fractional Calculus Operators of Distributed Order. Axioms; 2023; 12, 1075. [DOI: https://dx.doi.org/10.3390/axioms12121075]

50. Sandev, T.; Chechkin, A.; Kantz, H.; Metzler, R. Diffusion and Fokker-Planck-Smoluchowski equations with generalized memory kernel. Fract. Calc. Appl. Anal.; 2015; 18, pp. 1006-1038. [DOI: https://dx.doi.org/10.1515/fca-2015-0059]

51. Zhang, X.; Tong, D. A generalized Weber transform and its inverse formula. Appl. Math. Comput.; 2007; 193, pp. 116-126. [DOI: https://dx.doi.org/10.1016/j.amc.2007.03.054]

52. Dokoumetzidis, A.; Macheras, P. Fractional kinetics in drug absorption and disposition processes. J. Pharmacokinet. Pharmacodyn.; 2009; 36, pp. 165-178. [DOI: https://dx.doi.org/10.1007/s10928-009-9116-x]

53. Zola, R.S.; Lenzi, E.K.; Evangelista, L.R.; Barbero, G. Memory effect in the adsorption phenomena of neutral particles. Phys. Rev. E; 2007; 75, 042601. [DOI: https://dx.doi.org/10.1103/PhysRevE.75.042601]

54. Saxena, R.; Mathai, A.; Haubold, H. On fractional kinetic equations. Astrophys. Space Sci.; 2002; 282, pp. 281-287. [DOI: https://dx.doi.org/10.1023/A:1021175108964]

55. Saxena, R.; Mathai, A.; Haubold, H. On generalized fractional kinetic equations. Phys. A Stat. Mech. Its Appl.; 2004; 344, pp. 657-664. [DOI: https://dx.doi.org/10.1016/j.physa.2004.06.048]

56. Ökten, G.; Göncü, A. Generating low-discrepancy sequences from the normal distribution: Box–Muller or inverse transform?. Math. Comput. Model.; 2011; 53, pp. 1268-1281. [DOI: https://dx.doi.org/10.1016/j.mcm.2010.12.011]

57. Scott, D.W. Box–Muller transformation. WIREs Comput. Stat.; 2011; 3, pp. 177-179. [DOI: https://dx.doi.org/10.1002/wics.148]

58. Guimaraes, V.G.; Ribeiro, H.V.; Li, Q.; Evangelista, L.R.; Lenzi, E.K.; Zola, R.S. Unusual diffusing regimes caused by different adsorbing surfaces. Soft Matter; 2015; 11, pp. 1658-1666. [DOI: https://dx.doi.org/10.1039/C5SM00151J] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25633342]

59. 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]

60. L’Ecuyer, P. Maximally equidistributed combined Tausworthe generators. Math. Comput.; 1996; 65, pp. 203-213. [DOI: https://dx.doi.org/10.1090/S0025-5718-96-00696-5]

61. Available online: https://www.boost.org/ (accessed on 2 May 2024).