Non-Markovian Diffusion and Adsorption–Desorption Dynamics: Analytical and Numerical Results

Abstract

The interplay of diffusion with phenomena like stochastic adsorption–desorption, absorption, and reaction–diffusion is essential for life and manifests in diverse natural contexts. Many factors must be considered, including geometry, dimensionality, and the interplay of diffusion across bulk and surfaces. To address this complexity, we investigate the diffusion process in heterogeneous media, focusing on non-Markovian diffusion. This process is limited by a surface interaction with the bulk, described by a specific boundary condition relevant to systems such as living cells and biomaterials. The surface can adsorb and desorb particles, and the adsorbed particles may undergo lateral diffusion before returning to the bulk. Different behaviors of the system are identified through analytical and numerical approaches.

1. Introduction

Describing the physical world is complex and necessitates employing various methods, both experimental and theoretical, to grasp and interpret the behavior of systems. Diffusion is prevalent in various situations that can be usual or anomalous and may appear combined with different processes, such as adsorption–desorption [1] and reaction–diffusion [2]. In the case of normal diffusion, the system exhibits Markovian characteristics [3], particularly the mean square displacement with a linear dependence on time, i.e., $\langle {(r-\langle r\rangle )}^2\rangle \sim t$. In contrast, anomalous diffusion results in stochastic processes that govern the system exhibiting non-Markovian features [4], yielding a nonlinear relationship for the mean square displacement, e.g., $\langle {(r-\langle r\rangle )}^2\rangle \sim t^{\alpha }$, where $\alpha <1$ or $\alpha >1$ corresponds to sub- or superdiffusion, respectively [5,6]. In conjunction with these situations, adsorption–desorption or reaction processes may occur, directly impacting various aspects of the system. These phenomena play a relevant role in many scenarios, such as antibody binding and coupling to a receptor in a cell [7], electrical impedance [5,8], polymer dynamics at solid–liquid interfaces [9], molecule traveling through a cellular membrane [10], the dynamics of loci in a chromosome [11], the movement of a tracer particle [12], catalytic kinetics [13,14], and, in particular, in living systems [15,16,17]. The particles (or molecules) adsorbed by the surface can diffuse or become immobile and, after some time, can be desorbed to the bulk. In this manner, the system can also present lateral diffusion associated with the adsorbed particles. This behavior has been observed in different systems; for example, lateral diffusion occurs in a cellular membrane by lipids and proteins, which may occur in different modes, such as homogeneous or hop diffusion [18,19], and, in addition, on a tubular membrane, which plays a vital role in neuronal axons by transporting signaling molecules and proteins [20]. Lateral diffusion is not restricted to biological systems and is also found in graphene oxide sheets dispersed in solvents such as water [21], and charge carriers find their place in semiconductors through the lateral photoelectric effect [22].

This study investigates a heterogeneous diffusion process for a three-dimensional system subjected to an adsorption–desorption process on a surface. The diffusion equation governs the bulk particles with spatial dependence on the diffusion coefficient in the z-direction. The adsorption–desorption process occurs on the surface located at the point $z=0$. The surface is made up of the plane $(x,y)$, and the kinetic terms explain the interaction with the particles in the bulk. These terms represent the adsorption–desorption process of the particles from the bulk to the surface (adsorption) and from the surface to the bulk (desorption). Adsorbed particles can diffuse on the surface and, after some time, are desorbed to the bulk. However, some particles will be desorbed in the next step since the process is stochastic. We consider that the particles can diffuse on the surface prior to desorption. We only consider inhomogeneity in the bulk, represented by a spatial dependence on the diffusion coefficient in the z-direction. These features imply a coupling between the processes in the bulk and on the surface, where each process influences the other. This analysis explores the relationship between bulk processes and surface processes, particularly examining the role of bulk heterogeneity in regulating surface behavior. This is performed in Section 2, where the analytical and numerical calculations are presented and discussed. In Section 3, we discuss the results and offer some conclusions about the behavior exhibited in the problems analyzed here.

2. The Problem: Diffusion and Kinetics

Let us initiate our analysis by considering that the subsequent equation controls the diffusion process in the bulk:

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_{\mathrm{bulk}}(\mathbf{r},t)={\mathcal{D}}_{b,\Vert }^{ }{\nabla }_{\Vert }^2{\rho }_{\mathrm{bulk}}(\mathbf{r},t)+{\nabla }_{\perp }·\left[{\mathcal{D}}_{b,\perp }\left(\mathbf{r}\right){\nabla }_{\perp }{\rho }_{\mathrm{bulk}}(\mathbf{r},t)\right].\end{array}$$

(1)

Equation (1) is the typical continuity equation applied to Fick’s law, i.e.,

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\rho }_{\mathrm{bulk}}(\mathbf{r},t)+\nabla ·{\mathcal{J}}_{\mathrm{bulk}}(\mathbf{r},t)=0\end{array}$$

(2)

