Stochastic resetting on comblike structures

.


I. INTRODUCTION
Combs are two-or three-dimensional branched structures with a backbone crossed by perpendicular fingers.These fingers may be one or two dimensional side structures.A random walker moving along the backbone may enter into a finger (or fingers) and move there for a time and return to the backbone to start the process again.As a result of a Brownian motion in a two-dimensional comb the mean squared displacement (MSD) shows a subdiffusive behavior depending on time as t 1/2 and was originally introduced to understand anomalous transport in percolation clusters and many other applications [1].Although three-dimensional combs have been less developed, they have modeled transport in spiny dendrites [2] or ultraslow diffusion in combs with circular fingers [3].
In this paper we investigate a diffusion process on a xyz-comb (see Fig. 1) with stochastic resetting.The diffusion on a three dimensional comb is governed by the following equation ∂ ∂t P (x, y, z, t) = L FP P (x, y, t), where is the Fokker-Planck (transport) operator, and D x δ(y)δ(z), D y δ(z) and D z are the diffusion coefficients along the x, y and z directions, respectively.The δ-functions δ(y)δ(z) in front of the second spatial derivative with respect to x, mean that diffusion along the backbone (x-axis) is allowed only at y = z = 0, while the δ-function δ(z) in front of the second spatial derivative with respect to y, means that the diffusion along the main fingers (or branches) (y-axis) is allowed only at z = 0.The z-axis is a secondary finger, an auxiliary direction along which the particle performs normal diffusion.
On the other hand, the diffusion process in one dimension under stochastic resetting was introduced by Evans and Majumdar [4].The corresponding equation is given by where the initial position reads P (x, t = 0|x 0 ) = δ(x − x 0 ), D is the diffusion coefficient, r is the rate of resetting arXiv:2005.11575v1[cond-mat.stat-mech]23 May 2020 FIG. 1: Three-dimensional comb-like structure, which is a discrete caricature of the continuous 3D comb model described by Eq. ( 1).It consists of the backbone along the x axis and continuously distributed side-branchesfingers along the y and z axes.
to the initial position, the second term on the right hand side represents the loss of probability from the position x due to reset to the initial position x 0 , and the third term is the gain of probability at x 0 due to resetting from all other positions.This equation represents a renewal process: each resetting event to the initial position x 0 renews the process at a rate r.Between two consecutive renewal events, the particle undergoes free diffusion [4].
Multi-dimensional diffusion has already been studied in the literature [6].There, diffusion with resets in an multi-dimensional, homogeneous, infinite media is studied.In this work we analyze the transport properties and the long time behavior of diffusion in an heterogeneous environment and determine the properties emerging from resetting in the different space coordinates.
The paper is organized as follows.In Section II we consider diffusion in a three-dimensional comb with global exponential (Markovian) resetting.We give exact results for the marginal probability density functions (PDFs), stationary distributions and MSDs along all three axes.
We also confirm the analytical results by numerical simulations by employing a Langevin equation approach for comb structure.Excellent agreement has been shown.Diffusion in three-dimensional comb with exponential resetting to the backbone is considered in Section III and the corresponding PDFs and MSDs are also found.In Section IV exponential resetting to the fingers is analyzed.We also discuss the resetting mechanisms in two-dimensional comb-structures in Section V.In Section VI we give detailed explanation of the topological constraint of the transport properties of both two-and threedimensional comb structures.Summary is provided in Section VII.

II. GLOBAL RESETTING A. Analytical results
We start our analysis by considering diffusion in a three dimensional comb with global resetting, represented by the equation with the initial position P (x, y, z, t = 0|x 0 ) = δ(x − x 0 )δ(y)δ(z).This equation can also be interpreted in terms of a renewal process: each resetting event to the initial position (x 0 , y 0 , z 0 ) = (x 0 , 0, 0) renews the process at a rate r.Between two consecutive renewal events, the particle undergoes diffusion on the xyz-comb structure.
To find the solution of Eq. ( 3) we apply the Fourier transformations * with respect to x, y and z, and the Laplace transformation † with respect to t.Therefore, for the PDF in the Fourier-Laplace domain we obtain, see Supplemental material 1 for details of calculations,

