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., $\langle {(r-\langle r\rangle )}^2\rangle \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., $\langle {(r-\langle r\rangle )}^2\rangle \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:

$$\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}$$

(1)

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:

$$\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}$$

(2)

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

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

(3)

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:

$$\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}$$

(4)

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

$$\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}$$

(5)

and

$$\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}$$

(6)

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

$$\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}$$

(7)

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

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

(8)

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

$$\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}$$

(9)

By performing the inverse integral transforms, we obtain that

$$\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}$$

(10)

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

$$\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}$$

(11)

and

$$\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}$$

(12)

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

$$\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}$$

(13)

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

$$\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}$$

(14)

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.,

$$\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}$$

(15)

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

$$\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}$$

(16)

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.,

$$\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}$$

(17)

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

$$\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}$$

(18)

and

$$\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}$$

(19)

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

$$\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}$$

(20)

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

$$\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}$$

(21)

and

$$\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}$$

(22)

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

$$\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}$$

(23)

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:

• ${\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);
• ${\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);
• ${\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

$$\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}$$

(24)

with

$$\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}$$

(25)

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

$$\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}$$

(26)

where

$$\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}$$

(27)

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

$$\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}$$

(28)

with

$$\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}$$

(29)

$$\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}$$

(30)

and

$$\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}$$

(31)

Finally, for the third case, we have

$$\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}$$

(32)

with

$$\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}$$

(33)

$$\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}$$

(34)

and

$$\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}$$

(35)

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:

$$\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}$$

(36)

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,

$$\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}$$

(37)

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

$$\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}$$

(38)

and

$$\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}$$

(39)

where

$$\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}$$

(40)

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

$$\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}$$

(41)

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

$$\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}$$

(42)

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:

$$\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}$$

(43)

and

$$\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}$$

(44)

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

$$\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}$$

(45)

and

$$\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}$$

(46)

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

$$\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}$$

(47)

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:

$$\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}$$

(48)

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 $(\langle {\zeta }_i\left(t\right)\rangle =0)$, $\langle {\zeta }_i\left(t\right){\zeta }_j\left(t^{\prime }\right)\rangle =0$, for $i\ne j$, and $\langle {\zeta }_i\left(t\right){\zeta }_i\left(t^{\prime }\right)\rangle \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.