combined with ${\mathcal{J}}_{\mathrm{bulk}}(\mathbf{r},t)$=$-{\mathcal{D}}_{b,\Vert }{\nabla }_{\Vert }{\rho }_{\mathrm{bulk}}(\mathbf{r},t)-{\mathcal{D}}_{b,\perp }\left(\mathbf{r}\right){\nabla }_{\perp }{\rho }_{\mathrm{bulk}}(\mathbf{r},t)$, which states that the time variation in the density of particles in the bulk is equal to the spatial change caused by concentration gradients. We assume that the spatial term is composed by a diffusion process taking place in the bulk, where there is heterogeneity in the $\widehat{\mathbf{z}}$ direction, but particles can diffuse in 3D. More specifically, ${\mathcal{D}}_{b,\perp }\left(\mathbf{r}\right)={\mathcal{D}}_{b,\perp }{r}_{\perp }^{-\eta }$ ($-1\le \eta$ gives the degree of heterogeneity for the z direction), ${\nabla }_{\perp }·\left[r_{\perp }^{-\eta }{\nabla }_{\perp }\rho (\mathbf{r},t)\right]={\partial }_z\left[{\left|z\right|}^{-\eta }{\partial }_z\rho (\mathbf{r},t)\right]$, and ${\nabla }_{\Vert }^2={\partial }_x^2+{\partial }_y^2$, where ${\mathbf{r}}_{\Vert }=x\widehat{\mathbf{x}}+y\widehat{\mathbf{y}}$ and ${\mathbf{r}}_{\perp }=z\widehat{\mathbf{z}}$. ${\mathcal{D}}_{b,\Vert }$ and ${\mathcal{D}}_{b,\perp }$ are the diffusion coefficients in the parallel and perpendicular directions, and ${\rho }_{\mathrm{bulk}}(\mathbf{r},t)$ represents the distribution of particles in bulk, in units of particles per volume. Note that the diffusion coefficients ${\mathcal{D}}_{b,\left|\right|}$ and ${\mathcal{D}}_{b,\perp }$ are connected to diffusion in the plane $(x,y)$ and in the perpendicular direction z, respectively. We should mention that the form of the diffusion coefficient poses a scale-dependent dispersivity for diffusion in the $\widehat{z}$ direction and thus creates heterogeneity in the bulk. These forms of the diffusion coefficients allow us to consider anisotropic diffusion and enable us to analyze scenarios related to anomalous diffusion, where the mean square displacement has a nonlinear time dependence. We underline that similar spatial dependence on the diffusion coefficient has successfully been used to investigate diffusion on fractals [23,24,25], turbulence [26], solute transport in fractal porous media [27,28], and atom deposition in a porous substrate [29]. For the isotropic situations, we have ${\mathcal{D}}_{b,\left|\right|}={\mathcal{D}}_{b,\perp }$ in the bulk.

For the processes that occur on the surface, we assume that the bulk particles can be adsorbed by the surface located in $z=0$. Once adsorbed, particles can diffuse within the surface and a reaction process may take place; this is intimately related to several technologies, such as catalysts for the production of biochemical sensors [30], fuel [31], energy devices [32], and many others. Hence, the following equation is considered:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial }{\partial t}{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)& ={\mathcal{D}}_s{\nabla }_{\Vert }^2{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)+{\int }_0^{t}d{t}^{\prime }{k}_{\mathrm{ads}}(t-t^{\prime }){\rho }_{\mathrm{bulk}}(\mathbf{r},t^{\prime }){|}_{z=0}\hfill \\ & -{\int }_0^{t}d{t}^{\prime }{k}_{\mathrm{total}}(t-t^{\prime }){\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t),\hfill \end{array}\end{array}$$

(3)

where ${\mathcal{D}}_s$ is the diffusion coefficient for the particles’ concentration (or distribution) on the surface, ${\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)$ represents the density of particles on the surface (particles per area), $k_{\mathrm{total}}\left(t\right)$ is a kernel that can be connected to the desorption and reaction processes, i.e., $k_{\mathrm{total}}\left(t\right)=k_{\mathrm{surf},r}\left(t\right)+k_{\mathrm{desor}}\left(t\right)$, which may be present on the surface. $k_{\mathrm{desor}}\left(t\right)$ is related to the desorption rate, and $k_{\mathrm{surf},r}\left(t\right)$ is related to the reaction processes on the surfaces during the diffusion process on the surface. Equation (3) can be obtained by considering the continuity with additional terms related to the interaction between the surface and bulk, i.e.,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial }{\partial t}{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)& +{\nabla }_{\Vert }·{\mathcal{J}}_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)={\int }_0^{t}d{t}^{\prime }{k}_{\mathrm{ads}}(t-t^{\prime }){\rho }_{\mathrm{bulk}}(\mathbf{r},t^{\prime }){|}_{z=0}\hfill \\ & -{\int }_0^{t}d{t}^{\prime }{k}_{\mathrm{total}}(t-t^{\prime }){\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t^{\prime }),\hfill \end{array}\end{array}$$

(4)

