Ornstein–Uhlenbeck Process on Three-Dimensional Comb under Stochastic Resetting

: The Ornstein–Uhlenbeck (O-U) process with resetting is considered as the anomalous transport taking place on a three-dimensional comb. The three-dimensional comb is a comb inside a comb structure, consisting of backbones and ﬁngers in the following geometrical correspondence x –backbone → y –ﬁngers–backbone → z –ﬁngers. Realisation of the O-U process on the three-dimensional comb leads to anomalous (non-Markovian) diffusion. This speciﬁc anomalous transport in the presence of resets results in non-equilibrium stationary states. Explicit analytical expressions for the mean values and the mean squared displacements along all three directions of the comb are obtained and veriﬁed numerically. The marginal probability density functions for each direction are obtained numerically by Monte Carlo simulation of a random transport described by a system of coupled Langevin equations for the comb geometry.


Introduction
Nowadays, processes employing anomalous diffusion emerging from different types of approaches have been greatly studied. These variations are frequently caused by unique constraints placed on the diffusing particles, such as the properties of the substrate or the nature of interactions between them. Examples of anomalous diffusion include the diffusion on fractals [1], diffusion with random stopping [2], which temporarily stuns particles during their diffusive movements, diffusion in disordered materials [3], protein mobility [4], dynamics of ultracold atoms [5], telomere diffusion in the cell nucleus [6], and diffusion with geometric constraints (comb-like models) [7][8][9][10][11][12], to name but a few, which are characterised by a power-law dependence of the mean squared displacement (MSD) on time, x 2 (t) ∼ t α , where α is anomalous diffusion exponent. For 0 < α < 1, the corresponding process is subdiffusive, for α > 1 -superdiffusive, while for α = 1, it is a normal diffusion. Diffusion on comb-like structures is a convenient way to study anomalous diffusion, and it has been extensively explored in the scientific literature, both the two-dimensional comb [13][14][15][16][17][18] and the three-dimensional comb [19,20].
These models have been demonstrated to be helpful in describing anomalous transport through porous solid pellets with different porous geometries [21]. Comb models are used to model electron transport in disordered nanostructured semiconductors [17,22], Three-dimensional comb structure with illustrated diffusive movement (orange balls and arrows, simulating movement of the particle) in its structural features, comprised of the backbone (x-axis), main fingers (y-axis), and secondary fingers (z-axis), and depicting the respective processes: Brownian motion (BM) and Ornstein-Uhlenbeck (O-U) in these directions.

O-U Processes on Three-Dimensional Comb
Diffusion and random walks on three-dimensional comb structures have been investigated in the literature (see, for example, Refs. [14,15,19,20]). In this paper, we generalise the standard diffusion process on the three-dimensional comb by introducing the following Fokker-Planck equation

∂ ∂t
P(x, y, z, t) = δ(y)δ(z) L FP,x P(x, y, z, t) + δ(z) L FP,y P(x, y, z, t) + σ 2 z 2 ∂ 2 ∂z 2 P(x, y, z, t), with the initial condition at the backbone P(x, y, z, t = 0) = δ(x − x 0 )δ(y)δ(z), and zero boundary conditions at infinity. The solution of this three-dimensional equation is the probability density function (PDF) P(x, y, z, t), which gives the position (x, y, z) of the particle in the three-dimensional, geometrically constrained space at time t. The Dirac δ-functions δ(y)δ(z) and δ(z) mean that the movement in the x-direction is allowed only at y = z = 0 and in the y-direction only at z = 0, respectively, while in the z-direction, the particle performs normal diffusion. The equation consists of three terms represented by (σ 2 i /2) ∂ 2 ∂i 2 , i = {x, y, z} for the three diffusing directions, respectively, responsible for the fluctuating nature of the diffusion, as well as two deterministic terms represented by λ j ∂ ∂j j, j = {x, y}, in control of the mean-reversion in the xand y-directions. Consequently, the Fokker-Planck operators for the backbone and for the main fingers read L FP,y ≡ λ y ∂ ∂y y + σ 2 respectively, where [λ x δ(y)δ(z)] is the rate of mean-reversion along the backbone, µ is the long-term mean reversion value, towards which the process along the backbone is driven to, λ y δ(z) is the rate of mean-reversion along the fingers, with y = 0 as the long-term mean reversion value, (σ 2 x /2) δ(y)δ(z) is the diffusion coefficient along the backbone, (σ 2 y /2) δ(z) is the diffusion coefficient along the main fingers, and (σ 2 z /2) is the diffusion coefficient along the secondary fingers. We note that we have omitted the mean-reversion term λ z along the secondary fingers (z-direction) because of the complexity of the analytical calculations that arise if that term is present. We leave it as an interesting problem that can be investigated in the future.

