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

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.


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., ⟨(r − ⟨r⟩) 2 ⟩ ∼ 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., ⟨(r − ⟨r⟩) 2 ⟩ ∼ t α , where α < 1 or α > 1 corresponds to subor superdiffusion, respectively [5,6].In conjunction with these situations, adsorptiondesorption 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 solidliquid 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.

The Problem: Diffusion and Kinetics
Let us initiate our analysis by considering that the subsequent equation controls the diffusion process in the bulk: Equation ( 1) is the typical continuity equation applied to Fick's law, i.e., , 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 ẑ direction, but particles can diffuse in 3D.More specifically, D b,⊥ (r) = D b,⊥ r , where r ∥ = x x + y ŷ and r ⊥ = zẑ.D b,∥ and D b,⊥ are the diffusion coefficients in the parallel and perpendicular directions, and ρ bulk (r, t) represents the distribution of particles in bulk, in units of particles per volume.Note that the diffusion coefficients D b,|| and D b,⊥ 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 ẑ 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 D b,|| = D b,⊥ 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: where D s is the diffusion coefficient for the particles' concentration (or distribution) on the surface, ρ surf (r ∥ , t) represents the density of particles on the surface (particles per area), k total (t) is a kernel that can be connected to the desorption and reaction processes, i.e., k total (t) = k surf,r (t) + k desor (t), which may be present on the surface.k desor (t) is related to the desorption rate, and k surf,r (t) 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., combined with Fick's law, i.e., J surf (r (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 ads (t) 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 ads (t) ∝ δ(t), 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 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 For the initial condition, we consider ρ bulk (r, 0) = φ bulk (r) (so any initial configuration is allowed in the bulk) and, for simplicity, ρ surf (r || , 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., where dy, and, consequently, for k r (t) = 0, we have 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 ρ bulk (r, t) and ρ surf (r ∥ , 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 and its inverse, 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: subjected to the boundary condition |ψ(∞, k z )| < ∞.The eigenfunctions obtained from Equation (10) are given by where J ν (x) denotes the Bessel function [37] with order ν = (1 + η)/(2 + η).Equation (11) allows us to define the following integral transform: 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 and it inverse where and for the initial condition previously defined, which assumes that the particles are initially in bulk.The Green function is given by where I −ν (x) is the Bessel function of the modified argument [37].From the inverse of the Fourier transform, we obtain the Green function in the form It is illustrative to look at the behavior of the Green function (Equation ( 19)) exhibited in Figure 1a,b.These are depicted for η = −0.5, while Figure 2a,b illustrate the behavior of the Green function for η = 0.5.The spatial dependence on the diffusion coefficient obtained for different values of η 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 subjected to the boundary conditions taking into account that Ḡ(k ∥ , z, z ′ , 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 and 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 k total (s) = 0 and k r (s) = 0, at a constant absorption rate, that is, k ads (s) = k ads , with the particles initially in bulk.In this case, Equations ( 22) and ( 23) are reduced, respectively, to and 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 the bulk, namely, with where and 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., and evaluate Both equations are related to the probability of finding particles in bulk (S bulk (t)) or on the surface (S surf (t)) and, consequently, correspond to the fractions of particles present in bulk and on the surface.As described in this section, the black solid lines correspond to the numerical simulations using the Langevin equation, and the red dashed-dotted lines correspond to the analytical model.For the analytic model, we consider, with illustrative purposes, k desor = 0.0, k ads → ∞ (total absorption process), D b,∥ = D b,⊥ = 2.0, and z ′ = 4.0, in arbitrary units.In the numerical simulations, we considered 200 k particles with the same starting conditions as in the analytical model; the desorption process is absent (rh = 0%/step), D = 2.0, and h = 0.01.
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: and, consequently, we have information on r ∥ = x(t) x + y(t) ŷ.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 In these stochastic equations, ζ i (t) (i = x, y, and z) is white Gaussian noise with a normalized deviation generated using the Box-Muller method [42].In addition, ⟨ζ i (t 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).Figures 4 and 5 illustrate the results obtained by using the previous Langevin equations for the directions x, y, and z when η = −0.5 and η = 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).Note that the spatial dependence of the diffusion coefficient, as in Figure 1a,b directly influences the system's behavior.This was a short run with 250 steps and 50M particles with the initial position on (x 0 = 0, y 0 = 0, z 0 = 4.0).In this simulation, h = 0.01, D = 2.0, and the desorption process is absent.The second case is characterized by a surface that can adsorb and desorb at some rates, i.e., k total (s) = k desor , k r (s) = 0, and k ads (s) = k ads .In this case, for the particles initially in bulk, from Equations ( 22) and ( 23), we obtain and After obtaining the inverses of the Fourier and Laplace transforms, we obtain the following analytical results: with and 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 z t+h = 0.5 -desorption probability (surface to the bulk) -p des = rh z t+h = z t = 0 -non-desorption probability -p stay = 1 − rh (for more details, see Appendix A).The black solid lines correspond to the numerical simulations, as described in the text, whereas the red dashed-dotted lines correspond to the analytical model.We have used the model to adjust the data obtained from the numerical simulations with the Langevin equations.Thus, the kinetic parameters, i.e., k ads and k desor , adjusted both the model and the numerical simulations.The parameter values for η = 0.5 were k desor = 1.5 and k ads = 0.9; for η = 0.0 were k desor = 0.45 and k ads = 2.9; and for η = −0.5 were k desor = 0.2 and k ads = 4.5.For simplicity, we consider D b,∥ = D b,⊥ = 2.0 and z ′ = 4.0, in arbitrary units.These numerical simulations utilized the parameters shown in Figure 3, except for rh, which is equal to 5%/step.