combined with Fick’s law, i.e., ${\mathcal{J}}_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)=-{\mathcal{D}}_s{\nabla }_{\Vert }{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)$. In this manner, Equation (3) (or Equation (4)) has terms that represent an interaction between the surface and the particles after the adsorption processes, where they can be desorbed or promote the formation of other particles. In the latter case, the kernel may also be connected to the presence of intermediate processes during the reaction. The kernel $k_{\mathrm{ads}}\left(t\right)$ governs the adsorption process, i.e., it tells us what the interaction is between the surface particles and the bulk ones. Thus, depending on the nature of this interaction, the adsorption–desorption phenomena may or may not follow nonexponential decay behavior. For example, when we consider pure physisorption processes, the subsequent state of a particle depends on its preceding state, which is a way to characterize a non-Markovian process, i.e., a process involving a memory effect [33]. For $k_{\mathrm{ads}}\left(t\right)\propto \delta \left(t\right)$, the kernel is a localized function of time, implying that the preceding state does not matter to the actual state. These important features show that the kernel related to adsorption–desorption may account for a large class of effects that can be short- or long-ranged, depending on the processes occurring on the surface.

Equations (1) and (3) are coupled by the following condition

$$\begin{array}{c}\hfill {\left.\mathbf{n}·\left[{\mathcal{D}}_{b,\perp }{r}_{\perp }^{-\eta }{\nabla }_{\perp }{\rho }_{\mathrm{bulk}}(\mathbf{r},t)\right]\right|}_{z=0}=\frac{d}{dt}{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)+{\int }_0^{t}d{t}^{\prime }{k}_r(t-t^{\prime }){\left.{\rho }_{\mathrm{bulk}}(\mathbf{r},t)\right|}_{z=0}.\end{array}$$

(5)

In Equation (5), the adsorption–desorption processes are represented by the first term on the right-hand side; it couples the processes present on the surface with the ones in the bulk. The other term, i.e., the second term, represents a reaction process (see, for example, Refs. [34,35]), where the particles are removed from the bulk to the surface. The additional boundary conditions to analyze the problem defined above are

$$\begin{array}{ccc}& & {\left.{\partial }_z{\rho }_{\mathrm{bulk}}(\mathbf{r},t)\right|}_{z\to \infty }=0,{\left.{\partial }_x{\rho }_{\mathrm{bulk}}(\mathbf{r},t)\right|}_{x\to \pm \infty }=0,{\left.{\partial }_y{\rho }_{\mathrm{bulk}}(\mathbf{r},t)\right|}_{y\to \pm \infty }=0,\hfill \\ & & {\left.{\partial }_x{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\left|\right|},t)\right|}_{x\to \pm \infty }=0,\mathrm{and}{\left.{\partial }_y{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\left|\right|},t)\right|}_{y\to \pm \infty }=0.\hfill \end{array}$$

For the initial condition, we consider ${\rho }_{\mathrm{bulk}}(\mathbf{r},0)={\phi }_{\mathrm{bulk}}\left(\mathbf{r}\right)$ (so any initial configuration is allowed in the bulk) and, for simplicity, ${\rho }_{\mathrm{surf}}({\mathbf{r}}_{\left|\right|},0)=0$, which implies that, initially, there is no concentration of particles on the surfaces.

Performing some calculations, we can show that the processes on the surface modify the bulk, i.e.,

$$\begin{array}{c}\hfill \frac{d}{dt}\left[\int d{\mathbf{r}}_{\perp }\int d{\mathbf{r}}_{\Vert }{\rho }_{\mathrm{bulk}}(\mathbf{r},t)+\int d{\mathbf{r}}_{\Vert }{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)\right]=-{\int }_0^{t}d{t}^{\prime }{k}_r(t-t^{\prime }){\left.{\rho }_{\mathrm{bulk}}(\mathbf{r},t)\right|}_{z=0},\end{array}$$

(6)

where $\int d{\mathbf{r}}_{\perp }\equiv {\int }_0^{\infty }dz$, $\int d{\mathbf{r}}_{\Vert }\equiv {\int }_{-\infty }^{\infty }dx{\int }_{-\infty }^{\infty }dy$, and, consequently, for $k_r\left(t\right)=0$, we have

$$\begin{array}{c}\hfill \int d{\mathbf{r}}_{\perp }\int d{\mathbf{r}}_{\Vert }{\rho }_{\mathrm{bulk}}(\mathbf{r},t)+\int d{\mathbf{r}}_{\Vert }{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)=\mathrm{constant},\end{array}$$

(7)

which is a direct consequence of conserving the total number of particles in the system.

We now investigate the spatial and time behavior of the density of particles on the surface and in the bulk from analytical and numerical points of view. Thus, we start by obtaining the analytical solution for this problem, that is, closed expressions for ${\rho }_{\mathrm{bulk}}(\mathbf{r},t)$ and ${\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)$. To do this, we use integral transforms and the Green function approach [36,37]. One of the integral transforms is the Laplace transform, defined as

$$\mathcal{L}\{{\rho }_{\mathrm{bulk}}(\mathbf{r},t);s\}={\int }_0^{\infty }dt{e}^{-st}{\rho }_{\mathrm{bulk}}(\mathbf{r},t)={\widehat{\rho }}_{\mathrm{bulk}}(\mathbf{r},s),$$

(8)

and its inverse,