Fokker-Planck Equations for the Marginal PDFs
In what follows, it is convenient to introduce the marginal PDFs, which are defined by The Laplace transform, We find the solution to Equation (7) using the approach presented in Appendix B. In the first step, we derive the differential equation for the marginal PDF for two directions P 1,2 (i, j, t) where i, j = {x, y, z}, i = j, and the outcome of the second step is the derivation of the differential equations for the marginal PDFs for all three directions, presented in Equations (4)- (6). Therefore, for the marginal PDF P 1,2 (x, y, t), defined by in Laplace space, see Equation (A34) with 1 η(s) = s, we have Multiplying Equation (10) by s 1/2 and applying the inverse Laplace transform, one obtains the following time-fractional diffusion equation where RL D µ t is the Riemann-Liouville fractional derivative (A3) of order µ = 1/2. Therefore, the fingers along the z-direction are the traps where the particle hinders before it comes back to the main fingers and to the backbone. Since the motion in the z-direction is Brownian, the returning time probability to the main fingers scales as t −3/2 , which is the waiting time PDF, as well. According to a continuous time random walk theory, such behaviour results in the fractional Fokker-Planck equation with the time-fractional derivative of order 1/2 [2]. The next step is to present the solution ofP 1,2 (x, y, s) as followŝ Following the procedure presented in Appendix B, we arrive at the differential equation in Laplace space for the diffusion on the backbone, Equation (A44), Multiplying Equation (13) by s 1/2 − λ 2 1/2 and then by s 1/2 = s × s −1/2 , we get where we define the modified mean-reverting rates the corresponding modified Fokker-Planck operators read By the inverse Laplace transform of Equation (14), we receive the fractional diffusion equation for diffusion on the backbone Here RL D µ t f (t) is the Riemann-Liouville fractional derivative (A3) of order µ = 1/2, and the second term on the right side of the equation contains the Prabhakar derivative in Riemann-Liouville form, see Equation (A14) with r = 0.
For diffusion along the main fingers, we get the Fokker-Planck equation for the marginal PDF p 2 (y, t); see Ref. [62], and for the transport along the secondary fingers, we have a standard diffusion equation for Brownian motion, i.e.,

First Moment and MSD
The mean value and the MSD along the backbone are calculated by multiplying Equation (16) by x and x 2 , respectively, and integrating with respect to x, from −∞ to +∞. The following equations are obtained, which, by using the integrals at the end of Appendix B, yield and Applying the Laplace transform to Equation (22), we obtain and then its Laplace inversion yields Here the Laplace transform formula for the three-parameter M-L function is used; see Equation (A13) in Appendix A. For the MSD in Equation (21), in Laplace space, we have and by the inverse Laplace transform, we obtain

Global Resetting
Next, we consider a general approach to the O-U process by resetting the threedimensional comb. To describe this phenomenon, we introduce the PDF P(x, y, z, t, |x r , y r , z r ), which describes two independent processes: O-U transport and resetting with the rate r to the fixed position R r = (x r , y r , z r ). The latter can be any point on the comb. Here, we consider a global resetting when the point to which we reset is the initial point R r = (x 0 , 0, 0). This process can be described by both a system of Langevin equations and a diffusion equation. The system of Langevin equations, considered in this section, stands for numerical (Monte Carlo) simulations. Therefore, for the analytical discussion of the problem, we first consider the following partial differential equation with initial condition P(x, y, z, t = 0|x 0 , 0, 0) = δ(x − x 0 )δ(y)δ(z) and zero boundary conditions at infinity. The only difference with the Fokker-Planck equation in the case of no resetting is that in Equation (27) we have two additional terms. The first one represents the decrease in the probability of the latest position of the particle before the resetting happens −rP(x, y, z, t), and the second additional term represents the increase in the probability for the particle of being found at the initial position, rδ(x − x 0 )δ(y)δ(z), where the particle is reset.