B. Marginal PDFs
In order to analyze the motion along all three directions, we analyze the marginal PDFs, In the Fourier-Laplace space, the marginal PDFs are p2 (k y , s|0) = P (k x = 0, k y , k z = 0, s|x 0 , 0, 0), ( 9) Therefore, from Eqs. ( 4) and ( 8), for the marginal PDF along the backbone we have where . By applying the inverse Fourier transform we obtain From Eq. (11), by inverse Fourier-Laplace transforms we arrive at the generalized (non-Markovian) diffusion equation along the backbone where Here is the three parameter Mittag-Leffler function [40]  ‡ and (δ) k = Γ(δ + k)/Γ(δ) is the Pochhammer symbol.The ‡ The Laplace transform of the three parameter Mittag-Leffler function reads [40] L Its asymptotic behaviors are given by [38,39] initial condition is p 1 (x, t = 0|x 0 ) = δ(x − x 0 ).The equation can also be written in the form where is a so-called regularized Prabhakar derivative [41], which has many applications nowadays [42][43][44].The stationary PDF along the backbone, obtained in the long time limit, becomes From Eqs. ( 4) and ( 9), for the PDF along the fingers we find where Dz .From here, by applying the inverse Fourier transform we obtain p2 (y, s|0) = (s + r) The marginal PDF p 2 (y, t|0) provides the transport equation along the main finger and it is governed by the equation with the initial condition p 2 (y, t = 0|0) = δ(y).Here the kernel ζ(t) is Eq. ( 21) can also be presented by means of the regularized Prabhakar derivative (17).It reads or equivalently where erf(z) = 2 √ π z 0 e −t 2 dt is the error function, while is the tempered Caputo derivative with the exponential truncation, where b > 0 is the truncation parameter [37,42].For the stationary PDF along the y direction, we find For the z direction, we have that yields where D 3 = D z .The corresponding equation for the transport along secondary fingers reads where and it can be rewritten in the equivalent form which is a diffusion-absorption equation.The stationary PDF along the z direction is It is interesting to note that the obtained stationary PDFs along each axis are exponential functions, so that the global resetting does not modify the stationary state with respect to the case of standard diffusion in one dimension under resetting, see Eq. ( 2).The global resetting affects the transient dynamics towards the stationary state only.For a given value of D for all three components, the width of the PDF varies.This can be seen from the analytical formulas, in which the width of the exponential stationary PDF depends on r for the z component, on √ r for the y component and on 4 √ r for the x component, see Eqs. ( 18), ( 26), (32).The obtained analytical results are verified by stochastic simulations with the Langevin equation approach, see Section II D.

C. Mean squared displacements
In the section, we analyze the MSDs along all three directions, Taking into account corresponding solutions for the marginal PDFs, we have where t n dt is the exponential integral function.This corresponds to the transition from subdiffusion to localization For the y fingers we have that corresponds to the transition from subdiffusion to localization as well, Eventually, the MSD for z-fingers reads that corresponds to saturation in the long time limit, Therefore, unlike an initial transient behaviour all the MSDs saturate towards a constant value (exhibiting stochastic localization) as in the scenario of onedimensional diffusion with resets.This confirms the existence of a non-equilibrium stationary state, which has been recently observed for many different dynamics under constant-rate resets [22].This variety of cases has been also obtained from stochastic simulations of the process based on a Langevin equation approach, see Section II D.
FIG. 2: Trajectory along individual axes with global stochastic resetting to the initial position (x 0 , y 0 , z 0 ) = (0, 0, 0) with rate r = 0.002 obtained from a Langevin simulation of the process.The resetting events are represented by black dots.Dashed regions are introduced in order for these resetting events to be more visible.