$${\mathcal{L}}^{-1}\{{\widehat{\rho }}_{\mathrm{bulk}}(\mathbf{r},s);t\}=\frac{1}{2\pi i}{\int }_{-i\infty +c}^{i\infty +c}ds{e}^{st}{\widehat{\rho }}_{\mathrm{bulk}}(\mathbf{r},s)={\rho }_{\mathrm{bulk}}(\mathbf{r},t).$$

(9)

In addition, we use the special integral transform that can be constructed from the eigenfunctions of the Sturm–Liouville problem related to the following differential equation:

$$\frac{\partial }{\partial z}\left\{z^{-\eta }\frac{\partial }{\partial z}\psi \left(z,k_z\right)\right\}=-k_z^{2+\eta }\psi \left(z,k_z\right),$$

(10)

subjected to the boundary condition $\left|\psi \right(\infty ,k_z\left)\right|<\infty$. The eigenfunctions obtained from Equation (10) are given by

$$\psi (z,k_z)={\left(z{k}_z\right)}^{\frac{1}{2}(1+\eta )}{J}_{-\nu }\left[\frac{2}{2+\eta }{\left(k_{z}z\right)}^{\frac{1}{2}(2+\eta )}\right]$$

(11)

where $J_{\nu }\left(x\right)$ denotes the Bessel function [37] with order $\nu =(1+\eta )/(2+\eta )$. Equation (11) allows us to define the following integral transform:

$$\begin{array}{ccc}\hfill {\mathcal{F}}_{\eta }\left\{{\rho }_{\mathrm{bulk}}(\mathbf{r},t);z\right\}& =& {\int }_0^{\infty }dz\psi (z,k_z){\rho }_{\mathrm{bulk}}(\mathbf{r},t)={\tilde{\rho }}_{\mathrm{bulk}}({\mathbf{r}}_{\Vert },k_z,t),\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\mathcal{F}}_{\eta }^{-1}\left\{{\tilde{\rho }}_{\mathrm{bulk}}({\mathbf{r}}_{\Vert },k_z,t);k_z\right\}& =& {\int }_0^{\infty }d{k}_z\psi (z,k_z){\tilde{\rho }}_{\mathrm{bulk}}({\mathbf{r}}_{\Vert },k_z,t)={\rho }_{\mathrm{bulk}}(\mathbf{r},t).\hfill \end{array}$$

(13)

We observe that Equations (12) and (13) may be related to a generalized Hankel transform [38,39,40,41]. By using these integral transforms and Fourier transform

$${\mathcal{F}}_{x,y}\{{\rho }_{\mathrm{bulk}}(\mathbf{r},t);x,y\}={\int }_{-\infty }^{\infty }dx{e}^{-i{k}_{x}x}{\int }_{-\infty }^{\infty }dy{e}^{-i{k}_{y}y}{\rho }_{\mathrm{bulk}}(\mathbf{r},t)={\overline{\rho }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z,t)$$

(14)

and it inverse

$${\mathcal{F}}_{x,y}^{-1}\{{\overline{\rho }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z,t);k_x,k_y\}=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }d{k}_x{e}^{i{k}_{x}x}\frac{1}{2\pi }{\int }_{-\infty }^{\infty }d{k}_y{e}^{i{k}_{y}y}{\overline{\rho }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z,t)={\rho }_{\mathrm{bulk}}(\mathbf{r},t),$$

(15)

where ${\mathbf{k}}_{\Vert }=(k_x,k_y)$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\overline{\rho }}}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z,s)& =-{\int }_0^{\infty }d{z}^{\prime }{\tilde{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },z,z^{\prime },s)\hfill \\ & +\left[s{\widehat{\overline{\rho }}}_{\mathrm{surf}}({\mathbf{k}}_{\Vert },s)-{\widehat{k}}_r\left(s\right){\widehat{\overline{\rho }}}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },0,s)\right]\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z,s),\hfill \end{array}\end{array}$$

(16)

and

$$\begin{array}{c}\hfill {\widehat{\overline{\rho }}}_{\mathrm{surf}}({\mathbf{k}}_{\Vert },s)=\frac{{\widehat{k}}_{\mathrm{ads}}\left(s\right)}{s+{\mathcal{D}}_s{\left|{\mathbf{k}}_{\Vert }\right|}^2+{\widehat{k}}_{\mathrm{total}}\left(s\right)}{\widehat{\overline{\rho }}}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },0,s),\end{array}$$

(17)

for the initial condition previously defined, which assumes that the particles are initially in bulk. The Green function is given by

$$\begin{array}{c}\hfill \overline{\mathcal{G}}({\mathbf{k}}_{\Vert },z,z^{\prime },t)=-\frac{e^{-{\mathcal{D}}_{b,\Vert }{\left|{\mathbf{k}}_{\Vert }\right|}^{2}t}}{(2+\eta ){\mathcal{D}}_{b,\perp }t}{\left(z{z}^{\prime }\right)}^{\frac{1}{2}\left(1+\eta \right)}{e}^{-\frac{1}{{(2+\eta )}^2{\mathcal{D}}_{b,\perp }t}\left(z^{2+\eta }+z^{\prime 2+\eta }\right)}{\mathrm{I}}_{-\nu }\left[\frac{2{\left(z{z}^{\prime }\right)}^{\frac{1}{2}(2+\eta )}}{{(2+\eta )}^2{\mathcal{D}}_{b,\perp }t}\right],\end{array}$$

