Stochastic resetting on comb-like structures^*^**V. D. and A. M. contributed equally to this work.

Viktor Domazetoski v.domazetoski@manu.edu.mk Research Center for Computer Science and Information Technologies, Macedonian Academy of Sciences and Arts, Bul. Krste Misirkov 2, 1000 Skopje, Macedonia Axel Masó-Puigdellosas Axel.Maso@uab.cat Grup de Física Estadística, Departament de Física. Universitat Autònoma de Barcelona. Edifici Cc. E-08193 Cerdanyola (Bellaterra) Spain Trifce Sandev trifce.sandev@manu.edu.mk Research Center for Computer Science and Information Technologies, Macedonian Academy of Sciences and Arts, Bul. Krste Misirkov 2, 1000 Skopje, Macedonia Institute of Physics & Astronomy, University of Potsdam, D-14776 Potsdam-Golm, Germany Institute of Physics, Faculty of Natural Sciences and Mathematics, Ss. Cyril and Methodius University, Arhimedova 3, 1000 Skopje, Macedonia Vicenç Méndez Vicenc.Mendez@uab.cat Grup de Física Estadística, Departament de Física. Universitat Autònoma de Barcelona. Edifici Cc. E-08193 Cerdanyola (Bellaterra) Spain Alexander Iomin iomin@physics.technion.ac.il Department of Physics, Technion, Haifa 32000, Israel Ljupco Kocarev lkocarev@manu.edu.mk Research Center for Computer Science and Information Technologies, Macedonian Academy of Sciences and Arts, Bul. Krste Misirkov 2, 1000 Skopje, Macedonia Faculty of Computer Science and Engineering, Ss. Cyril and Methodius University,
P.O. Box 393, 1000 Skopje, Macedonia

(Submitted 12 Dec 2019, Revised 25 Apr 2020)

Abstract

We study a diffusion process on a three-dimensional comb under stochastic resetting. We consider three different types of resetting: global resetting from any point in the comb to the initial position, resetting from a finger to the corresponding backbone and resetting from secondary fingers to the main fingers. The transient dynamics along the backbone in all three cases is different due to the different resetting mechanisms, finding a wide range of dynamics for the mean squared displacement. For the particular geometry studied herein, we compute the stationary solution and the mean square displacement and find that the global resetting breaks the transport in the three directions. Regarding the resetting to the backbone, the transport is broken in two directions but it is enhanced in the main axis. Finally, the resetting to the fingers enhances the transport in the backbone and the main fingers but reaches a steady value for the mean squared displacement in the secondary fingers.

comb, anomalous diffusion, stochastic resetting

pacs:

87.19.L-, 05.40.Fb, 82.40.-g

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 book . Although three-dimensional combs have been less developed, they have modeled transport in spiny dendrites MeIo13 or ultraslow diffusion in combs with circular fingers IoMe16 .

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

\displaystyle\frac{\partial}{\partial t}P(x,y,z,t)=L_{\mathrm{FP}}P(x,y,t),

(1)

where

L_{\mathrm{FP}}=\mathcal{D}_{x}\delta(y)\delta(z)\frac{\partial^{2}}{\partial x^{2}}+\mathcal{D}_{y}\delta(z)\frac{\partial^{2}}{\partial y^{2}}+\mathcal{D}_{z}\frac{\partial^{2}}{\partial z^{2}}