D. Langevin equation approach. Numerical simulations
To verify the analytical solution obtained in the previous section, we perform numerical calculations, considering a system of Langevin equations [45,46], and where resets, as a renewal process, can be easily performed, see Ref. [47].The system of coupled Langevin equations reads where As a result, the noise in Eqs.(36a) and (36b) is multiplicative, however in Refs.[46,48] the authors verified that the value ε has no influence in the diffusive process, as long as ε and the noise amplitudes β 1 , β 2 , β 3 are of the same order of magnitude.In our simulations, we have set ε = β x = β y = β z = 0.1.The noises η x (t), η y (t), η z (t), were sampled from a Gaussian distribution N (0, ∆t).The time evolution of the diffusive particle is a renewal process, where each resetting event to (x 0 , y 0 , z 0 ) renews the process at a Poisson rate r.
This effect of stochastic resetting is modeled by sampling a resetting time from an exponential distribution with parameter r representing the time between two events in a Poisson point process.During this resetting time, the particle undergoes diffusion on the threedimensional comb and resets at (x 0 , y 0 , z 0 ) afterwards.Graphical representation of the simulations of particle trajectories along all directions is given in Fig. 2.
Regarding the simulation of marginal PDFs and temporal evolution of the variance, ensembles of 5 × 10 4 particle positions were simulated considering a time step of ∆t = 1 across a time span of 10 5 in order to observe convergence of the processes, with the MSD being calculated for each of the ensembles along all three directions σ 2 In Fig. 3 we give comparison of the analytical and simulation results for the marginal PDFs, where we use that , with ∆t = 1, see [45,47].In Fig. 4 we show the simulated time evolution of the MSDs in the three directions.From the simulation results one can verify that they are in a very good agreement with the saturation values obtained analytically, if one uses that 2 D i = β 2 i .For more results on the corresponding PDFs obtained by the numerical simulations we refer to the Supplemental material 1. FIG. 4: MSDs along all three axes for global resetting with rate r = 0.05.Each color corresponds to an axis, x (blue), y (green), z (red) with and without resetting.

III. RESETTING TO THE BACKBONE
Global resetting takes the particle to a particular position of the comb (with the coordinates (x 0 , 0, 0) in Eq. ( 3)).However, it is only one of many possible mechanisms of resetting.Here, we proceed with a slightly softer resetting procedure, which takes the particle to the backbone.This resetting is applied to the y and z directions only, taking a walker being at (x, y, z) to the point (x, 0, 0).In this case, the governing equation reads ∂ ∂t P (x, y, z, t|x0, 0, 0) = LFPP (x, y, z, t|x0, 0, 0) which differs from Eq. ( 3) in the last term only.This difference results from the difference between the global resetting and resetting to the backbone.In the former case the particle is taken at the particular position (x 0 , 0, 0) as stated by the δ(x−x 0 )δ(y)δ(z) term in Eq. ( 3).However, in the latter case, considered here, the particle appears at y = 0, z = 0 but the x position is not modified.Mathematically it can be written as the marginal distribution, the double integral term in Eq. ( 72).From the Fourier-Laplace transformations, we arrive at the following PDF in the Fourier-Laplace space, see Supplemental material 2 for details, From this equation, the marginal PDFs for the three axis can be straightforwardly obtained as done in the previous section: Here we note that the corresponding equations, Eqs. ( 40) and ( 41), for the marginal PDFs along the y and z directions are the same as in the case of global resetting, Eqs. ( 20) and (27), respectively.Along the backbone the PDF is which is governed by the equation where T C D α b f (t) is the tempered fractional derivative (25) of order 1/4.
The corresponding MSDs along y and z axes are the same as those for the case of global resetting, since the effect of the resetting in these two dimensions is equivalent for both scenarios.However, the dynamics on the x axis change substantially as reflected in the MSD FIG. 5: Trajectory along individual axes with stochastic resetting to the backbone with rate r = 0.002 obtained from a Langevin simulation of the process.The resetting events are represented by black dots.Note that there is no resetting along x axis but the dashed regions are given only for indication when the resetting events along y and z axis occur.
Its asymptotics reads Γ(5/4) , rt 1, . The resetting mechanism studied in this section enhances the transport since it returns particles to the x axis.Consequently, instead of the saturation of a stationary value for the MSD, one can see from Eq. (91) that in the long time limit, x 2 (t) ∼ t, i.e. it scales diffusively.The short time limit scales as x 2 (t) ∼ t 1/4 , as in the case of global resetting.This means that we observe accelerating transport along the backbone, ranging from subdiffsuion to normal diffusion.The numerical simulations of particle trajectories along all three directions, by using the Langevin equations approach, are shown in Fig. 5.We see that while in the y and z axis we observe recurrent returns to the origin, in the x axis the motion does not return.Instead, it freely moves away from the origin.In Fig. 6 we give comparison of the analytical and simulation results for the marginal PDFs from where one observes excellent agreement between both approaches.Same parameters for ε, β i and D i as in the case of global resetting are used.The analytical results has also been confirmed by simulation results of the MSDs given in Fig. 7, where the y and z components of the MSD reach a stationary value while the MSD in the x direction increases linearly.This is in agreement with the analytical results found above.More simulation results for the PDFs are given in the Supplemental material 2.