(18)

where $I_{-\nu }\left(x\right)$ is the Bessel function of the modified argument [37]. From the inverse of the Fourier transform, we obtain the Green function in the form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{G}({\mathbf{r}}_{\Vert },z,z^{\prime },t)& =-\frac{1}{4\pi {\mathcal{D}}_{b,\Vert }t}{e}^{\frac{|{\mathbf{r}}_{\Vert }{|}^{\mathbf{2}}}{4{\mathcal{D}}_{b,\Vert }t}}\hfill \\ & \times \frac{{\left(z{z}^{\prime }\right)}^{\frac{1}{2}\left(1+\eta \right)}}{(2+\eta ){\mathcal{D}}_{b,\perp }t}{e}^{-\frac{1}{{(2+\eta )}^2{\mathcal{D}}_{b,\perp }t}\left(z^{2+\eta }+z^{\prime 2+\eta }\right)}{\mathrm{I}}_{-\nu }\left[\frac{2{\left(z{z}^{\prime }\right)}^{\frac{1}{2}(2+\eta )}}{{(2+\eta )}^2{\mathcal{D}}_{b,\perp }t}\right].\hfill \end{array}\end{array}$$

(19)

It is illustrative to look at the behavior of the Green function (Equation (19)) exhibited in Figure 1a,b. These are depicted for $\eta =-0.5$, while Figure 2a,b illustrate the behavior of the Green function for $\eta =0.5$. The spatial dependence on the diffusion coefficient obtained for different values of $\eta$ is responsible for different behaviors of the Green function and, consequently, for the spread of the distributions when the system is characterized by this type of heterogeneity.

Before proceeding, we underline that Equation (18) is obtained by solving the following equation

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{D}}_{b,\perp }\frac{\partial }{\partial z}\left\{z^{-\eta }\frac{\partial }{\partial z}\overline{\mathcal{G}}({\mathbf{k}}_{\Vert },z,z^{\prime },t)\right\}& -{\mathcal{D}}_{b,\Vert }{\left|{\mathbf{k}}_{\Vert }\right|}^2\overline{\mathcal{G}}({\mathbf{k}}_{\Vert },z,z^{\prime },t)\hfill \\ \hfill & -\frac{\partial }{\partial t}\overline{\mathcal{G}}({\mathbf{k}}_{\Vert },z,z^{\prime },t)=\delta (z-z^{\prime })\delta \left(t\right),\hfill \end{array}\end{array}$$

(20)

subjected to the boundary conditions

$$\begin{array}{c}\hfill {\left.{\mathcal{D}}_{b,\perp }{z}^{-\eta }\frac{\partial }{\partial z}\overline{\mathcal{G}}({\mathbf{k}}_{\Vert },z,z^{\prime },t)\right|}_{z=0}=0,{\left.{\mathcal{D}}_{b,\perp }{z}^{-\eta }\frac{\partial }{\partial z}\overline{\mathcal{G}}({\mathbf{k}}_{\Vert },z,z^{\prime },t)\right|}_{z=\infty }=0,\end{array}$$

(21)

taking into account that $\overline{\mathcal{G}}({\mathbf{k}}_{\Vert },z,z^{\prime },t)=0$ for $t<0$.