is the Fokker-Planck (transport) operator, and $\mathcal{D}_{x}\delta(y)\delta(z)$ , $\mathcal{D}_{y}\delta(z)$ and $\mathcal{D}_{z}$ are the diffusion coefficients along the $x$ , $y$ and $z$ directions, respectively. The $\delta$ -functions $\delta(y)\delta(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 $\delta$ -function $\delta(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.

Refer to caption — Figure 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-branches - fingers along the $y$ and $z$ axes.

On the other hand, the diffusion process in one dimension under stochastic resetting was introduced by Evans and Majumdar maj1 . The corresponding equation is given by

\displaystyle\frac{\partial}{\partial t}P(x,t|x_{0})=\mathcal{D}\frac{\partial^{2}}{\partial x^{2}}P(x,t|x_{0})-rP(x,t|x_{0})+r\delta(x-x_{0}),

(2)

where the initial position reads $P(x,t=0|x_{0})=\delta(x-x_{0})$ , $\mathcal{D}$ is the diffusion coefficient, $r$ is the rate of resetting 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 maj1 . It is known that this equation has a stationary solution in the long time limit given by

P_{\mathrm{st}}(x|x_{0})=\frac{1}{\sqrt{4\mathcal{D}/r}}e^{-\frac{|x-x_{0}|}{\sqrt{\mathcal{D}/r}}}.

Afterwards, other types of motion and resetting mechanisms have been studied by introducing the two resetting terms to the Fokker-Planck equation of the corresponding process maj2 ; EvMa14 ; maj3 ; DuHe14 ; Pa15 ; ChSc15 ; Ma19 ; Gu19 ; lenzi_front . For instance, space-dependent reset rates maj2 or diffusion in a potential landscape Pa15 and the telegraphers equation Ma19 have been analyzed under this perspective. Other works have studied motion with resetting by employing a renewal equation MoVi13 ; mendez ; pal1 ; NaGu16 ; EuMe16 ; MoMaVi17 ; Sh18 ; EvMa18 ; axel ; MaCaMe19p ; BoChSo19 ; BoChSo19p ; MaCaMe19pp ; MaCaMe19ppp , which has also been used to study the completion time of search processes with resetting CaMe15 ; reu ; pal reu ; chechkin ; pal pre ; pal prr ; masoliver2019 .

Multi-dimensional diffusion has already been studied in the literature EvMa14 . 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 three-dimensional comb structures. Summary is provided in Section VII.

II Global resetting

II.1 Analytical results

We start our analysis by considering diffusion in a three dimensional comb with global resetting, represented by the equation

		$\displaystyle\frac{\partial}{\partial t}P(x,y,z,t\|x_{0},0,0)=L_{\mathrm{FP}}P(x,y,z,t\|x_{0},0,0)$
		$\displaystyle-rP(x,y,z,t\|x_{0},0,0)+r\delta(x-x_{0})\delta(y)\delta(z),$		(3)

with the initial position $P(x,y,z,t=0|x_{0})=\delta(x-x_{0})\delta(y)\delta(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. (II.1) we apply the Fourier transformations^†^††The Fourier transform of a function $f(\xi)$ is given by $f(k)=\mathcal{F}[f(x)](k)=\int_{-\infty}^{\infty}f(\xi)\,e^{\imath k\xi}\,d\xi$ . The inverse Fourier transform then reads $f(\xi)=\mathcal{F}^{-1}[f(k)](x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}f(k)\,e^{-\imath k\xi}\,dk$ . with respect to $x$ , $y$ and $z$ , and the Laplace transformation^‡^‡‡The Laplace transform of a function $f(t)$ reads $\hat{f}(s)=\mathcal{L}[f(t)](s)=\int_{0}^{\infty}f(t)\,e^{-st}\,dt$ . with respect to $t$ . Therefore, for the PDF in the Fourier-Laplace domain we obtain, see Supplemental material 1 for details of calculations,

		$\displaystyle\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)=\frac{1}{s}\times\frac{(s+r)^{1/4}}{(s+r)^{1/4}+\frac{\mathcal{D}_{x}}{2\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}k_{x}^{2}}$
		$\displaystyle\times\frac{(s+r)^{1/2}}{(s+r)^{1/2}+\frac{\mathcal{D}_{y}}{2\sqrt{\mathcal{D}_{z}}}k_{y}^{2}}\times\frac{(s+r)}{(s+r)+\mathcal{D}_{z}k_{z}^{2}}\times e^{\imath k_{x}x_{0}}.$

II.2 Marginal PDFs

In order to analyze the motion along all three directions, we analyze the marginal PDFs,

	$\displaystyle p_{1}(x,t\|x_{0})=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P(x,y,z,t\|x_{0},0,0)\,dy\,dz$		(5)
	$\displaystyle p_{2}(y,t\|0)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P(x,y,z,t\|x_{0},0,0)\,dx\,dz,$		(6)
	$\displaystyle p_{3}(z,t\|0)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P(x,y,z,t\|x_{0},0,0)\,dx\,dy.$		(7)

In the Fourier-Laplace space, the marginal PDFs are

\displaystyle\hat{p}_{1}(k_{x},s|x_{0})=\hat{P}(k_{x},k_{y}=0,k_{z}=0,s|x_{0},0,0),

(8)

\displaystyle\hat{p}_{2}(k_{y},s|0)=\hat{P}(k_{x}=0,k_{y},k_{z}=0,s|x_{0},0,0),

(9)

\displaystyle\hat{p}_{3}(k_{z},s|0)=\hat{P}(k_{x}=0,k_{y}=0,k_{z},s|x_{0},0,0).

(10)

Therefore, from Eqs. (II.1) and (8), for the marginal PDF along the backbone we have

\displaystyle\hat{p}_{1}(k_{x},s|x_{0})=\frac{s^{-1}(s+r)^{1/4}}{(s+r)^{1/4}+\mathcal{D}_{1}\,k_{x}^{2}}\,e^{\imath k_{x}x_{0}},

(11)

where $\mathcal{D}_{1}=\frac{\mathcal{D}_{x}}{2\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}$ . By applying the inverse Fourier transform we obtain

\displaystyle\hat{p}_{1}(x,s|x_{0})=\frac{1}{2\sqrt{\mathcal{D}_{1}}}s^{-1}(s+r)^{1/8}e^{-\frac{(s+r)^{1/8}}{\sqrt{\mathcal{D}_{1}}}|x-x_{0}|}.

(12)

From Eq. (11), by inverse Fourier-Laplace transforms we arrive at the generalized (non-Markovian) diffusion equation along the backbone

\displaystyle\int_{0}^{t}\gamma(t-t^{\prime})\frac{\partial}{\partial t^{\prime}}p_{1}(x,t^{\prime}|x_{0})\,dt^{\prime}=\mathcal{D}_{1}\frac{\partial^{2}}{\partial x^{2}}p_{1}(x,t|x_{0}),

(13)

where

\displaystyle\gamma(t)=t^{1/4}E_{1,5/4}^{1/4}(-r\,t).

(14)

Here

E_{\alpha,\beta}^{\delta}(z)=\sum_{k=0}^{\infty}\frac{(\delta)_{k}}{\Gamma(\alpha k+\beta)}\frac{z^{k}}{k!}

(15)

is the three parameter Mittag-Leffler function prabhakar ^§^§§The Laplace transform of the three parameter Mittag-Leffler function reads prabhakar $\mathcal{L}\left[t^{\beta-1}E_{\alpha,\beta}^{\delta}(\pm at^{\alpha})\right](s)=\frac{s^{\alpha\delta-\beta}}{\left(s^{\alpha}\mp a\right)^{\delta}},\quad\Re(s)>|a|^{1/\alpha}.$ Its asymptotic behaviors are given by garrappa ; pre2015 $E_{\alpha,\beta}^{\gamma}(-z^{\alpha})\simeq\left\{\begin{array}[]{ll}\frac{1}{\Gamma(\beta)}\exp\left(-\gamma\frac{\Gamma(\beta)}{\Gamma(\alpha+\beta)}z^{\alpha}\right),\quad z\ll 1.\\ \frac{z^{-\alpha\gamma}}{\Gamma(\beta-\alpha\gamma)},\quad z\gg 1.\end{array}\right.$ and $(\delta)_{k}=\Gamma(\delta+k)/\Gamma(\delta)$ is the Pochhammer symbol. The initial condition is $p_{1}(x,t=0|x_{0})=\delta(x-x_{0})$ . The equation can also be written in the form

\displaystyle{{}_{C}}\mathbf{D}_{1,-r,0+}^{1/4,1/4}p_{1}(x,t|x_{0})=\mathcal{D}_{1}\frac{\partial^{2}}{\partial x^{2}}p_{1}(x,t|x_{0}),

(16)

where

\displaystyle{{}_{C}}\mathbf{D}_{\rho,-\nu,0+}^{\delta,\mu}f(t)=\int_{0}^{t}(t-t^{\prime})^{-\mu}E_{\rho,1-\mu}^{-\delta}\left(-\nu t^{\rho}\right)\frac{df(t^{\prime})}{dt^{\prime}}\,dt^{\prime},

(17)

is a so-called regularized Prabhakar derivative polito , which has many applications nowadays sandev mathematics ; epl ; weron . The stationary PDF along the backbone, obtained in the long time limit, becomes

\displaystyle p_{1,\mathrm{st}}(x|x_{0})=\frac{1}{\sqrt{4\mathcal{D}_{1}/\sqrt[4]{r}}}e^{-\frac{|x-x_{0}|}{\sqrt{\mathcal{D}_{1}/\sqrt[4]{r}}}}.

(18)

From Eqs. (II.1) and (9), for the PDF along the fingers we find

\displaystyle\hat{p}_{2}(k_{y},s|0)

\displaystyle=\frac{s^{-1}(s+r)^{1/2}}{(s+r)^{1/2}+\mathcal{D}_{2}\,k_{y}^{2}},

(19)

where $\mathcal{D}_{2}=\frac{\mathcal{D}_{y}}{2\sqrt{\mathcal{D}_{z}}}$ . From here, by applying the inverse Fourier transform we obtain

\hat{p}_{2}(y,s|0)=\frac{(s+r)^{1/4}}{2s\sqrt{\mathcal{D}_{2}}}e^{-\frac{(s+r)^{1/4}}{\sqrt{\mathcal{D}_{2}}}|y|}.

(20)

The marginal PDF $p_{2}(y,t|0)$ provides the transport equation along the main finger and it is governed by the equation

\displaystyle\int_{0}^{t}\zeta(t-t^{\prime})\frac{\partial}{\partial t^{\prime}}p_{2}(y,t^{\prime}|0)\,dt^{\prime}=\mathcal{D}_{2}\frac{\partial^{2}}{\partial y^{2}}p_{2}(y,t|0),

(21)

with the initial condition $p_{2}(y,t=0|0)=\delta(y)$ . Here the kernel $\zeta(t)$ is

\displaystyle\zeta(t)=t^{-1/2}E_{1,1/2}^{-1/2}(-rt)=\frac{1}{\sqrt{\pi t}}e^{-rt}+\sqrt{r}\,\mathrm{erf}\left(\sqrt{rt}\right).

(22)

Eq. (21) can also be presented by means of the regularized Prabhakar derivative (17). It reads

\displaystyle{{}_{C}}\mathbf{D}_{1,-r,0+}^{1/2,1/2}p_{2}(y,t|0)=\mathcal{D}_{2}\frac{\partial^{2}}{\partial y^{2}}p_{2}(y,t|0),

(23)

or equivalently

	$\displaystyle{{}_{TC}}D_{r}^{1/2}p_{2}(y,t\|0)=\mathcal{D}_{2}\frac{\partial^{2}}{\partial y^{2}}p_{2}(y,t\|0)$
	$\displaystyle-\sqrt{r}\int_{0}^{t}\mathrm{erf}\left(\sqrt{r(t-t^{\prime})}\right)\frac{\partial}{\partial t^{\prime}}p_{2}(y,t^{\prime}\|0)\,dt^{\prime},$		(24)

where $\mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}\,dt$ is the error function, while

\displaystyle{{}_{TC}}D_{b}^{\alpha}f(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}e^{-b(t-t^{\prime})}(t-t^{\prime})^{-\alpha}\frac{d}{dt^{\prime}}f(t^{\prime})\,dt^{\prime}

(25)

is the tempered Caputo derivative with the exponential truncation, where $b>0$ is the truncation parameter fcaa2015 ; sandev mathematics . For the stationary PDF along the $y$ direction, we find

p_{2,\mathrm{st}}(y|0)=\frac{1}{\sqrt{4\mathcal{D}_{2}/\sqrt{r}}}e^{-\frac{|y|}{\sqrt{\mathcal{D}_{2}/\sqrt{r}}}}.

(26)

For the $z$ direction, we have

\displaystyle\hat{p}_{3}(k_{z},s|0)=\frac{s^{-1}(s+r)}{(s+r)+\mathcal{D}_{z}k_{z}^{2}},

(27)

that yields

\displaystyle\hat{p}_{3}(z,s|0)=\frac{1}{2\sqrt{\mathcal{D}_{z}}}s^{-1}(s+r)^{1/2}e^{-\frac{(s+r)^{1/2}}{\sqrt{\mathcal{D}_{3}}}|z|},

(28)

where $\mathcal{D}_{3}=\mathcal{D}_{z}$ . The corresponding equation for the transport along secondary fingers reads

\displaystyle\int_{0}^{t}\xi(t-t^{\prime})\frac{\partial}{\partial t^{\prime}}p_{3}(z,t^{\prime}|0)\,dt^{\prime}=\mathcal{D}_{3}\frac{\partial^{2}}{\partial z^{2}}p_{3}(z,t|0).

(29)

where

\displaystyle\xi(t)=\delta(t)+r,

(30)

and it can be rewritten in the equivalent form

\displaystyle\frac{\partial}{\partial t}p_{3}(z,t|0)=\mathcal{D}_{3}\frac{\partial^{2}}{\partial z^{2}}p_{3}(z,t|0)-rp_{3}(z,t|0),

(31)

which is a diffusion-absorption equation. The stationary PDF along the $z$ direction is

\displaystyle p_{3,\mathrm{st}}(z|0)=\frac{1}{\sqrt{4\mathcal{D}_{3}/r}}e^{-\frac{|z|}{\sqrt{\mathcal{D}_{3}/r}}}.

(32)

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 $\mathcal{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 $\sqrt{r}$ for the $y$ component and on $\sqrt[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.4.

II.3 Mean squared displacements

In the section, we analyze the MSDs along all three directions,

\left\langle x^{2}(t)\right\rangle=\int_{-\infty}^{\infty}x^{2}\,p_{1}(x,t|x_{0})\,dx,

\left\langle y^{2}(t)\right\rangle=\int_{-\infty}^{\infty}y^{2}\,p_{2}(y,t|0)\,dy,

\left\langle z^{2}(t)\right\rangle=\int_{-\infty}^{\infty}z^{2}\,p_{3}(z,t|0)\,dz.

Taking into account corresponding solutions for the marginal PDFs, we have

	$\displaystyle\left\langle x^{2}(t)\right\rangle$	$\displaystyle=x_{0}^{2}+2\,\mathcal{D}_{1}\,t^{1/4}E_{1,5/4}^{1/4}(-rt)$
		$\displaystyle=x_{0}^{2}+2\,\mathcal{D}_{1}\,t^{1/4}\left[\frac{1}{\sqrt[4]{rt}}-\frac{\textrm{E}_{3/4}(rt)}{\Gamma(1/4)}\right],$		(33)

where $\textrm{E}_{n}(z)=\int_{1}^{\infty}\frac{e^{-zt}}{t^{n}}\,dt$ is the exponential integral function. This corresponds to the transition from subdiffusion to localization

\left\langle x^{2}(t)\right\rangle\sim x_{0}^{2}+2\,\mathcal{D}_{1}\left\{\begin{array}[]{l l}&\frac{t^{1/4}}{\Gamma(5/4)},\quad rt\ll 1\\ \\ &\frac{1}{\sqrt[4]{r}},\quad rt\gg 1.\end{array}\right.

For the $y$ fingers we have

\displaystyle\left\langle y^{2}(t)\right\rangle=2\,\mathcal{D}_{2}\,\frac{\mathrm{erf}(\sqrt{rt})}{\sqrt{r}},

(34)

that corresponds to the transition from subdiffusion to localization as well,

\left\langle y^{2}(t)\right\rangle\sim 2\,\mathcal{D}_{2}\left\{\begin{array}[]{l l}&\frac{t^{1/2}}{\Gamma(3/2)},\quad rt\ll 1,\\ \\ &\frac{1}{\sqrt{r}},\quad rt\gg 1,\end{array}\right.

Eventually, the MSD for $z$ - fingers reads

\displaystyle\left\langle z^{2}(t)\right\rangle=2\,\mathcal{D}_{3}\,\frac{1-e^{-rt}}{r},

(35)

that corresponds to saturation in the long time limit,

\left\langle z^{2}(t)\right\rangle\sim 2\,\mathcal{D}_{3}\left\{\begin{array}[]{l l}&t,\quad rt\ll 1\\ \\ &\frac{1}{r},\quad rt\gg 1.\end{array}\right.

Therefore, unlike an initial transient behaviour all the MSDs saturate towards a constant value (exhibiting stochastic localization) as in the scenario of one-dimensional 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 axel . This variety of cases has been also obtained from stochastic simulations of the process based on a Langevin equation approach, see Section II.4.

II.4 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 iomin jstat ; lenzi njp , and where resets, as a renewal process, can be easily performed, see Ref. pal_pre2015 . The system of coupled Langevin equations reads


	$\displaystyle x(t+\Delta t)=\left\{\begin{array}[]{ll}x(0),\,\textrm{with prob.}\>r\,\Delta t,\\ \\ x(t)+\beta_{1}A(y)B(z)\eta_{x}(t),\,\textrm{with prob.}\>(1-r\,\Delta t),\end{array}\right.$		(36d)

	$\displaystyle y(t+\Delta t)=\left\{\begin{array}[]{ll}y(0),\,\textrm{with prob.}\>r\,\Delta t,\\ \\ y(t)+\beta_{2}B(z)\eta_{y}(t),\,\textrm{with prob.}\>(1-r\,\Delta t),\end{array}\right.$		(36h)

	$\displaystyle z(t+\Delta t)=\left\{\begin{array}[]{ll}z(0),\,\textrm{with prob.}\>r\,\Delta t,\\ \\ z(t)+\beta_{3}\eta_{z}(t),\,\textrm{with prob.}\>(1-r\,\Delta t),\end{array}\right.$		(36l)

where $\beta_{1}$ , $\beta_{2}$ , $\beta_{3}$ are constants related to the diffusion coefficients $\mathcal{D}_{1}$ , $\mathcal{D}_{2}$ , $\mathcal{D}_{3}$ , $\eta_{x}(t)$ , $\eta_{y}(t)$ , $\eta_{z}(t)$ are zero mean Gaussian noises ( $\langle\eta_{x}(t)\rangle=0,\>\langle\eta_{y}(t)\rangle=0,\>\langle\eta_{z}(t)\rangle=0$ ), $A(y)$ and $B(z)$ are functions introduced to mimic $\delta$ -functions (see Refs. iomin jstat ; jstat2020 ), and $r$ is the parameter of the Poisson process. To replicate the Dirac $\delta$ -function, diffusion across the $x$ and $y$ directions is permitted in a narrow band of thickness $2\varepsilon$ along the $x$ and $y$ axes. As a result, the noise in Eqs. (36d) and (36h) is multiplicative, however in Refs. lenzi njp ; jstat2020 the authors verified that the value $\varepsilon$ has no influence in the diffusive process, as long as $\varepsilon$ and the noise amplitudes $\beta_{1}$ , $\beta_{2}$ , $\beta_{3}$ are of the same order of magnitude. In our simulations, we have set $\varepsilon=\beta_{x}=\beta_{y}=\beta_{z}=0.1$ . The noises $\eta_{x}(t)$ , $\eta_{y}(t)$ , $\eta_{z}(t)$ , were sampled from a Gaussian distribution $N(0,\Delta 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 three-dimensional 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\times 10^{4}$ particle positions were simulated considering a time step of $\Delta 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 $\sigma^{2}_{x}(t)=\langle(x(t)-\langle x(t)\rangle)^{2}\rangle$ , $\sigma^{2}_{y}(t)=\langle(y(t)-\langle y(t)\rangle)^{2}\rangle$ , $\sigma^{2}_{z}(t)=\langle(z(t)-\langle z(t)\rangle)^{2}\rangle$ . In Fig. 3 we give comparison of the analytical and simulation results for the marginal PDFs, where we use that $\beta_{i}=\sqrt{2\,\mathcal{D}_{i}}$ , $i=\{1,2,3\}$ , with $\Delta t=1$ , see iomin jstat ; pal_pre2015 . 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\,\mathcal{D}_{i}=\beta_{i}^{2}$ . For more results on the corresponding PDFs obtained by the numerical simulations we refer to the Supplemental material 1.

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. (II.1)). 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

		$\displaystyle\frac{\partial}{\partial t}P(x,y,z,t\|x_{0},0,0)=L_{\mathrm{FP}}P(x,y,z,t\|x_{0},0,0)$
		$\displaystyle+r\delta(y)\delta(z)\int_{-\infty}^{\infty}dy^{\prime}\int_{-\infty}^{\infty}dz^{\prime}P(x,y^{\prime},z^{\prime},t\|x_{0},0,0),$

which differs from Eq. (II.1) 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 $\delta(x-x_{0})\delta(y)\delta(z)$ term in Eq. (II.1). 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. (81). From the Fourier-Laplace transformations, we arrive at the following PDF in the Fourier-Laplace space, see Supplemental material 2 for details,

	$\displaystyle\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)$
	$\displaystyle=\frac{1}{s}\times\frac{s\,(s+r)^{-3/4}}{s\,(s+r)^{-3/4}+\mathcal{D}_{1}\,k_{x}^{2}}\,e^{ik_{x}x_{0}}$
	$\displaystyle\times\frac{(s+r)^{1/2}}{(s+r)^{1/2}+\mathcal{D}_{2}\,k_{y}^{2}}\times\frac{s+r}{s+r+\mathcal{D}_{3}\,k_{z}^{2}}.$		(38)

From this equation, the marginal PDFs for the three axis can be straightforwardly obtained as done in the previous section:

\hat{p}_{1}(k_{x},s)=\frac{(s+r)^{-3/4}}{s\,(s+r)^{-3/4}+\mathcal{D}_{1}\,k_{x}^{2}}\,e^{ik_{x}x_{0}},

(39)

\hat{p}_{2}(k_{y},s)=\frac{s^{-1}\,(s+r)^{1/2}}{(s+r)^{1/2}+\mathcal{D}_{2}\,k_{y}^{2}},

(40)

\hat{p}_{3}(k_{z},s)=\frac{s^{-1}\,(s+r)}{(s+r)+\mathcal{D}_{3}\,k_{z}^{2}}.

(41)

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

\hat{p}_{1}(x,s)=\frac{1}{2\sqrt{\mathcal{D}_{1}}}s^{-1/2}(s+r)^{-3/8}\,e^{-\frac{s^{1/2}(s+r)^{-3/8}}{\sqrt{\mathcal{D}_{1}}}|x-x_{0}|},

(42)

which is governed by the equation

\displaystyle{{}_{TC}}D_{r}^{1/4}p_{1}(x,t|x_{0})=\mathcal{D}_{1}\frac{\partial^{2}}{\partial x^{2}}p_{1}(x,t|x_{0}),

(43)

where ${{}_{TC}}D_{b}^{\alpha}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

\langle x^{2}(t)\rangle=x_{0}^{2}+2\,\mathcal{D}_{1}\,t^{1/4}E_{1,5/4}^{-3/4}(-r\,t).

(44)

Its asymptotics reads

\displaystyle\left\langle x^{2}(t)\right\rangle\sim x_{0}^{2}+2\,\mathcal{D}_{1}\left\{\begin{array}[]{l l}&\frac{t^{1/4}}{\Gamma(5/4)},\quad rt\ll 1,\\ \\ &r^{3/4}\,t,\quad rt\gg 1.\end{array}\right.

(48)

. 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. (100) that in the long time limit, $\langle x^{2}(t)\rangle\sim t$ , i.e. it scales diffusively. The short time limit scales as $\langle x^{2}(t)\rangle\sim 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 $\varepsilon$ , $\beta_{i}$ and $\mathcal{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

	$\displaystyle\frac{\partial}{\partial t}P(x,y,z,t\|x_{0},0,0)=L_{\mathrm{FP}}P(x,y,z,t\|x_{0},0,0)$
	$\displaystyle-rP(x,y,z,t\|x_{0},0,0)+r\delta(z)\int_{-\infty}^{\infty}dz^{\prime}P(x,y,z^{\prime},t\|x_{0},0,0),$		(49)

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,

	$\displaystyle\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)$
	$\displaystyle=\frac{1}{s}\times\frac{s^{1/2}\,(s+r)^{-1/4}}{s^{1/2}\,(s+r)^{-1/4}+\mathcal{D}_{1}\,k_{x}^{2}}\,e^{ik_{x}x_{0}}$
	$\displaystyle\times\frac{s\,(s+r)^{-1/2}}{s\,(s+r)^{-1/2}+\mathcal{D}_{2}\,k_{y}^{2}}\times\frac{s+r}{s+r+\mathcal{D}_{3}\,k_{z}^{2}}.$		(50)

This also yields the images of the marginal PDFs for the different axis. Eventually, we have

\hat{p}_{1}(k_{x},s)=\frac{s^{-1/2}\,(s+r)^{-1/4}}{s^{1/2}\,(s+r)^{-1/4}+\mathcal{D}_{1}\,k_{x}^{2}}\,e^{ik_{x}x_{0}},

(51)

\hat{p}_{2}(k_{y},s)=\frac{(s+r)^{-1/2}}{s\,(s+r)^{-1/2}+\mathcal{D}_{2}\,k_{y}^{2}},

(52)

\hat{p}_{3}(k_{z},s)=\frac{s^{-1}\,(s+r)}{(s+r)+\mathcal{D}_{3}\,k_{z}^{2}}.

(53)

From these expressions one obtains the corresponding MSDs, which are

\langle x^{2}(t)\rangle=x_{0}^{2}+2\,\mathcal{D}_{1}\,t^{1/4}E_{1,5/4}^{-1/4}(-r\,t),

(54)

	$\displaystyle\langle y^{2}(t)\rangle$	$\displaystyle=2\,\mathcal{D}_{2}\,t^{1/2}E_{1,3/2}^{-1/2}(-r\,t)$
		$\displaystyle=2\,\mathcal{D}_{2}\left[e^{-rt}\frac{t^{1/2}}{\Gamma(1/2)}+\frac{2rt+1}{2\sqrt{r}}\mathrm{erf}(\sqrt{rt})\right],$		(55)

and the MSD along the $z$ axis is the same as in the previous cases. Their asymptotics read

\displaystyle\left\langle x^{2}(t)\right\rangle\sim x_{0}^{2}+2\,\mathcal{D}_{1}\left\{\begin{array}[]{l l}&\frac{t^{1/4}}{\Gamma(5/4)},\quad rt\ll 1,\\ \\ &\frac{r^{1/4}\,t^{1/2}}{\Gamma(3/2)},\quad rt\gg 1,\end{array}\right.

(59)

\displaystyle\left\langle y^{2}(t)\right\rangle\sim 2\,\mathcal{D}_{2}\left\{\begin{array}[]{l l}&\frac{t^{1/2}}{\Gamma(3/2)},\quad rt\ll 1,\\ \\ &r^{1/2}\,t,\quad rt\gg 1,\end{array}\right.

(63)

In this case, the MSD in the $x$ -axis behaves subdiffusively with $\langle x^{2}(t)\rangle\sim t^{1/4}$ as in the case with no resetting, and then it turns to $\langle x^{2}(t)\rangle\sim t^{1/2}$ , which means an accelerating subdiffusive transport. Along the $y$ -axis, the MSD scales as $\langle y^{2}(t)\rangle\sim t^{1/2}$ in the short time limit, and then it turns to linear dependence in time $\langle y^{2}(t)\rangle\sim 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 $\varepsilon$ , $\beta_{i}$ and $\mathcal{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 two-dimensional $xy$ -comb. The $y$ and $z$ axes in the three-dimensional 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 $\delta(y)$ and $\delta(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 non-zero flux along the $x$ backbone and $y$ fingers, and for the two dimensional case it was discussed in Refs. dvork2009 ; LZ98 ; IoZaPf16 ; SaIoMe16 . Here, we extend the arguments of Refs. IoZaPf16 ; SaIoMe16 for the three dimensional case of Eq. (1). Let us consider the Liouville equation

\frac{\partial}{\partial t}P+\mathrm{div}\,\mathbf{j}=0,

(64)

where the three dimensional current $\mathbf{j}=(j_{x},\,j_{y},\,j_{z})$ describes Markov processes in Eq. (1). In this case, the three-dimensional current reads


	$\displaystyle j_{x}=-D_{x}(y,z)\,\frac{\partial}{\partial x}P(x,y,z,t),$		(65a)
	$\displaystyle j_{y}=-D_{y}(z)\,\frac{\partial}{\partial y}P(x,y,z,t),$		(65b)
	$\displaystyle j_{z}=-D_{z}\,\frac{\partial}{\partial z}P(x,y,z,t).$		(65c)

Here, we take a general diffusivity function in the $x$ and $y$ directions $D_{x}(y,z)$ and $D_{y}(z)$ , respectively (instead of $\mathcal{D}_{x}\,\delta(y)\delta(z)$ and $\mathcal{D}_{y}\,\delta(z)$ in Eq. (1)). Therefore, Eq. (64) together with Eqs. (65), can be regarded as the three-dimensional non-Markovian master equation.

Integrating Eq. (64) with respect to $y$ and $z$ from $-\epsilon/2$ to $\epsilon/2$ : $\int_{-\epsilon/2}^{\epsilon/2}dy\dots$ and $\int_{-\epsilon/2}^{\epsilon/2}dz\dots$ , one obtains for the l.h.s. of the equation, after application of the middle point theorem, $\epsilon^{2}\,\frac{\partial}{\partial t}P(x,y=0,z=0,t)$ , which is exact in the limit $\epsilon\rightarrow 0$ . This term can be neglected in this limit $\epsilon\rightarrow 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 $\epsilon$ 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. (65c) we obtain that this term responsible for the transport in the $z$ direction reads

\mathcal{D}_{z}\Big{[}P^{\prime}(x,y,z,t)\big{|}_{z=\tfrac{\epsilon}{2}}-\\ P^{\prime}(x,y,z,t)\big{|}_{z=-\tfrac{\epsilon}{2}}\Big{]}\,,

where prime means derivative with respect to $z$ . This corresponds to the two outgoing fluxes from the $y$ fingers in the $\pm z$ directions: $F_{z}(+)+F_{z}(-)$ . The transport in the $y$ direction, after integration, is

\epsilon\,D_{y}(z\rightarrow 0)\,\partial_{y}^{2}P(x,y,z=0,t)\equiv F_{y}\,.

It should be stressed that the second derivative over $y$ , presented in the form

\epsilon\,\frac{\partial^{2}}{\partial y^{2}}P=\left[\frac{\partial}{\partial y}P(y+\epsilon/2)-\frac{\partial}{\partial y}P(y-\epsilon/2)\right]\,,

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}\neq 0$ , as containing an incoming flux. Therefore, $\epsilon\,D_{y}(z\rightarrow 0)\neq 0$ . Taking $D_{y}(z)=\frac{1}{\pi}\frac{\epsilon\mathcal{D}_{y}}{y^{2}+\epsilon^{2}}$ , one obtains in the limit $\epsilon\rightarrow 0$ a nonzero flux $F_{y}$ with $D_{y}(z)=\mathcal{D}_{y}\,\delta(z)$ , which is the diffusion coefficient in the $y$ direction in Eq. (64).

Now we perform integration in the $\epsilon$ vicinity of the $x$ backbone, where we take into account the singularity of the $y$ component of diffusion coefficient, which is $\mathcal{D}_{y}\,\delta(z)$ . We also admit that integration of the $j_{z}$ current in Eq. (65c) with respect to $y$ yields zero. Therefore, integration with respect to $y$ and $z$ yields from Eq. (65b) the term responsible for the transport in the $y$ direction as follows

\mathcal{D}_{y}\Big{[}P^{\prime}(x,y,z=0,t)\big{|}_{y=\tfrac{\epsilon}{2}}-P^{\prime}(x,y,z=0,t)\big{|}_{y=-\tfrac{\epsilon}{2}}\Big{]}\,.

Here prime means derivative with respect to $y$ . This corresponds to the two outgoing fluxes from the backbone in the $\pm y$ directions: $F_{y}(+)+F_{y}(-)$ . The transport along the $x$ direction, after integration of Eq. (65a), is

\epsilon^{2}\,D(y\rightarrow 0,z\rightarrow 0)\,\frac{\partial^{2}}{\partial x^{2}}P(x,y=0,z=0,t)=\\ F_{x}(x+\epsilon)+F_{x}(x-\epsilon)\,.

In complete analogy with the $y$ coordinate, the second derivative with respect to $x$ , presented in the form

\epsilon\,\frac{\partial^{2}}{\partial x^{2}}P=\left[\frac{\partial}{\partial x}P(x+\epsilon/2)-\frac{\partial}{\partial x}P(x-\epsilon/2)\right]

as $\epsilon\rightarrow 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}\neq 0$ and correspondingly the flux $F_{x}\equiv F_{x}(+)+F_{x}(-)$ has to be nonzero as well: $F_{x}(\pm)\neq 0$ . Therefore, $\epsilon^{2}\,D(y\rightarrow 0,z\rightarrow 0)\neq 0$ . Now, taking the diffusion coefficient in the form $D(y,z)=\frac{1}{\pi}\frac{\epsilon\mathcal{D}}{y^{2}+\epsilon^{2}}\cdot\frac{1}{\pi}\frac{\epsilon\mathcal{D}}{z^{2}+\epsilon^{2}}$ , one obtains in the limit $\epsilon\rightarrow 0$ a nonzero flux $F_{x}$ with $D(y,z)=\mathcal{D}_{x}\delta(y)\delta(z)$ , which is the diffusion coefficient in the $x$ direction in Eqs. (1), (64) and (65a).

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 maj1 ; axel . 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 $\left\langle x^{2}(t)\right\rangle\sim 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: $\langle x^{2}(t)\rangle\sim t^{1/2}$ and $\langle y^{2}(t)\rangle\sim 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.

Acknowledgment

This research was partially supported by Grant No. CGL2016-78156- C2-2-R (VM and AM). TS was supported by the Alexander von Humboldt Foundation.

References

(1) A. Iomin, V. Méndez and W. Horsthemke, Fractional Dynamics in Comb-like Structures (World Scientific, 2018).
(2) V. Méndez and A. Iomin, Chaos, Solitons & Fractals 53, 46 (2013).
(3) A. Iomin and V. Méndez, Chaos, Solitons & Fractals 82, 142 (2016).
(4) M. R. Evans and S. N. Majumdar, Phys. Rev. Lett. 106, 160601 (2011).
(5) M. R. Evans and S. N. Majumdar, J. Phys. A: Math. Theor. 44, 435001 (2011).
(6) M. R. Evans and S. N. Majumdar, J. Phys. A: Math. Theor. 47, 285001 (2014).
(7) S. Gupta, S. N. Majumdar, and G. Schehr, Phys. Rev. Lett. 112, 220601 (2014).
(8) X. Durang, M. Henkel, and H. Parl, J. Phys. A: Math. Theor. 47, 045002 (2014).
(9) A. Pal, Phys. Rev. E 91, 012113 (2015).
(10) C. Christou and A. Schadschneider, J. Phys. A: Math. Theor. 48, 285003 (2015).
(11) J. Masoliver, Phys. Rev. E 99, 012121 (2019).
(12) D. Gupta, J. Phys. A: Math. Theor. 47, 033212 (2019).
(13) A. A. Tateishi, H. V. Ribeiro, and E. K. Lenzi, Front. Phys. 5, 52 (2017).
(14) M. Montero and J. Villarroel, Phys. Rev. E 87, 012116 (2013).
(15) V. Méndez and D. Campos, Phys. Rev. E 93, 022106 (2016).
(16) A. Pal, A. Kundu, and M. R. Evans, J. Phys. A: Math. Theor. 49, 225001 (2016).
(17) A. Nagar and S. Gupta, Phys. Rev. E 93, 060102 (2016).
(18) S. Eule and J. J. Metzger, New J. Phys. 18, 033006 (2016).
(19) M. Montero, A. Masó-Puigdellosas and J. Villarroel, Euro. Phys. J. B 90, 176 (2017).
(20) V. P. Shkilev, Phys. Rev. E 96, 012126 (2017).
(21) M. R. Evans and S. N. Majumdar, J. Phys. A: Math. Theor. 52, 01LT01 (2018).
(22) A. Masó-Puigdellosas, D. Campos, and V. Méndez, Phys. Rev. E 99, 012141 (2019).
(23) A. Masó-Puigdellosas, D. Campos, and V. Méndez, J. Stat. Mech.: Theor. Exp. 2019, 033201.
(24) A. S. Bodrova, A. V. Chechkin, and I. M. Sokolov, Phys. Rev. E 100, 012119 (2019).
(25) A. S. Bodrova, A. V. Chechkin, and I. M. Sokolov, Phys. Rev. E 100, 012120 (2019).
(26) A. Masó-Puigdellosas, D. Campos, and V. Méndez, Front. Phys. 2019 00112.
(27) A. Masó-Puigdellosas, D. Campos, and V. Méndez, Phys. Rev. E 100, 042104 (2019).
(28) D. Campos, and V. Méndez, Phys. Rev. E 92, 062115 (2015).
(29) S. Reuveni, Phys. Rev. Lett. 116 170601 (2016).
(30) A. Pal and S. Reuveni, Phys. Rev. Lett. 118, 030603 (2017).
(31) A. Chechkin and I. M. Sokolov, Phys. Rev. Lett. 121, 050601 (2018).
(32) U. Basu, A. Kundu, and A. Pal, Phys. Rev. E 100, 032136 (2019); A. Pal and V. V. Prasad, Phys. Rev. E 99, 032123 (2019).
(33) J. Masoliver and M. Montero, Phys. Rev. E 100, 042103 (2019).
(34) A. Pal and V. V. Prasad, Phys. Rev. Research 1, 032001 (2019).
(35) I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products (Academic Press, San Diego, 2007).
(36) A. M. Mathai, R. K. Saxena and H. J. Haubold, The $H$ -function: Theory and Applications ( Springer, 2010).
(37) T. Sandev, A. Chechkin, H. Kantz, R. Metzler, Fract. Calc. Appl. Anal. 18, 1006 (2015).
(38) R. Garra and R. Garrappa, Commun. Nonlinear Sci. Numer. Simul. 56, 314 (2018).
(39) T. Sandev, A. V. Chechkin, N. Korabel, H. Kantz, I. M. Sokolov, and R. Metzler, Phys. Rev. E 92, 042117 (2015).
(40) T. R. Prabhakar, Yokohama Math. J. 19, 7 (1971).
(41) M. D’Ovidio and F. Polito, Theory Probab. Appl. 62, 552 (2018); arXiv:1307.1696 (2013); R. Garra, R. Gorenflo, F. Polito, and Z. Tomovski, Appl. Math. Comput. 242 576 (2014).
(42) T. Sandev, Mathematics 5, 66 (2017); T. Sandev, W. Deng, and P. Xu, J. Phys. A: Math. Theor. 51, 405002 (2018).
(43) T. Sandev and A. Iomin, Europhys. Lett. 124, 20005 (2018).
(44) A. Stanislavsky and A. Weron, Phys. Rev. Research 1, 023006 (2019); A Stanislavsky, A Weron, J. Chem. Phys. 149, 044107 (2018).
(45) V. Méndez, A. Iomin, W. Horsthemke and D. Campos, J. Stat. Mech. 2017, 063205 (2017).
(46) H. V. Ribeiro, A. A. Tateishi, L. G. A. Alves, R. S. Zola and E. K. Lenzi, New J. Phys. 16, 093050 (2014).
(47) A. Pal, Phys. Rev. E 91, 012113 (2015).
(48) E. K. Lenzi, T. Sandev, H. V. Ribeiro, P. Jovanovski, A. Iomin, and L. Kocarev, J. Stat. Mech. doi: 10.1088/1742-5468/ab7af4 (2020).
(49) O.A. Dvoretskaya, P.S. Kondratenko, and L.V. Matveev, JETP 110, 58 (2010).
(50) I.A. Lubashevskii and A.A. Zemlyanov, JETP 87, 700 (1998).
(51) A. Iomin, V. Zaburdaev and T. Pfohl, Chaos, Solitons & Fractals 92, 115 (2016).
(52) T. Sandev, A. Iomin, and V. Méndez, J. Phys. A: Math. Theor. 49, 355001 (2016).

Supplemental material: Stochastic resetting on comb-like structures, V. Domazetoski, A. Masó-Puigdellosas, T. Sandev, V. Méndez, A. Iomin and L. Kocarev

.1 Solutions for diffusion with global resetting

The diffusion equation in a three-dimensional comb with global resetting can be represented by

	$\displaystyle\frac{\partial}{\partial t}P(x,y,z,t\|x_{0},0,0)$	$\displaystyle=\mathcal{D}_{x}\delta(y)\delta(z)\frac{\partial^{2}}{\partial x^{2}}P(x,y,z,t\|x_{0},0,0)+\mathcal{D}_{y}\delta(z)\frac{\partial^{2}}{\partial y^{2}}P(x,y,z,t\|x_{0},0,0)$
		$\displaystyle+\mathcal{D}_{z}\frac{\partial^{2}}{\partial z^{2}}P(x,y,z,t\|x_{0},0,0)-rP(x,y,z,t\|x_{0},0,0)+r\delta(x-x_{0})\delta(y)\delta(z).$		(66)

Here we consider the initial condition $P(x,y,z,t=0|x_{0})=\delta(x-x_{0})\delta(y)\delta(z)$ . The Laplace transformation of Eq. (.1) reads

	$\displaystyle s\hat{P}(x,y,z,s\|x_{0},0,0)-\delta(x-x_{0})\delta(y)\delta(z)$	$\displaystyle=\mathcal{D}_{x}\delta(y)\delta(z)\frac{\partial^{2}}{\partial x^{2}}\hat{P}(x,y,z,s\|x_{0},0,0)+\mathcal{D}_{y}\delta(z)\frac{\partial^{2}}{\partial y^{2}}\hat{P}(x,y,z,s\|x_{0},0,0)$
		$\displaystyle+\mathcal{D}_{z}\frac{\partial^{2}}{\partial z^{2}}\hat{P}(x,y,z,s\|x_{0},0,0)-r\hat{P}(x,y,z,s\|x_{0},0,0)+\frac{r}{s}\delta(x-x_{0})\delta(y)\delta(z).$		(67)

Then, by Fourier transform with respect to $x$ , $y$ and $z$ we obtain

	$\displaystyle s\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)-e^{\imath k_{x}x_{0}}$	$\displaystyle=-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s\|x_{0},0,0)-\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s\|x_{0},0,0)$
		$\displaystyle-\mathcal{D}_{z}k_{z}^{2}\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)-r\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)+\frac{r}{s}e^{\imath k_{x}x_{0}}.$		(68)

from where it follows

\displaystyle\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)=\frac{s^{-1}(s+r)e^{\imath k_{x}x_{0}}}{(s+r+\mathcal{D}_{z}k_{z}^{2})}-\frac{\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)}{(s+r+\mathcal{D}_{z}k_{z}^{2})}-\frac{\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)}{(s+r+\mathcal{D}_{z}k_{z}^{2})},

(69)

The inverse Fourier transform with respect to $k_{z}$ yields

	$\displaystyle\hat{P}(k_{x},k_{y},z,s\|x_{0},0,0)$	$\displaystyle=\frac{(s+r)^{-1/2}}{2\sqrt{\mathcal{D}_{z}}}\left[s^{-1}(s+r)e^{\imath k_{x}x_{0}}-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s\|x_{0},0,0)\right.$
		$\displaystyle\left.-\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s\|x_{0},0,0)\right]e^{-\sqrt{\frac{s+r}{\mathcal{D}_{z}}}\|z\|},$		(70)

from where we find $\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)$ ,

\displaystyle\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)