Fokker-Planck Equation for the Marginal PDF
The differential equation for the marginal PDF along the backbone can be found by using the same approach as in the case without resetting. Therefore, one finds from where the following equation forP 1,2 (x, y, s) is obtained The analysis of the functionP 1,2 (x, y, s) gives us the expression (A41), in the following form Therefore, we end up with the differential equation for the marginal PDF along the backbone in Laplace space, If we multiply both sides of Equation (31) with (s + r) 1/2 − λ 2 and then with (s + r) 1/2 ≡ (s + r)(s + r) −1/2 , we obtain where by its inverse Laplace space, we acquire the following temporal partial differential equation This is a fractional differential equation for the transport of the particles along the backbone from where we can see that tempered Riemann-Liouville derivative (A5) and tempered Prabhakar derivative (A6) occur in the equation. From here, one can show that the following equality holds truê which is actually the renewal equation for the marginal PDF; see also [63], From Equation (34), in the long time limit, one observes that the particle reaches a non-equilibrium stationary state (NESS), given by For the propagation of the particles along the main fingers in the presence of resetting, we get [62] ∂ ∂t where TRL D µ 0+ f (t) is the so-called tempered Riemann-Liouville fractional derivative (A5) of order µ = 1/2, with tempering parameter r. For the secondary fingers, we have describing a Brownian motion with resetting along the secondary fingers. From these equations for the marginal PDFs along the main fingers and secondary fingers, it can be found that the NESSs are reached in the long time limit; see Ref. [62].

First Moment and MSD
The mean value and the MSD are obtained from Equation (34) by multiplying it by x and x 2 , respectively, and integrating from −∞ to +∞ over x. Thus, we find for the mean value, and for the MSD, where x(s) c and x 2 (s) c are defined in Equations (23) and (25). Then, we have and From here, the mean and the MSD can be found by numerical inverse Laplace transform in MATHEMATICA [64].

Numerical Simulations: Coupled Langevin Equations
The diffusive motion with global resetting on the three-dimensional comb structure can be simulated by the following coupled Langevin equations [20,63] (in the case of no resetting, we refer to Refs. [65,66]) for the movement along the backbone, for the movement along the main fingers, and for the movement along the secondary fingers.
Here ξ i , i = {x, y, z}, is the white noise with zero mean, ξ i (τ∆t) = 0, and correlation function ξ i (τ∆t)ξ i (τ ∆t) = δ((τ − τ )∆t). The functions A(y) and B(z) are introduced to model Dirac δ functions, , j = {y, z} and ε = 10 −3 ; see Ref. [62]. The factor of 2 is introduced for the technical reason to compensate for the additional degree of freedom, z. It allows the particle to make the correct motion on the main fingers in the y-direction and still be able to enter the secondary fingers in the z-direction. This is illustrated in Figure 2, where the red dots and arrows represent successive movements of the particle inside the comb structure. For the simulations of the marginal PDF along the backbone, the diffusion coefficient along the backbone and the mean-reverting rate are renormalised by factor 1/[2 D y ]; see Refs. [20,63]. The simulated trajectories are presented in Figure 3.  [62]. The red arrows and points illustrate the consecutive movements of the particle in the comb structure.  (43)-(45), where the black dots stand for the resetting events and the shaded and white areas represent the time between two consecutive resetting events,

Resetting to the Backbone
As admitted above, the three-dimensional comb suggests different resetting algorithms due to the arbitrariness of the resetting point. Therefore, besides global resetting to the initial position, one can consider another possibility of the resets. In particular, here we consider resetting to the backbone when R r = (x, 0, 0). In this scenario, the O-U process takes place in the xand y-directions, while Brownian motion is in the z-direction, as depicted in Figure 1. This process is described by the following equation with the initial condition, P r (x, y, z, t = 0|x, 0, 0) = δ(x − x 0 )δ(y)δ(z), and zero boundary conditions at infinity, where the last term, rδ(y)δ(z) ∞ −∞ ∞ −∞ P r (x, y , z , t|x, 0, 0) dz dy , gives us the mechanism of resetting to the backbone.

Fokker-Planck Equation for the Marginal PDF
This equation can be solved by using the same approach as the one presented in Appendix B. We present the solution ofP(x, y, z, s) in the form where (s + r) = D zr 2 (x, y, s) →r(x, y, s) = s + r D z →ĝ(x, y, s) = 1 2 Then, in Laplace space, Equation (46) becomes where we useP 1,2 (x, y, s) = ∞ −∞P (x, y, z, s) dz. Next, from Equation (A41), we find that the expressions ofĥ(x, s) andq(x, s) are the same as in the case of global resetting, In Laplace space, for the marginal PDF in the xdirection, we obtain the equation If we multiply both sides of Equation (52), first with (s + r) 1/2 − λ 2 and then with (s + r) 1/2 = (s + r)(s + r) −1/2 , we have By inverse Laplace transform, we get the equation in temporal space in the following way where truncated power-law and truncated M-L memory kernels occur in the equation.