Using these results, it is possible to obtain the profile of the density of particles in the bulk and on the surface. In particular, before the inversion procedures, we may notice that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\overline{\rho }}}_{\mathrm{surf}}({\mathbf{k}}_{\Vert },s)& =-\frac{{\widehat{k}}_{\mathrm{ads}}(s)}{(s+{\mathcal{D}}_s|{\mathbf{k}}_{\Vert }{|}^2+{\widehat{k}}_{\mathrm{total}}(s))(1-k_r(s)\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,0,s))-s{\widehat{k}}_{\mathrm{ads}}(s)\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,0,s)}\hfill \\ & \times {\int }_0^{\infty }d{z}^{\prime }{\overline{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z^{\prime },s),\hfill \end{array}\end{array}$$

(22)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\overline{\rho }}}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z,s)& =-{\int }_0^{\infty }d{z}^{\prime }{\tilde{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },z,z^{\prime },s)\hfill \\ & -\frac{(s+{\mathcal{D}}_s|{\mathbf{k}}_{\Vert }{|}^2+{\widehat{k}}_{\mathrm{total}}(s))\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z,s)}{(s+{\mathcal{D}}_s|{\mathbf{k}}_{\Vert }{|}^2+{\widehat{k}}_{\mathrm{total}}(s))(1-k_r(s)\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,0,s))-s{\widehat{k}}_{\mathrm{ads}}(s)\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,0,s)}\hfill \\ & \times \left[\frac{s{k}_{\mathrm{ads}}(s)}{s+{\mathcal{D}}_s{|{\mathbf{k}}_{\Vert }|}^2+{\widehat{k}}_{\mathrm{total}}(s)}+k_r(s)\right]{\int }_0^{\infty }d{z}^{\prime }{\overline{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z^{\prime },s).\hfill \end{array}\end{array}$$

(23)

At this point, we have the tools in the Laplace–Fourier space that are needed to consider the following two particular cases of the previous equations.

The first refers to the absorption process of particles by the surface in the absence of reaction processes, that is, when ${\widehat{k}}_{\mathrm{total}}\left(s\right)=0$ and ${\widehat{k}}_r\left(s\right)=0$, at a constant absorption rate, that is, ${\widehat{k}}_{\mathrm{ads}}\left(s\right)=k_{\mathrm{ads}}$, with the particles initially in bulk. In this case, Equations (22) and (23) are reduced, respectively, to

$$\begin{array}{c}\hfill {\widehat{\overline{\rho }}}_{\mathrm{surf}}({\mathbf{k}}_{\Vert },s)=-\frac{k_{\mathrm{ads}}}{s+{\mathcal{D}}_s{\left|{\mathbf{k}}_{\Vert }\right|}^2-s{k}_{\mathrm{ads}}\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,0,s)}{\int }_0^{\infty }d{z}^{\prime }{\overline{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z^{\prime },s)\end{array}$$

(24)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\overline{\rho }}}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z,s)& =-{\int }_0^{\infty }d{z}^{\prime }{\tilde{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },z,z^{\prime },s)\hfill \\ & -\frac{s{k}_{\mathrm{ads}}\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z,s)}{s+{\mathcal{D}}_s{\left|{\mathbf{k}}_{\Vert }\right|}^2-s{k}_{\mathrm{ads}}\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,0,s)}{\int }_0^{\infty }d{z}^{\prime }{\overline{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z^{\prime },s).\hfill \end{array}\end{array}$$

(25)

After obtaining the inverses of the Laplace and Fourier transforms of both equations, we arrive at the final expressions for the density of particles on the surface and in the bulk, namely,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },s)& =-k_{\mathrm{ads}}\int d{\mathbf{r}}_{\Vert }^{\prime }{\int }_0^{t}d{t}^{\prime }{\mathfrak{G}}^{\left(1\right)}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },t-t^{\prime })\hfill \\ & \times \int d{\mathbf{r}}_{\Vert }^{\prime \prime }{\int }_0^{\infty }d{z}^{\prime }{\phi }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert }^{\prime }-{\mathbf{r}}_{\Vert }^{\prime \prime },z^{\prime })\mathcal{G}({\mathbf{r}}_{\Vert },0,z^{\prime },t^{\prime })\hfill \end{array}\end{array}$$

(26)

with

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{G}}^{\left(1\right)}({\mathbf{r}}_{\Vert },t)& =\sum _{m=0}^n\left(\begin{array}{c}n\\ m\end{array}\right){\int }_0^{t}d{t}^{\prime }\frac{{(t-t^{\prime })}^{\frac{1+\eta }{2+\eta }n-1}{t}^{\prime m}}{\mathrm{\Gamma }\left(\frac{1+\eta }{2+\eta }n\right)\mathrm{\Gamma }(1+m)}\hfill \\ & \times \int {\mathbf{r}}_{\Vert }^{\prime }{\mathcal{G}}_{{\mathcal{D}}_s,xy}({\mathbf{r}}_{\Vert }^{\prime },t-t^{\prime })\frac{{\partial }^m}{\partial t^{\prime m}}{\mathcal{G}}_{{\mathcal{D}}_{b,\Vert },xy}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },t^{\prime }),\hfill \end{array}\end{array}$$

(27)

where

$$\begin{array}{c}\hfill {\mathcal{G}}_{{\mathcal{D}}_i,xy}({\mathbf{r}}_{\Vert },t)=\frac{1}{4\pi {\mathcal{D}}_{i}t}exp\left(-\frac{{\mathbf{r}}_{\Vert }^2}{4{\mathcal{D}}_{i}t}\right),\end{array}$$

(28)

and

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert },z,s)& =& -\int d{\mathbf{r}}_{\Vert }^{\prime }{\int }_0^{\infty }d{z}^{\prime }{\phi }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert }^{\prime },z^{\prime })\mathcal{G}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },z,z^{\prime },t)\hfill \\ & -& k_{\mathrm{ads}}\int d{\mathbf{r}}_{\Vert }^{\prime }{\int }_0^{t}d{t}^{\prime }{\mathfrak{G}}^{\left(1\right)}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },t-t^{\prime }){\int }_0^{t^{\prime }}d{t}^{\prime \prime }\int d{\mathbf{r}}_{\Vert }^{\prime \prime }\mathcal{G}({\mathbf{r}}_{\Vert }^{\prime }-{\mathbf{r}}_{\Vert }^{\prime \prime },0,z,t^{\prime }-t^{\prime \prime })\hfill \\ & \times & \int d{\mathbf{r}}_{\Vert }^{\prime \prime \prime }{\int }_0^{\infty }d{z}^{\prime }{\phi }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert }^{\prime \prime }-{\mathbf{r}}_{\Vert }^{\prime \prime \prime },z^{\prime })\mathcal{G}({\mathbf{r}}_{\Vert }^{\prime \prime \prime },0,z^{\prime },t^{\prime \prime }).\hfill \end{array}$$

(29)