IV. RESETTING TO THE MAIN FINGERS
Finally, we study the dynamics of the system when the resetting applies to particle located at any secondary finger along the z axis and moves to main finger (axis y in Fig. 1).In this case, the governing equation reads ∂ ∂t P (x, y, z, t|x0, 0, 0) = LFPP (x, y, z, t|x0, 0, 0) where the last term is now the marginal distribution in the variables x and y.In the Fourier-Laplace space, the solution of the equation reads, see Supplemental material 3 for details, This also yields the images of the marginal PDFs for the different axis.Eventually, we have From these expressions one obtains the corresponding MSDs, which are and the MSD along the z axis is the same as in the previous cases.Their asymptotics read Γ(5/4) , rt 1, In this case, the MSD in the x-axis behaves subdiffusively with x 2 (t) ∼ t 1/4 as in the case with no resetting, and then it turns to x 2 (t) ∼ t 1/2 , which means an accelerating subdiffusive transport.Along the y-axis, the MSD scales as y 2 (t) ∼ t 1/2 in the short time limit, and then it turns to linear dependence in time y 2 (t) ∼ t.Along the z-axis the MSD from the normal diffusive behavior reaches a stationary value in the long time limit.We also performed numerical simulations by using the Langevin equations approach.Same parameters for ε, β i and D i as previous are used.The simulation results show very good agreement with the analytical results, see Figs. 8, 9 and 10.For more details on analytical computation and simulation results see also Supplemental material 3.

V. REMARKS ON TWO DIMENSIONAL COMB
Here we note that the results obtained for the three-dimensional xyz-comb can be used for the twodimensional xy-comb.The y and z axes in the threedimensional comb would correspond to the x and y axes in the two-dimensional comb.Therefore, the results obtained for the y and z directions in the three-dimensional comb with resetting in the backbone correspond to the results for the x and y directions in the two-dimensional comb with global resetting.Furthermore, the results obtained for the y and z directions in the three-dimensional comb with resetting in the main fingers correspond to the results for the x and y directions in the two-dimensional comb with resetting in the fingers.