First Moment and MSD
Now the mean value and the MSD are defined by Equation (54), and the standard procedure described in the previous section results in for the mean value, and for the MSD, correspondingly. The analytical transformation of the mean value and MSD in time representation is a complex operation and is left for further deliberation and analysis, and instead, in this paper, we use the numerical inverse Laplace transform of Equations (55) and (56) in MATHEMATICA [64].

Numerical Simulations: Coupled Langevin Equations
In complete analogy with Section 3.3, diffusion with resetting to the backbone can be simulated by the following coupled Langevin equations for the movement along the backbone, for the movement along the main fingers, and for the movement along the secondary fingers. Here, ξ i , i = {x, y, z}, is the white noise with zero mean, ξ i (τ∆t) = 0, and correlation function ξ i (τ∆t)ξ i (τ ∆t) = δ((τ − τ )∆t), and A(y) and B(z) are the approximations of the Dirac δ functions as defined above in Section 3.3. The essential difference from the case of global resetting is in the movement along the main backbone, which reflects the difference between the Langevin Equations (43) and (57), where resets to the initial position x = x 0 are absent. The simulated trajectories are reflected in Figure 4.

Resetting to the Main Fingers
Finally, we studied the third mechanism of resetting the particle diffusing on comblike structures, which is the mechanism of resetting to the main fingers with coordinates (x, y, z 0 ), letting the process evolve in the xand y-directions without interruption. The partial differential equation for this process is given by the equation ∂ ∂t P r (x, y, z, t|x, y, 0) =δ(y)δ(z)L FP,x P r (x, y, z, t|x, y, 0) + δ(z)L FP,y P r (x, y, z, t|x, y, 0) with initial condition P r (x, y, z, t = 0|x, y, 0) = δ(x − x 0 )δ(y)δ(z) and zero boundary conditions at infinity.

Fokker-Planck Equation for the Marginal PDF
Again, we utilise the approach elaborated upon in Appendix B. The only difference from the other modes of resetting is the last term, given by rδ(z) ∞ −∞ P r (x, y, z , t|x, y, 0) dz , where the Dirac δ function indicates that the particle is reset only on the main fingers in the z-direction. Therefore, in Laplace space, we obtain 2ĝ(x, y, s) r(x, y, s) .
As in the previous mechanisms of resetting, the expressions (A31) and (A33) have the form Following the same procedure as before, we end up with an equation in Laplace space for the marginal PDFP 1,2 (x, y, s), D z δ(y)L FP,x + L FP,y P 1,2 (x, y, s) + r (s + r) 1/2P 1,2 (x, y, s).
In order to get the equation for the marginal PDF in the x-direction, we present the solutionP 1,2 (x, y, s) in the form From the approach discussed in Appendix B, we obtainq(x, s) in the form of Equation (A41) as followŝ Substituting forĥ(x, s) = 1 2q (x, s)p 1 (x, s), we obtain If we multiply Equation (66) by With the inverse Laplace transform of Equation (67), we obtain the Fokker-Planck equation for the marginal PDF along the backbone, where is the memory kernel.

First Moment and MSD
We calculate the mean value and the MSD by means of the standard procedure described above, from where we obtain for the mean value, and for the MSD. The final MSD can be obtained by numerical inverse Laplace transform in MATHEMATICA [64].

Numerical Simulations: Coupled Langevin Equations
For simulation purposes, here we also employed the coupled Langevin equations approach. Thus, we have for the movement along the backbone, for the movement along the main fingers, and for the movement along the secondary fingers. The typical trajectories along all three directions are given in Figure 5.  (73), where the black dots stand for the resetting events and the shaded and white areas represent the time between two consecutive resetting events, for ∆t = 0.01, r = 1, µ = 3, σ x = σ y = σ z = 1, x 0 = 0, λ x = 1, λ y = 0.001.