Figure 3 illustrates an important analytical result, the survival probability, which can be constructed from the solution that accounts for the particles absorbed by the surface. Indeed, from Equation (29), we obtain the survival probability, i.e.,

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{bulk}}\left(t\right)=\int d{\mathbf{r}}_{\Vert }{\int }_0^{\infty }dz{\rho }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert },z,t),\end{array}$$

(30)

and evaluate

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{surf}}\left(t\right)=\int d{\mathbf{r}}_{\Vert }{\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },t)=1-{\mathcal{S}}_{\mathrm{bulk}}\left(t\right).\end{array}$$

(31)

Both equations are related to the probability of finding particles in bulk (${\mathcal{S}}_{\mathrm{bulk}}\left(t\right)$) or on the surface (${\mathcal{S}}_{\mathrm{surf}}\left(t\right)$) and, consequently, correspond to the fractions of particles present in bulk and on the surface.

In Figure 3, we also show the result obtained by numerical simulation using the Langevin equation for the particles on the surface and in the bulk. They are defined as follows for the $x-y$ direction:

$$\begin{array}{c}\hfill \begin{array}{cc}{x}_{t+h}=x_t+\sqrt{Dh}{\zeta }_x\left(t\right),& -\infty <x<\infty \\ y_{t+h}=y_t+\sqrt{Dh}{\zeta }_y\left(t\right),& -\infty <y<\infty ,\end{array}\end{array}$$

(32)

and, consequently, we have information on ${\mathbf{r}}_{\Vert }=x\left(t\right)\widehat{\mathbf{x}}+y\left(t\right)\widehat{\mathbf{y}}$. Note that the Langevin equations for the x and y directions were used to simulate the particles on the surface and in bulk. For the z direction, the Langevin equation for the particles, in bulk, is given by

$$\begin{array}{c}\hfill z_{t+h}=z_t-\left[\left(\frac{\eta Dh}{2}sign\left(z_t\right)\left(\right|z_t{\left|\right)}^{-\eta -1}\right)-(\sqrt{Dh}\left(\right|z_t{\left|\right)}^{-\frac{\eta }{2}}){\zeta }_z\left(t\right)\right],0\le z<\infty .\end{array}$$

(33)

In these stochastic equations, ${\zeta }_i\left(t\right)$ ($i=x,y,$ and z) is white Gaussian noise with a normalized deviation generated using the Box–Muller method [42]. In addition, $\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 })$. The pseudorandom number generator utilized in this work was the maximally equidistributed combined Tausworthe generator [43] implemented via the tauss88 function from the C++ library Boost [44] (for more details, see Appendix A). For this case, we consider that the surface only absorbs particles governed by the Langevin equations, i.e., Equations (32) and (33). Figure 4 and Figure 5 illustrate the results obtained by using the previous Langevin equations for the directions x, y, and z when $\eta =-0.5$ and $\eta =0.5$. These figures correspond to the projections $(x,z)$ and $(x,y)$ to show how heterogeneity influences the dynamics aspects of the system. An animated version of the $(x,z)$ projections can be seen at https://youtu.be/lF6WpIQI-c4 (accessed on 6 March 2024).

The second case is characterized by a surface that can adsorb and desorb at some rates, i.e., ${\widehat{k}}_{\mathrm{total}}\left(s\right)=k_{\mathrm{desor}}$, ${\widehat{k}}_r\left(s\right)=0$, and ${\widehat{k}}_{\mathrm{ads}}\left(s\right)=k_{\mathrm{ads}}$. In this case, for the particles initially in bulk, from Equations (22) and (23), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\overline{\rho }}}_{\mathrm{surf}}({\mathbf{k}}_{\Vert },s)& =-\frac{k_{\mathrm{ads}}}{s+{\mathcal{D}}_s{\left|{\mathbf{k}}_{\Vert }\right|}^2+k_{desorb}-s{k}_{\mathrm{ads}}\widehat{\overline{\mathcal{G}}}(0,0,z^{\prime },s)}\hfill \\ & \times {\int }_0^{\infty }d{z}^{\prime }{\overline{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z^{\prime },s),\hfill \end{array}\end{array}$$

(34)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{\overline{\rho }}}_{\mathrm{bulk}}& ({\mathbf{k}}_{\Vert },z,s)=-{\int }_0^{\infty }d{z}^{\prime }{\tilde{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },z,z^{\prime },s)\hfill \\ & -\frac{s{k}_{\mathrm{ads}}\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z,s)}{s+{\mathcal{D}}_s{\left|{\mathbf{k}}_{\Vert }\right|}^2+k_{\mathrm{desor}}-s{k}_{\mathrm{ads}}\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,0,s)}{\int }_0^{\infty }d{z}^{\prime }{\overline{\phi }}_{\mathrm{bulk}}({\mathbf{k}}_{\Vert },z^{\prime })\widehat{\overline{\mathcal{G}}}({\mathbf{k}}_{\Vert },0,z^{\prime },s).\hfill \end{array}\end{array}$$

(35)