Discussion and Conclusions
This work investigates a diffusion process in a heterogeneous medium with adsorptiondesorption 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 adsorptiondesorption process is described by the kinetic functions k ads (t) and k desor (t), 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.Figures 3 and 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.If the particle is not on the surface, it is in the bulk.The code will compute the diffusion in z and check if the particle was absorbed.

Figure 1 .
Figure 1.A three-dimensional density plot for the Green function (Equation (19)) when η = −0.5.In (a) a cut is represented through the center in the y direction, while (b) shows the behavior of the Green function around the point r = (0, 0, 1), i.e., how the particles are distributed.For simplicity, we consider D b,∥ t = D b,⊥ t = 0.5 and z ′ = 1.0, in arbitrary units.

Figure 2 .
Figure 2. A three-dimensional density plot for the Green function (Equation (19)) when η = 0.5.In (a) a cut is represented through the center in the y direction, while (b) shows the behavior of the Green function around the point r = (0, 0, 1), i.e., how the particles are distributed.For simplicity, we consider D b,∥ t = D b,⊥ t = 0.5 and z ′ = 1.0, in arbitrary units.

5 Figure 3 .
Figure 3. Trends of the survival probability for the adsorption-desorption case for different values of η.As described in this section, the black solid lines correspond to the numerical simulations using the Langevin equation, and the red dashed-dotted lines correspond to the analytical model.For the analytic model, we consider, with illustrative purposes, k desor = 0.0, k ads → ∞ (total absorption process), D b,∥ = D b,⊥ = 2.0, and z ′ = 4.0, in arbitrary units.In the numerical simulations, we considered 200 k particles with the same starting conditions as in the analytical model; the desorption process is absent (rh = 0%/step), D = 2.0, and h = 0.01.

Figure 4 .
Figure 4. Probability density maps for simulations with η = −0.5.(a,b) represent the (x, z) and (x, y) of the bulk and surface densities at a time equal to 1.0.(c,d) correspond to (a,b) at the time 2.5.Note that the spatial dependence of the diffusion coefficient, as in Figure1a,b directly influences the system's behavior.This was a short run with 250 steps and 50M particles with the initial position on (x 0 = 0, y 0 = 0, z 0 = 4.0).In this simulation, h = 0.01, D = 2.0, and the desorption process is absent.

Figure 5 .
Figure 5. Probability density maps for simulations with η = 0.5.(a,b) represent the (x, z) and (x, y) projections of the bulk and surface densities at a time equal to 1.0.(c,d)advance the time to 2.5.Apart from the value of η, all coefficients were the same as in Figure4and the desorption process is absent.

5 Figure 6 .
Figure 6.Trends of the survival probability for the adsorption-desorption case with different values of η.The black solid lines correspond to the numerical simulations, as described in the text, whereas the red dashed-dotted lines correspond to the analytical model.We have used the model to adjust the data obtained from the numerical simulations with the Langevin equations.Thus, the kinetic parameters, i.e., k ads and k desor , adjusted both the model and the numerical simulations.The parameter values for η = 0.5 were k desor = 1.5 and k ads = 0.9; for η = 0.0 were k desor = 0.45 and k ads = 2.9; and for η = −0.5 were k desor = 0.2 and k ads = 4.5.For simplicity, we consider D b,∥ = D b,⊥ = 2.0 and z ′ = 4.0, in arbitrary units.These numerical simulations utilized the parameters shown in Figure3, except for rh, which is equal to 5%/step.

Figure A2 .
Figure A2.This part of the code applies Equation (A2) to model stochastic desorption.RNG means random number generator.

Figure A3 .
Figure A3.Here, we calculated diffusion in z and checked for adsorption.Then, "nsurf" can be used to calculate the ratio of particles on the surface.Finally, we can generate the Gaussian distribution.