\displaystyle=\frac{1}{\mathcal{D}_{y}}\frac{s^{-1}(s+r)e^{\imath k_{x}x_{0}}-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)}{\frac{2\sqrt{\mathcal{D}_{z}}}{\mathcal{D}_{y}}(s+r)^{1/2}+k_{y}^{2}}.

(71)

By applying the inverse Fourier transform with respect to $k_{y}$ , we have

\displaystyle\hat{P}(k_{x},y,z=0,s|x_{0},0,0)

\displaystyle=\frac{(s+r)^{-1/4}}{2\sqrt{2\mathcal{D}_{y}\sqrt{D}_{z}}}\left[s^{-1}(s+r)e^{\imath k_{x}x_{0}}-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)\right]e^{-\sqrt{\frac{2\sqrt{\mathcal{D}_{z}}}{\mathcal{D}_{y}}}},

(72)

and thus

\displaystyle\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)

\displaystyle=\frac{1}{\mathcal{D}_{x}}\frac{s^{-1}(s+r)e^{\imath k_{x}x_{0}}}{\frac{2\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}{\mathcal{D}_{x}}(s+r)^{1/4}+k_{x}^{2}}.

(73)

Therefore, we finally obtain

\displaystyle\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)=s^{-1}\frac{(s+r)^{1/4}}{(s+r)^{1/4}+\frac{\mathcal{D}_{x}}{2\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}k_{x}^{2}}\times\frac{(s+r)^{1/2}}{(s+r)^{1/2}+\frac{\mathcal{D}_{y}}{2\sqrt{\mathcal{D}_{z}}}k_{y}^{2}}\times\frac{(s+r)}{(s+r)+\mathcal{D}_{z}k_{z}^{2}}\times e^{\imath k_{x}x_{0}}.