Discussion
The results of the numerical calculation of the MSDs and the marginal PDFs obtained in the previous sections are collected in Figures 6-8. The main conclusion is as follows. Resetting strongly affects the O-U process on the three-dimensional comb. In particular, there is a localisation of the particles along the backbone in the presence of resetting, which leads to saturation of the MSDs for all of the scenarios considered in  In what follows, we discuss these scenarios.
From Figures 3-5, one can conclude that the waiting times for the movement in xdirection are longer than those for the y-direction in every type of resetting mechanism, as it is expected, since the time spent by the particle while moving along the fingers and secondary fingers can be considered as a waiting time for the movement along the backbone, while the time spent by the particle in the secondary fingers can be considered as a waiting time for the movement along the main fingers. Of great interest is the observation that the duration of these waiting times for the movement in the x-direction is shortest for the case of global resetting, longer in the case of resetting to the backbone, and longest in the case of resetting to the main fingers. This means that the type of resetting mechanism have significant impact on the time the particle spends being trapped in the main and secondary fingers.
In the left panels of Figure 6, the analytical results for the MSDs are presented, which are confirmed with Monte Carlo simulations, and in the right panels, the Monte Carlo simulations for the marginal PDFs are depicted, using the same parameters as for the MSDs. The analytical results of the MSD in temporal form are computed using the numerical (Stehfest) inverse Laplace transform [64]. The results displayed are for the three different modes of resetting, namely: global resetting, Figure 6a,b, resetting to the backbone, Figure 6c,d, and resetting to the main fingers, Figure 6e,f, for three different rates of resetting r = {0, 1, 5}.
From the results for the MSDs and the PDFs depicted in Figure 6, we can conclude the following. In the case of global resetting to the initial position (x 0 , y 0 , z 0 ), the increase in the resetting rate is characterised by a decrease in the overall MSD, and the particle achieves the stationary state much faster, and from its PDF, it is evident that with a rise in the resetting rate r, the probability of the particle of being found around the initial position also increases. In the event of resetting to the backbone and the main fingers, the opposite is the case. The MSDs of the processes climb with larger resetting rates, but the behaviour of achieving the stationary state remains the same; namely, the speed of establishment of the stationary state is proportional to the increase in the resetting rate, even though the value of the stationary state remains the same regardless of the value of the resetting rate. We can also say that the stationary state is established faster in the case of resetting to the backbone than the one with resetting to the main fingers, as is evident in Figure 6c,e. As can be seen in Figure 6d,f, the marginal PDF for the diffusion on the backbone is not strongly dependent on whether we choose resetting to the backbone or main fingers and exhibits very similar behaviour with both modes of resetting. In both cases, with an increase in the resetting rate r, the probability of the particle of being found near the long-term mean value µ is larger, and the probability of being found in the vicinity of the initial position decreases, something that is completely opposite to the behaviour when we have global resetting.
From Figure 7, it can be concluded that the saturation value of the MSD increases with the increase in the mean reverting rate λ y in the y-direction in the case of global resetting of the system. In Figure 8, we see that the PDFs for the particles to be found at the long-term mean value µ also rise with the increase in the mean-reverting rate λ y for all three modes of resetting.   We set x 0 = 0, σ x = σ y = σ z = 1, µ = 5, λ x = 3, r = 1, N = 10 3 , dt = 0.001, t = 10 3 .

Summary
In this work, we analyse the O-U processes on a three-dimensional comb without and with resetting. Three different resetting mechanisms were used: global resetting, resetting to the backbone, and resetting to the fingers. From the corresponding Fokker-Planck equation, analytical results of the MSD are obtained, and the PDF is analysed in Laplace space from where one can obtain the NESS, which is reached in the long time limit. Numerical simulations of the MSD and PDF are performed with the help of coupled Langevin equations for the comb geometry. The influence of the comb geometry, resetting rate and the bounding potential on the particle dynamics is analysed. Due to the application of the Ornstein-Uhlenbeck process in the description of diffusing particles in harmonic potential in physics, as well as in the description of interest and currency exchange rates in economics, we believe that the considered models on comb structures can be of interest in the description of the corresponding anomalous diffusive processes in complex systems that occur as a result of particle trapping, long-tailed waiting times, and resetting mechanisms.

Funding:
The Authors acknowledge financial support by the German Science Foundation (DFG, Grant number ME 1535/12-1). This work is also supported by the Alliance of International Science Organizations (Project No. ANSO-CR-PP-2022-05). TS was supported by the Alexander von Humboldt Foundation.

Data Availability Statement:
No new data have been produced for this work.

Acknowledgments:
The authors thank Alexander Iomin for helpful suggestions and critical reading of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: where E γ α,β (z) is a three-parameter M-L function (A7). For r = 0 it becomes the Prabhakar derivative in Riemann-Liouville form.

Appendix B. Solving the Corresponding Equations
Let us analyse the following differential equation.