After obtaining the inverses of the Fourier and Laplace transforms, we obtain the following analytical results:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_{\mathrm{surf}}({\mathbf{r}}_{\Vert },s)& =-k_{\mathrm{ads}}\int d{\mathbf{r}}_{\Vert }^{\prime }{\int }_0^{t}d{t}^{\prime }{\mathfrak{G}}^{\left(2\right)}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },t-t^{\prime })\hfill \\ & \times \int d{\mathbf{r}}_{\Vert }^{\prime \prime }{\int }_0^{\infty }d{z}^{\prime }{\phi }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert }^{\prime }-{\mathbf{r}}_{\Vert }^{\prime \prime },z^{\prime })\mathcal{G}({\mathbf{r}}_{\Vert },0,z^{\prime },t^{\prime })\hfill \end{array}\end{array}$$

(36)

with

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{G}}^{\left(2\right)}({\mathbf{r}}_{\Vert },t)& =\sum _{m=0}^n\left(\begin{array}{c}n\\ m\end{array}\right){\int }_0^{t}d{t}^{\prime }\frac{{(t-t^{\prime })}^{\frac{1+\eta }{2+\eta }n-1}{t}^{\prime m}}{\Gamma \left(\frac{1+\eta }{2+\eta }n\right)\Gamma (1+m)}\hfill \\ & \times \int d{\mathbf{r}}_{\Vert }^{\prime }{e}^{-k_{\mathrm{desor}}(t-t^{\prime })}{\mathcal{G}}_{{\mathcal{D}}_s,xy}({\mathbf{r}}_{\Vert }^{\prime },t-t^{\prime })\frac{{\partial }^m}{\partial t^{\prime m}}\left[e^{-k_{\mathrm{desor}}{t}^{\prime }}{\mathcal{G}}_{{\mathcal{D}}_{b,\Vert },xy}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },t^{\prime })\right],\hfill \end{array}\end{array}$$

(37)

and

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert },z,s)& =& -\int d{\mathbf{r}}_{\Vert }^{\prime }{\int }_0^{\infty }d{z}^{\prime }{\phi }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert }^{\prime },z^{\prime })\mathcal{G}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },z,z^{\prime },t)\hfill \\ & -& k_{\mathrm{ads}}\int d{\mathbf{r}}_{\Vert }^{\prime }{\int }_0^{t}d{t}^{\prime }{\mathfrak{G}}^{\left(2\right)}({\mathbf{r}}_{\Vert }-{\mathbf{r}}_{\Vert }^{\prime },t-t^{\prime }){\int }_0^{t^{\prime }}d{t}^{\prime \prime }\int d{\mathbf{r}}_{\Vert }^{\prime \prime }\mathcal{G}({\mathbf{r}}_{\Vert }^{\prime }-{\mathbf{r}}_{\Vert }^{\prime \prime },0,z,t^{\prime }-t^{\prime \prime })\hfill \\ & \times & \int d{\mathbf{r}}_{\Vert }^{\prime \prime \prime }{\int }_0^{\infty }d{z}^{\prime }{\phi }_{\mathrm{bulk}}({\mathbf{r}}_{\Vert }^{\prime \prime }-{\mathbf{r}}_{\Vert }^{\prime \prime \prime },z^{\prime })\mathcal{G}({\mathbf{r}}_{\Vert }^{\prime \prime \prime },0,z^{\prime },t^{\prime \prime }).\hfill \end{array}$$

(38)

In Figure 6, the number of particles adsorbed, obtained from the analytical and numerical points of view, is represented. We underline that the numerical simulations were performed considering that the desorption process of the particles, from the surface to the bulk, follows the conditions

$$\begin{array}{c}\hfill \begin{array}{cc}{z}_{t+h}=0.5& —\mathrm{desorption}\mathrm{probability}(\mathrm{surface}\mathrm{to}\mathrm{the}\mathrm{bulk})—p_{des}=rh\\ z_{t+h}=z_t=0& —\mathrm{non-desorption}\mathrm{probability}—p_{stay}=1-rh\end{array}\end{array}$$

(39)

(for more details, see Appendix A).

3. Discussion and Conclusions

This work investigates a diffusion process in a heterogeneous medium with adsorption–desorption occurring on a surface. The particles adsorbed can diffuse on the surface and may desorb back into the bulk after some time. The bulk process is governed by a diffusion equation with a spatially dependent diffusion coefficient. The surface adsorption–desorption process is described by the kinetic functions $k_{\mathrm{ads}}\left(t\right)$ and $k_{\mathrm{desor}}\left(t\right)$, which represent the kernels in Equation (3).

We performed numerical simulations using the Langevin equation with multiplicative noise to gain additional insight into the diffusion process from a different perspective. This approach complements the analytical model. Figure 3 and Figure 6 compare the analytical and numerical results, showing good agreement for the cases analyzed. The first case was considered to have a purely adsorbing surface, which means that desorption was absent and all particles were adsorbed. The second case incorporated desorption, allowing particles to return to the bulk after some time. For this case, we adjusted the parameters related to adsorption–desorption in the analytical model to match the numerical results.

These findings extend previous work in [45,46] by incorporating bulk inhomogeneity and surface diffusion. We believe the results presented here can contribute to the study of diffusion processes with surface adsorption–desorption and the possibility of surface diffusion.