(74)

After we find the marginal PDFs we will analyze the mean squared displacements along all three directions, $\left\langle x^{2}(t)\right\rangle=\int_{-\infty}^{\infty}x^{2}\,p_{1}(x,t|x_{0})\,dx$ , $\left\langle y^{2}(t)\right\rangle=\int_{-\infty}^{\infty}y^{2}\,p_{2}(y,t|0)\,dy$ , and $\left\langle y^{2}(t)\right\rangle=\int_{-\infty}^{\infty}z^{2}\,p_{3}(z,t|0)\,dz$ . Therefore, we have

	$\displaystyle\left\langle x^{2}(t)\right\rangle$	$\displaystyle=2\left(\frac{\mathcal{D}_{x}}{2\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}\right)\mathcal{L}^{-1}\left[\frac{s^{-1}}{(s+r)^{1/4}}\right]=2\left(\frac{\mathcal{D}_{x}}{2\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}\right)t^{1/4}E_{1,5/4}^{1/4}(-rt)$
		$\displaystyle=2\left(\frac{\mathcal{D}_{x}}{2\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}\right)t^{1/4}\left[\frac{1}{\sqrt[4]{rt}}-\frac{\textrm{E}_{3/4}(rt)}{\Gamma(1/4)}\right],$		(75)