VI. REMARK ON THE THREE DIMENSIONAL COMB GEOMETRY
The topological (comb) constraint of the transport properties of both two-and three-dimensional combs should be discussed as well.To that end, let us understand the role of the δ(y) and δ(z) functions in the highly inhomogeneous diffusion coefficients in Eq. ( 1).One should recognize that the singularity of the x and y components of the diffusion tensor is the intrinsic transport property of the comb model (1).Note that this singularity of the diffusion coefficients relates to a nonzero flux along the x backbone and y fingers, and for the two dimensional case it was discussed in Refs.[49][50][51][52].
Here, we extend the arguments of Refs.[51,52] for the three dimensional case of Eq. ( 1).Let us consider the Liouville equation where the three dimensional current j = (j x , j y , j z ) describes Markov processes in Eq. ( 1).In this case, the three-dimensional current reads − /2 dz . . ., one obtains for the l.h.s. of the equation, after application of the middle point theorem, 2 ∂ ∂t P (x, y = 0, z = 0, t), which is exact in the limit → 0. This term can be neglected in this limit → 0. Considering integration of the r.h.s. of the equation, one should bear in mind that this procedure is artificial and its implementation needs some care.First we consider the currents outside of the vicinity of the x backbone.In this case, according to the comb geometry, j x = 0 and we consider a two-dimensional y − z comb.Therefore, we perform integration with respect to z only.From Eq. (56c) we obtain that this term responsible for the transport in the z direction reads where prime means derivative with respect to z.This corresponds to the two outgoing fluxes from the y fingers in the ±z directions: F z (+) + F z (−).The transport in the y direction, after integration, is It should be stressed that the second derivative over y, presented in the form ensures both incoming and outgoing fluxes for F y along the y direction at a point y.Following the Kirchhoff's law, we have F y + F z (+) + F z (−) = 0 for every point y and at z = 0. Function F y contains both incoming and outgoing fluxes of the probability, while F z (+) and F z (−) are both outgoing probability fluxes.If the latter outgoing fluxes are not zero, the flux F y has to be nonzero as well: F y = 0, as containing an incoming flux.Therefore, D y (z → 0) = 0. Taking D y (z) = 1 π Dy y 2 + 2 , one obtains in the limit → 0 a nonzero flux F y with D y (z) = D y δ(z), which is the diffusion coefficient in the y direction in Eq. ( 55).Now we perform integration in the vicinity of the x backbone, where we take into account the singularity of the y component of diffusion coefficient, which is D y δ(z).We also admit that integration of the j z current in Eq. (56c) with respect to y yields zero.Therefore, integration with respect to y and z yields from Eq. (56b) the term responsible for the transport in the y direction as follows D y P (x, y, z = 0, t) y= 2 − P (x, y, z = 0, t) y=− 2 .
Here prime means derivative with respect to y.This corresponds to the two outgoing fluxes from the backbone in the ±y directions: F y (+)+F y (−).The transport along the x direction, after integration of Eq. (56a), is In complete analogy with the y coordinate, the second derivative with respect to x, presented in the form as → 0, ensures both incoming and outgoing fluxes for F x along the x direction at a point x.Again, after the integration, the Liouville equation is a kind of the Kirchhoff's law: F x (+) + F x (−) + F y (+) + F y (−) = 0 for each point x and at y = 0. Note, that the flax in the z direction is zero due to the integration with respect to y.Since outgoing fluxes are not zero, j x = 0 and correspondingly the flux F x ≡ F x (+) + F x (−) has to be nonzero as well: F x (±) = 0. Therefore, 2 D(y → 0, z → 0) = 0. Now, taking the diffusion coefficient in the form 2 , one obtains in the limit → 0 a nonzero flux F x with D(y, z) = D x δ(y)δ(z), which is the diffusion coefficient in the x direction in Eqs. ( 1), (55) and (56a).