where

E_{\alpha,\beta}^{\delta}(z)=\sum_{k=0}^{\infty}\frac{(\delta)_{k}}{\Gamma(\alpha k+\beta)}\frac{z^{k}}{k!}

(76)

is the three parameter Mittag-Leffler function, and $(\delta)_{k}=\Gamma(\delta+k)/\Gamma(\delta)$ is the Pochhammer symbol, $\textrm{E}_{n}(z)=\int_{1}^{\infty}\frac{e^{-zt}}{t^{n}}\,dt$ is the exponential integral function,

\displaystyle\left\langle y^{2}(t)\right\rangle=2\left(\frac{\mathcal{D}_{y}}{2\sqrt{\mathcal{D}_{z}}}\right)\mathcal{L}^{-1}\left[\frac{s^{-1}}{(s+r)^{1/2}}\right]=2\left(\frac{\mathcal{D}_{y}}{2\sqrt{\mathcal{D}_{z}}}\right)t^{1/2}E_{1,3/2}^{1/2}(-rt)=2\left(\frac{\mathcal{D}_{y}}{2\sqrt{\mathcal{D}_{z}}}\right)\frac{\mathrm{erf}(\sqrt{rt})}{\sqrt{r}},

(77)

where $\mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}\,dt$ gives the error function, and

\displaystyle\left\langle z^{2}(t)\right\rangle=2\,\mathcal{D}_{z}\,\mathcal{L}^{-1}\left[\frac{s^{-1}}{s+r}\right]=\frac{2\,\mathcal{D}_{z}}{r}\left[1-e^{-rt}\right].

(78)

The Laplace transform of the three parameter Mittag-Leffler function (76) reads

\mathcal{L}\left[t^{\beta-1}E_{\alpha,\beta}^{\delta}(\pm at^{\alpha})\right](s)=\frac{s^{\alpha\delta-\beta}}{\left(s^{\alpha}\mp a\right)^{\delta}},\quad\Re(s)>|a|^{1/\alpha},

(79)

and its asymptotic behaviors are given by

E_{\alpha,\beta}^{\gamma}(-z^{\alpha})\simeq\left\{\begin{array}[]{ll}\frac{1}{\Gamma(\beta)}\exp\left(-\gamma\frac{\Gamma(\beta)}{\Gamma(\alpha+\beta)}z^{\alpha}\right),\quad z\ll 1.\\ \frac{z^{-\alpha\gamma}}{\Gamma(\beta-\alpha\gamma)},\quad z\gg 1.\end{array}\right.

(80)

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.4, with different values of the ensemble size of $2\times 10^{4}$ and the time span of $5^{4}$ .

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

\begin{split}\frac{\partial}{\partial t}P(x,y,z,t|x_{0},0,0)&=\mathcal{D}_{x}\delta(y)\delta(z)\frac{\partial^{2}}{\partial x^{2}}P(x,y,z,t|x_{0},0,0)+\mathcal{D}_{y}\delta(z)\frac{\partial^{2}}{\partial y^{2}}P(x,y,z,t|x_{0},0,0)+\mathcal{D}_{z}\frac{\partial^{2}}{\partial z^{2}}P(x,y,z,t|x_{0},0,0)\\ &-rP(x,y,z,t|x_{0},0,0)+r\delta(y)\delta(z)\int_{-\infty}^{\infty}dy^{\prime}\int_{-\infty}^{\infty}dz^{\prime}P(x,y^{\prime},z^{\prime},t|x_{0},0,0).\end{split}

(81)

Here, instead of appearing at the particular position $(x_{0},0,0)$ as stated by $\delta(x-x_{0})\delta(y)\delta(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. (81)). 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

\begin{split}s\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)-e^{ik_{x}x_{0}}&=-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)-\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)\\ &-\mathcal{D}_{z}k_{z}^{2}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)-r\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)\\ &+r\hat{P}(k_{x},k_{y}=0,k_{z}=0,s|x_{0},0,0).\end{split}