VII. SUMMARY
We have studied the dynamics of a particle, which diffusing in a three-dimensional heterogeneous comb-like structure performs different types of resets.The hierarchical structure of the three-dimensional comb allows to study different resetting mechanisms that generate a wide variety of dynamics depending on the strength of the resetting mechanism.In particular, we studied thee types of resets in the three-dimensional comb and their influence on the dynamics of the MSD, and we found that at the short time there is no influence on the transport exponents, which remain the same as in the case without resetting.However, at the long time limit, the system is strongly affected by the resetting process which leads a change in the transport exponents.
We studied three kinds of resetting: global resetting of a particle from any point on the comb to a fixed point at (x, y, z) = (x 0 , 0, 0) and two kinds of softer resetting, where two coordinates (y = 0, z = 0) and one coordinate (z = 0) are fixed.When resets are global, the MSDs in x, y and z directions reach constant values exhibiting stochastic localization, i.e. a non-equilibrium steady state is reached.This result is in complete agreement with the results observed for the dynamics of walkers with constant rate resetting recently studied in the literature [4,22].For a softer version of resetting consisting with two fixed coordinates, the walker returns to any positions at the backbone.It means that if the position of the walker before the resetting is (x, y, z), then after the reset it is (x, 0, 0).In this case, the dynamics for the y and z axis remains the same as in the global resetting case, since the effect of the resetting to these two coordinates is the same.However, in the x direction, the resetting enhances the motion: it becomes subdiffusive x 2 (t) ∼ t 1/4 for short times and then normal diffusion takes place for the long time scale.The latter regime results from the fact that the mean waiting time to stay in (y, z) fingers becomes finite due to resetting.Indeed, the reset time PDF now plays the role of a waiting time PDF for the motion along the backbone (x direction).Since the reset time PDF is exponential (i.e., constant rate resetting or Markovian resetting process), the motion in the x direction becomes diffusive in the long time limit.For the softer resetting with one fixed z-coordinate, a stationary regime takes place in the z fingers only.In the x and y directions the transport is enhanced with the MSDs behaving as follows: x 2 (t) ∼ t 1/2 and y 2 (t) ∼ t.The obtained diffusion equations for the marginal PDFs shed light on the physical relevance of usage of the regularized Prabhakar derivative in diffusion theory.They could describe diffusion processes on comb structures with stochastic resetting.
Therefore, we finally obtain (65) After we find the marginal PDFs we will analyze the mean squared displacements along all three directions, ) dy, and y 2 (t) = ∞ −∞ z 2 p 3 (z, t|0) dz.Therefore, we have where is the three parameter Mittag-Leffler function, and (δ t n dt is the exponential integral function, where erf(z) = 2 √ π z 0 e −t 2 dt gives the error function, and The Laplace transform of the three parameter Mittag-Leffler function (67) reads and its asymptotic behaviors are given by (71) In Figs.11, and 12, we give graphical representations of the PDFs of the process without resetting and under global resetting.The results are obtained by numerical simulations similar to the ones elaborated in Sec.II D, with different values of the ensemble size of 2 × 10 4 and the time span of 5 4 .

Solutions for diffusion with reseting to the backbone
Previously we have studied a resetting mechanism which takes the particle to a particular position of the comb, the coordinate (x 0 , 0, 0).Here, we proceed to analyse a slightly different resetting procedure, consisting on taking the particle to the backbone.This is, if the particle is at any position (x, y, z), the resetting consists on taking it to the position (x, 0, 0) of the state space.In this case, the Fokker-Planck equation reads ∂ ∂t P (x, y, z, t|x 0 , 0, 0) = D x δ(y)δ(z) ∂ 2 ∂x 2 P (x, y, z, t|x 0 , 0, 0) + D y δ(z) ∂ 2 ∂y 2 P (x, y, z, t|x 0 , 0, 0) + D z ∂ 2 ∂z 2 P (x, y, z, t|x 0 , 0, 0) Here, instead of appearing at the particular position (x 0 , 0, 0) as stated by δ(x − x 0 )δ(y)δ(z), the particle appears at y = 0, z = 0 but the x position is not modified and, therefore, it can be mathematically written as the marginal distribution (double integral term in Eq. ( 72)).In order to solve the equation for the propagator P (x, y, z, t|x 0 , 0, 0), we apply the Fourier-Laplace transform as done in the previous case to obtain s P (k x , k y , k z , s|x 0 , 0, 0) − e ikxx0 = −D x k 2 x P (k x , y = 0, z = 0, s|x 0 , 0, 0) − D y k 2 y P (k x , k y , z = 0, s|x 0 , 0, 0) − D z k 2 z P (k x , k y , k z , s|x 0 , 0, 0) − r P (k x , k y , k z , s|x 0 , 0, 0) + r P (k x , k y = 0, k z = 0, s|x 0 , 0, 0).
Let us proceed to find the three unknown terms on the right hand side.Taking k y = k z = 0 one can isolate P (k x , k y = 0, k z = 0, s|x 0 , 0, 0) = e ikxx0 s − D x k 2 x P (k x , y = 0, z = 0, s|x 0 , 0, 0) s , In Fig. 13, we give graphical representations of the PDFs of the process under resetting to the backbone.

Solutions for diffusion with reseting to the main fingers
Finally, we study the dynamics of the system when, instead of taking the particle to a fixed position in the x axis or to the backbone, the resetting only applies to the z axis, taking the particle to the corresponding main finger.In In Fig. 14, we give graphical representations of the PDFs of the process under resetting to the backbone.
(t + ∆t) =    y(0), with prob.r ∆t, y(t) + β2B(z)ηy(t), with prob.(1 − r ∆t), (36b) y) and B(z) are functions introduced to mimic δfunctions (see Refs.[45,48]), and r is the parameter of the Poisson process.To replicate the Dirac δ-function, diffusion across the x and y directions is permitted in a narrow band of thickness 2ε along the x and y axes.

FIG. 8 :FIG. 9 :
FIG.8: Trajectory along individual axes with stochastic resetting to the fingers with rate r = 0.002 obtained from a Langevin simulation of the process.The resetting events are represented by black dots.Note that there is no resetting along x and y axis but the dashed regions are given only for indication when the resetting events along z axis occur.

FIG. 14 :z 2
FIG. 14: Non-normalized PDF along (a) x direction, (b) y direction, and (c) z direction, for resetting to the fingers with rate r = 0.002.