(82)

Now, isolating the propagator one gets

\begin{split}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)&=\frac{e^{ik_{x}x_{0}}}{s+r+\mathcal{D}_{z}k_{z}^{2}}-\frac{\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)}{s+r+\mathcal{D}_{z}k_{z}^{2}}-\frac{\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)}{s+r+\mathcal{D}_{z}k_{z}^{2}}\\ &+\frac{r\hat{P}(k_{x},k_{y}=0,k_{z}=0,s|x_{0},0,0)}{s+r+\mathcal{D}_{z}k_{z}^{2}}.\end{split}

(83)

Let us proceed to find the three unknown terms on the right hand side. Taking $k_{y}=k_{z}=0$ one can isolate

\hat{P}(k_{x},k_{y}=0,k_{z}=0,s|x_{0},0,0)=\frac{e^{ik_{x}x_{0}}}{s}-\frac{\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)}{s},

(84)

with which the equation for the propagator reads

\begin{split}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)&=\frac{s+r}{s+r+\mathcal{D}_{z}k_{z}^{2}}\frac{e^{ik_{x}x_{0}}}{s}-\frac{s+r}{s+r+\mathcal{D}_{z}k_{z}^{2}}\frac{\mathcal{D}_{x}k_{x}^{2}}{s}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)\\ &-\frac{\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)}{s+r+\mathcal{D}_{z}k_{z}^{2}}.\end{split}

(85)

The other two terms can be obtained by taking the inverse Fourier transform of $k_{y}$ and $k_{z}$ and evaluating at $y=0$ and $z=0$ respectively. Starting by the last, one gets that

\begin{split}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)&=\frac{s+r}{2\sqrt{D_{z}(s+r)}+\mathcal{D}_{y}k_{y}^{2}}\frac{1}{s}\left(e^{ik_{x}x_{0}}-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)\right).\end{split}

(86)

Proceeding similarly for the $y$ variable one gets

\begin{split}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)&=\frac{e^{ik_{x}x_{0}}}{2\frac{s}{s+r}\sqrt{2D_{y}\sqrt{D_{z}(s+r)}}+\mathcal{D}_{x}k_{x}^{2}}.\end{split}

(87)

Putting all the elements into Eq. (85) and simplifying some terms one gets the following explicit expression for the propagator:

\begin{split}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)&=\frac{e^{ik_{x}x_{0}}4\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}(s+r)}}\sqrt{\mathcal{D}_{z}(s+r)}}{\left(s+r+\mathcal{D}_{z}k_{z}^{2}\right)\left(2\sqrt{\mathcal{D}_{z}(s+r)}+\mathcal{D}_{y}k_{y}^{2}\right)\left(2\frac{s}{s+r}\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}(s+r)}}+\mathcal{D}_{x}k_{x}^{2}\right)}.\end{split}

(88)

Finally, from these expressions one can get the corresponding MSD,

\langle x^{2}(t)\rangle=x_{0}^{2}+\frac{\mathcal{D}_{x}}{\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}\,\mathcal{L}^{-1}\left[\frac{(s+r)^{3/4}}{s^{2}}\right]=x_{0}^{2}+\frac{\mathcal{D}_{x}}{\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}\,t^{1/4}\,E_{1,5/4}^{-3/4}(-r\,t),

(89)

\langle y^{2}(t)\rangle=\frac{\mathcal{D}_{y}}{\sqrt{\mathcal{D}_{z}}}\,\mathcal{L}^{-1}\left[\frac{1}{s(s+r)^{1/2}}\right]=\frac{\mathcal{D}_{y}}{\sqrt{\mathcal{D}_{z}}}\frac{\mathrm{erf}(\sqrt{rt})}{\sqrt{r}}

(90)

\langle z^{2}(t)\rangle=2\,\mathcal{D}_{z}\,\mathcal{L}^{-1}\left[\frac{1}{s(s+r)}\right]=2\,\mathcal{D}_{z}\,\frac{1-e^{-rt}}{r}.

(91)

In Fig. 13, we give graphical representations of the PDFs of the process under resetting to the backbone.

.3 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 this case, the governing equation reads

\begin{split}\frac{\partial}{\partial t}P(x,y,z,t|x_{0},0,0)&=\mathcal{D}_{x}\delta(y)\delta(z)\frac{\partial^{2}}{\partial x^{2}}P(x,y,z,t|x_{0},0,0)+\mathcal{D}_{y}\delta(z)\frac{\partial^{2}}{\partial y^{2}}P(x,y,z,t|x_{0},0,0)+\mathcal{D}_{z}\frac{\partial^{2}}{\partial z^{2}}P(x,y,z,t|x_{0},0,0)\\ &-rP(x,y,z,t|x_{0},0,0)+r\delta(z)\int_{-\infty}^{\infty}dz^{\prime}P(x,y^{\prime},z^{\prime},t|x_{0},0,0),\end{split}

(92)

where the last term is now the marginal distribution in the variables $x$ and $y$ instead of only the $x$ variable. The corresponding equation in the Fourier-Laplace space reads

\begin{split}s\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)-e^{ik_{x}x_{0}}&=-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)-\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)\\ &-\mathcal{D}_{z}k_{z}^{2}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)-r\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)\\ &+r\hat{P}(k_{x},k_{y},k_{z}=0,s|x_{0},0,0),\end{split}

(93)

from which the propagator can be isolated to give

\begin{split}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)&=\frac{e^{ik_{x}x_{0}}}{s+r+\mathcal{D}_{z}k_{z}^{2}}-\frac{\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)}{s+r+\mathcal{D}_{z}k_{z}^{2}}-\frac{\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)}{s+r+\mathcal{D}_{z}k_{z}^{2}}\\ &+\frac{r\hat{P}(k_{x},k_{y},k_{z}=0,s|x_{0},0,0)}{s+r+\mathcal{D}_{z}k_{z}^{2}}.\end{split}

(94)

and, by taking $k_{z}=0$ , we can isolate the following expression for the marginal PDF:

\begin{split}\hat{P}(k_{x},k_{y},k_{z}=0,s|x_{0},0,0)&=\frac{e^{ik_{x}x_{0}}}{s}-\frac{\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)}{s}-\frac{\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)}{s}.\end{split}

(95)

Introducing this expression to Eq. (94) we get

\begin{split}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)&=\frac{s+r}{s(s+r+\mathcal{D}_{z}\,k_{z}^{2})}\left(e^{ik_{x}x_{0}}-\mathcal{D}_{x}\,k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)-\mathcal{D}_{y}\,k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)\right).\end{split}

(96)

As done for the previous cases, the two unknown terms remaining on the right hand side can be found by applying the inverse Fourier transform with respect to the $y$ and $z$ variables to Eq. (96) and evaluating at $y=z=0$ . Starting by the last, one finds that

\begin{split}\hat{P}(k_{x},k_{y},z=0,s|x_{0},0,0)&=\frac{1}{2\frac{s}{s+r}\sqrt{\mathcal{D}_{z}(s+r)}+\mathcal{D}_{y}k_{y}^{2}}\left(e^{ik_{x}x_{0}}-\mathcal{D}_{x}\,k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)\right)\end{split}

(97)

and proceeding equally for the $y$ variable:

\begin{split}\hat{P}(k_{x},y=0,z=0,s|x_{0},0,0)&=\frac{e^{ik_{x}x_{0}}}{2\sqrt{2\mathcal{D}_{y}\frac{s}{s+r}\sqrt{\mathcal{D}_{z}(s+r)}}+\mathcal{D}_{x}\,k_{x}^{2}}.\end{split}

(98)

Putting all the elements together one finds a formal expression for the propagator:

\begin{split}\hat{P}(k_{x},k_{y},k_{z},s|x_{0},0,0)&=\frac{4e^{ik_{x}x_{0}}\sqrt{\mathcal{D}_{z}(s+r)}\sqrt{2\mathcal{D}_{y}\frac{s}{s+r}\sqrt{\mathcal{D}_{z}(s+r)}}}{(s+r+\mathcal{D}_{z}\,k_{z}^{2})\left(2\frac{s}{s+r}\sqrt{\mathcal{D}_{z}(s+r)}+\mathcal{D}_{y}\,k_{y}^{2}\right)\left(2\sqrt{2\mathcal{D}_{y}\frac{s}{s+r}\sqrt{\mathcal{D}_{z}(s+r)}}+\mathcal{D}_{x}\,k_{x}^{2}\right)}.\end{split}

(99)

Finally, the corresponding MSD in the Laplace space read

\langle x^{2}(t)\rangle=x_{0}^{2}+\frac{\mathcal{D}_{x}}{\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}\,\mathcal{L}^{-1}\left[\frac{(s+r)^{1/4}}{s^{3/2}}\right]=x_{0}^{2}+\frac{\mathcal{D}_{x}}{\sqrt{2\mathcal{D}_{y}\sqrt{\mathcal{D}_{z}}}}\,t^{1/4}\,E_{1,5/4}^{-1/4}(-r\,t),

(100)

\langle y^{2}(t)\rangle=\frac{\mathcal{D}_{y}}{\sqrt{\mathcal{D}_{z}}}\,\mathcal{L}^{-1}\left[\frac{(s+r)^{1/2}}{s^{2}}\right]=\frac{\mathcal{D}_{y}}{\sqrt{\mathcal{D}_{z}}}\,t^{1/2}\,E_{1,3/2}^{-1/2}(-r\,t)=\frac{\mathcal{D}_{y}}{\sqrt{\mathcal{D}_{z}}}\left[e^{-rt}\frac{t^{1/2}}{\Gamma(1/2)}+\frac{2rt+1}{2\sqrt{r}}\,\mathrm{erf}(\sqrt{rt})\right],

(101)

\langle z^{2}(t)\rangle=2\,\mathcal{D}_{z}\,\mathcal{L}^{-1}\left[\frac{1}{s(s+r)}\right]=\frac{2\,\mathcal{D}_{z}}{r}\left[1-e^{-rt}\right].

(102)

In Fig. 14, we give graphical representations of the PDFs of the process under resetting to the backbone.

	$\displaystyle\frac{\partial}{\partial t}P(x,y,z,t\|x_{0},0,0)$	$\displaystyle=\mathcal{D}_{x}\delta(y)\delta(z)\frac{\partial^{2}}{\partial x^{2}}P(x,y,z,t\|x_{0},0,0)+\mathcal{D}_{y}\delta(z)\frac{\partial^{2}}{\partial y^{2}}P(x,y,z,t\|x_{0},0,0)$
		$\displaystyle+\mathcal{D}_{z}\frac{\partial^{2}}{\partial z^{2}}P(x,y,z,t\|x_{0},0,0)-rP(x,y,z,t\|x_{0},0,0)+r\delta(x-x_{0})\delta(y)\delta(z).$		(66)

	$\displaystyle s\hat{P}(x,y,z,s\|x_{0},0,0)-\delta(x-x_{0})\delta(y)\delta(z)$	$\displaystyle=\mathcal{D}_{x}\delta(y)\delta(z)\frac{\partial^{2}}{\partial x^{2}}\hat{P}(x,y,z,s\|x_{0},0,0)+\mathcal{D}_{y}\delta(z)\frac{\partial^{2}}{\partial y^{2}}\hat{P}(x,y,z,s\|x_{0},0,0)$
		$\displaystyle+\mathcal{D}_{z}\frac{\partial^{2}}{\partial z^{2}}\hat{P}(x,y,z,s\|x_{0},0,0)-r\hat{P}(x,y,z,s\|x_{0},0,0)+\frac{r}{s}\delta(x-x_{0})\delta(y)\delta(z).$		(67)

	$\displaystyle s\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)-e^{\imath k_{x}x_{0}}$	$\displaystyle=-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s\|x_{0},0,0)-\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s\|x_{0},0,0)$
		$\displaystyle-\mathcal{D}_{z}k_{z}^{2}\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)-r\hat{P}(k_{x},k_{y},k_{z},s\|x_{0},0,0)+\frac{r}{s}e^{\imath k_{x}x_{0}}.$		(68)

	$\displaystyle\hat{P}(k_{x},k_{y},z,s\|x_{0},0,0)$	$\displaystyle=\frac{(s+r)^{-1/2}}{2\sqrt{\mathcal{D}_{z}}}\left[s^{-1}(s+r)e^{\imath k_{x}x_{0}}-\mathcal{D}_{x}k_{x}^{2}\hat{P}(k_{x},y=0,z=0,s\|x_{0},0,0)\right.$
		$\displaystyle\left.-\mathcal{D}_{y}k_{y}^{2}\hat{P}(k_{x},k_{y},z=0,s\|x_{0},0,0)\right]e^{-\sqrt{\frac{s+r}{\mathcal{D}_{z}}}\|z\|},$		(70)

Stochastic resetting on comb-like structures***V. D. and A. M. contributed equally to this work.

Abstract

pacs:

I Introduction

II Global resetting

II.1 Analytical results

II.2 Marginal PDFs

II.3 Mean squared displacements

II.4 Langevin equation approach. Numerical simulations

III Resetting to the backbone

IV Resetting to the main fingers

V Remarks on two dimensional comb

VI Remark on the three dimensional comb geometry

VII Summary

Acknowledgment

References

Supplemental material: Stochastic resetting on comb-like structures, V. Domazetoski, A. Masó-Puigdellosas, T. Sandev, V. Méndez, A. Iomin and L. Kocarev

.1 Solutions for diffusion with global resetting

.2 Solutions for diffusion with reseting to the backbone

.3 Solutions for diffusion with reseting to the main fingers

Stochastic resetting on comb-like structures^*^**V. D. and A. M. contributed equally to this work.