Power-law relaxation of a confined diffusing particle subject to resetting with memory

We study the relaxation of a Brownian particle with long range memory under confinement in one dimension. The particle diffuses in an arbitrary confining potential and resets at random times to previously visited positions, chosen with a probability proportional to the local time spent there by the particle since the initial time. This model mimics an animal which moves erratically in its home range and returns preferentially to familiar places from time to time, as observed in nature. The steady state density of the position is given by the equilibrium Gibbs-Boltzmann distribution, as in standard diffusion, while the transient part of the density can be obtained through a mapping of the Fokker-Planck equation of the process to a Schr\"odinger eigenvalue problem. Due to memory, the approach at late times toward the steady state is critically self-organised, in the sense that it always follows a sluggish power-law form, in contrast to the exponential decay that characterises Markov processes. The exponent of this power-law depends in a simple way on the resetting rate and on the leading relaxation rate of the Brownian particle in the absence of resetting. We apply these findings to several exactly solvable examples, such as the harmonic, V-shaped and box potentials.


Introduction
Let us consider a single Brownian particle which, in addition to diffusion, undergoes a resetting process with rate r to previously visited positions.Resetting follows a stochastic rule that incorporates memory effects: the particle that resets at some time t chooses uniformly at random a preceding time 0 ≤ t ′ ≤ t, i.e., with probability density 1/t, and occupies the position at which it was located at t ′ again.Thus the resetting move tends to bring the particle near the most visited positions in the past.This model was introduced in the context of random walks on lattices and solved exactly in Ref. [1].The model was further extended to a free Brownian particle on the line [2].The position distribution at time t was found to approach a Gaussian form at late times with the variance growing anomalously slowly as (2D/r) ln t for large t, where D is the diffusion constant of the particle.This slow dynamics emerges from the resetting induced memory effect: the particle is sluggish to move away from its familiar territory, where the most visited sites are located.Nevertheless, the position distribution is always time-dependent and does not approach a stationary state at late times, contrary to diffusive processes subject to resetting to a single point, which generically admit non-equilibrium steady states (NESSs) [3,4,5].
Various other generalisations of this simple model have been studied in the recent past, for instance by considering a decaying memory [6], Lévy flights [7] or active particles [8].A central limit theorem and an anomalous large deviation principle have also been established for a class of memory walks of this type, for a broad range of memory kernels [9,10].The rigorous proofs of these latter results are based on a connection between the resetting process to previous sites and the growth of weighted random recursive trees.
In a more applied context, random walks with long range memory have become increasingly useful in ecology for the description and analysis of animal mobility.There is mounting evidence that animals do not follow pure Markov processes but use their memory and tend to revisit preferred places during ranging [11,12,13,14,15,16].The model above was able to describe quantitatively the movement patterns of Capuchin monkeys in the wild [1], as well as of individual elks released in an unknown environment [17].Another feature shared by many animal trajectories is spatial localisation, by opposition to unbounded diffusion, as home ranges are often clearly identifiable.The range-resident behaviour of animals can be modelled by incorporating in the simple random walk a central place (such as a den) toward which the walker is attracted, producing an effective |x| potential and therefore a stationary 'equilibrium' state at large times [18].Recent studies have considered harmonic potentials and Ornstein-Uhlenbeck (OU) processes as the starting point of models of home range movements of increasing complexity [19,20,21].
In this paper, we study a Brownian particle in one dimension under the combined effects of confinement by an external potential U (x) and of memory, as described in the rule above.In this situation, we expect the position distribution P r (x, t) to approach a stationary form at late times, where the subscript r in P r (x, t) denotes the nonzero resetting rate.We first recall that in the absence of resetting, i.e., for r = 0, the system approaches the equilibrium stationary state characterised by the Gibbs-Boltzmann distribution where the friction coefficient of the particle is set to 1 from now on and the partition function Z = ∞ −∞ e −U (x)/D dx normalises the distribution to unity.We assume that U (x) is sufficiently confining so that the integral Z is finite.Moreover, the position distribution relaxes generically to the equilibrium measure in Eq. ( 1) exponentially fast.More precisely, writing P 0 (x, t) = P st 0 (x) + P tr 0 (x, t) where the superscript 'tr' denotes the time-dependent transient part of the distribution, one finds that for large t where λ 1 denotes the first nonzero eigenvalue of the diffusion operator in the confining potential U (x).The precise value of λ 1 depends on the potential U (x).The x dependence of the leading transient is implicit in the amplitude on the right hand side (rhs) of Eq. ( 2).In fact, the decay exponent λ 1 , up to a multiplicative constant, is precisely the gap between the first excited state E 1 and the ground state E 0 = 0 of the associated Schrödinger operator [22] (see Section II).
Here, we investigate the behaviour of the position distribution P r (x, t), when one switches memory on with a non-zero rate r.In the presence of a finite resetting rate r > 0, we again expect the position distribution to approach a stationary form P st r (x) = P r (x, t → ∞) at late times.Our first goal is to characterise this stationary position distribution.Moreover, it is also important to study how the system relaxes to this stationary state for nonzero r, i.e., how the result in Eq. ( 2) gets modified with a nonzero r.Our main results are twofold: (i) We show that the stationary position distribution P st r (x) for nonzero r remains of the Gibbs-Boltzmann form independently of r, i.e., Like in the r = 0 case, this stationary state for r > 0 is actually an equilibrium state.This property markedly differs from the well-studied problem of a Brownian particle in a potential U (x) and subject to resetting to a fixed position [23,24,25,26,27,28,29,30,31], a case that manifestly violates detailed balance and generates a NESS.
(ii) Even more interestingly, we find that in complete contrast to the exponential decay in Eq. ( 2) for r = 0, the relaxation to the stationary state for r > 0 is always algebraic, or "critical".More precisely, the analogue of Eq. ( 2) now reads for large t where the dependence of the amplitude on x and r is implicit.The exponent θ 1 depends on r and also on the confining potential U (x).In fact we show that the exponent θ 1 can be expressed in terms of the decay rate λ 1 of the resetting-free case via the exact relation Thus as r increases, the exponent θ 1 decreases and the relaxation dynamics becomes slower.
The rest of the paper is organised as follows.In Section 2, to set the stage and our notations, we recall well-known results and the method to study diffusion in a confining potential in the absence of resetting, i.e., for r = 0.In Section 3, we derive our main results for the stationary state and the transient part in the presence of a finite resetting rate r > 0. In Section 4, after recovering previous results for the unbounded case U (x) = 0, we exemplify the calculations for the harmonic and |x| potentials.Section 5 is devoted to the case of a particle in a box, x ∈ [0, L] with reflecting boundary conditions at the two boundaries, where we derive the full probability distribution P r (x, t) at all times.A comparison of the theoretical results with numerical simulations is given in Section 6 and we conclude in Section 7.

Diffusion in a confining potential without resetting
In this section, we recall the basic results and the general method to compute the time-dependent position distribution of a Brownian particle diffusing in the presence of a confining potential U (x) [22].The position x(t) of the particle evolves via the over-damped Langevin equation where η(t) is a Gaussian white noise with zero mean ⟨η(t)⟩ = 0 and a delta correlator ⟨η(t)η(t ′ )⟩ = 2 D δ(t − t ′ ), with D the diffusion constant.(We recall that the friction coefficient is set to unity.)The associated Fokker-Planck equation reads starting from the initial condition P 0 (x, 0) = δ(x − x 0 ) and satisfying the boundary conditions: P 0 (x → ±∞, t) = 0.
To solve this partial differential equation, one can use the method of separation of variables and use the ansatz Substituting ( 8) in (7) and dividing both sides by ϕ(x) f (t) one gets ḟ (t) where λ is a constant independent of x and t.Solving for f (t) gives f (t) = A 0 e −λt where A 0 is a constant and one must have λ ≥ 0 so that the solution does not diverge as t → ∞.The function ϕ(x) satisfies a second order differential equation where λ plays the role of an eigenvalue, yet to be determined.One can further reduce this eigenvalue equation to a more familiar Schrödinger form via the substitution It is then easy to see that ψ(x) satisfies the time-independent Schrödinger equation (setting m = ℏ = 1) where the quantum potential V Q (x) is expressed in terms of the classical potential U (x) via the relation Thus one needs to solve the Schrödinger equation ( 12), find its spectrum, i.e., the eigenvalues E n = λ n /(2D) and the associated eigenfunctions ψ λn (x).Note that in general, the spectrum of V Q (x) may consist of both bound and scattering states, even though U (x) may be fully confining.Once we have the full spectrum, one can then write down the most general solution of Eq. ( 7) as a linear combination of these eigenfunctions where the coefficients {a λn } can be determined from the initial condition.Note that the quantum potential V Q (x) is such that its ground state corresponds to zero energy λ 0 /2D = E 0 = 0, with eigenfunction ψ 0 (x) ∝ e −U (x)/2D as one can easily check.Thus separating the ground state from the rest of the spectrum in Eq. ( 14), we get One then immediately identifies the first term on the rhs of ( 15) as the stationary Gibbs distribution The second term in (15), summing over all nonzero modes, represents the transient part of the distribution P tr 0 (x, t).If the first excited state λ 1 > 0 is strictly positive, i.e., well separated from the ground state λ 0 = 0, then it governs the leading late time decay of the transient and one obtains Note that the result in Eq. ( 17) holds provided that a λ 1 ̸ = 0. Sometimes, the initial condition may have a special symmetry that renders a λ 1 = 0.In that case, the leading decay will be governed by the next nonzero excited state.
Let us remark that the result in Eq. ( 17) is valid for any confining potential U (x), such that (i) the integral Z = ∞ −∞ e −U (x)/D dx exists, (ii) the associated quantum potential V Q (x) in Eq. ( 13) is such that its first excited state has eigenvalue λ 1 > 0 (recall that the ground state has zero energy if a steady state exists).In addition, the leading decay in Eq. ( 17) will be a pure exponential if λ 1 < λ 2 , i.e., there is a nonzero gap between the first and the second excited state.There are some potentials U (x) for which λ 1 > 0 may be gapped from the ground state, but then there is a continuum of scattering states beyond λ 1 .In this case, depending on how the density of scattering states vanishes as one approaches λ 1 from above, one may have an additional sub-leading algebraic correction to the leading exponential decay ∼ e −λ 1 t .We discuss below two examples that can be explicitly solved to illustrate these features.
Ornstein-Uhlenbeck (OU) process.As a first simple example, consider the OU process with stiffness µ > 0, for which U (x) = µx 2 /2.In this case, the quantum potential in Eq. ( 13) reads Thus, one has a quantum harmonic oscillator of mass m = 1 and frequency ω = µ/(2D), with the height of the potential shifted by µ/(4D).Consequently, the energy levels are well separated from each other and are given by The eigenvalues are given by λ n = 2DE n = n µ and the leading transient in Eq. ( 17) decays as where ψ 1 (x) represents the eigenfunction associated with the first excited state of the harmonic oscillator.Note that if the initial condition is symmetric, e.g., when P 0 (x) = δ(x), this symmetry is preserved at all times, indicating that only the even eigenfunctions contribute to the spectral decomposition in Eq. ( 14).This means that all the odd coefficients vanish: a λ 2m+1 = 0 for all m = 0, 1, 2, . ... In that case, the leading transient corresponds to the second excited state of the harmonic oscillator and P tr 0 (x, t) decays as ∼ e −2µt at late times.
The potential U (x) = F |x| with F > 0. This is an interesting example since, in this case, the quantum potential in Eq. ( 13) reads This quantum potential has a single bound state at E = 0 and a continuum of excited states with energies and we expect that the transient in Eq. ( 17) will decay, to leading order for large t, exponentially as e −λ 1 t ∼ e −F 2 t/(4D) .However, unlike in the harmonic oscillator example above, the eigenvalue λ 1 is not isolated here and there is a continuum of scattering states above E 1 .Moreover, the density of such states vanishes as √ E − E 1 as E → E 1 from above (see Section 4) and the leading decay of the transient in Eq. ( 17) now behaves as e −F 2 t/(4D) t −3/2 .The sub-leading t −3/2 decay emerges from replacing the sum in Eq. ( 15) by an integral above E 1 .Since the density of states vanishes as √ E − E 1 , this leads to the power law correction t −3/2 multiplying the leading exponential e −F 2 t/(4D) [32].

Diffusion in a confining potential in the presence of resetting
We now allow the particle to perform resetting moves, in addition to the diffusion in the confining external potential.Let us recall the memory-driven resetting dynamics.At time t, with rate r the particle chooses any previous time in the past 0 ≤ t ′ < t with uniform probability density 1/t and resets to the position that it occupied at that previous instant t ′ .Let P r (x, t) denote the position distribution at time t.To take into account the memory effect, we now need to define the two-point function P r (x, t ; x ′ , t ′ ), which is the joint probability density for the particle to be near x ′ at t ′ and near x at t, with t ′ ≤ t.Clearly, if we integrate over one of the positions, we recover the marginal one-point probability density With these ingredients at hand, we can now write down the Fokker-Planck equation for the evolution of the one-point position distribution P r (x, t) as The first two terms on the rhs describe, as before, the diffusion in the external potential U (x).The third term describes the loss of probability density from position x at time t due to resetting to other positions.The last term on the rhs describes the gain in the probability density at x at time t due to resetting from other positions labelled by x ′ .If the particle has to reset to x from x ′ at time t, it must have been at x at a previous time t ′ ≤ t and the probability density of this event is simply the two-point function P r (x ′ , t; x, t ′ ).Finally, we need to integrate over all positions x ′ from which the particle may arrive at x by resetting.The reason for the solvability for the one-point position distribution in this model can then traced back to the second relation in Eq. (22), which allows us to write a closed equation for the one-point function (without involving multiple-point functions) provided that the resetting rate r is independent of the position.We obtain It is easy to check that the total probability ∞ −∞ P r (x, t) dx remains conserved and equals unity due to the initial condition, e.g., P r (x, 0) = δ(x − x 0 ).One also needs to impose the vanishing boundary conditions as x → ±∞, i.e., P r (x → ±∞, t) = 0 for all t.Evidently, for r = 0, Eq. ( 24) reduces to the standard Fokker-Planck equation (7).
To solve the Fokker-Planck equation ( 24), we follow the same route as in the r = 0 case, namely, we use the method of separation of variables with the ansatz Substituting (25) in Eq. ( 24), dividing both sides by ϕ(x)f r (t) and assembling the only-time dependent and only-space dependent parts separately gives ḟr (t) where λ ≥ 0 is a constant independent of x and t.Note that the space-dependent function ϕ(x) is independent of r and satisfies the same eigenvalue equation (10) as the r = 0 case.Consequently, the analysis of the spatial part is exactly as in the r = 0 case discussed in the previous section.In particular, we get where ψ(x) satisfies the Schrödinger equation ( 12), with the same eigenvalues λ n 's as in the r = 0 case.Thus, the r-dependence enters only in the time-dependent part f r (t) of the solution.For a given λ, we get from Eq. ( 26) a second order ordinary differential equation for f r (t) Multiplying this equation by t and taking the time derivative, one obtains Through the change of variable z = −(λ + r)t with g(z) ≡ f (t), Eq. ( 29) can be recast as a confluent hypergeometric equation, with b = 1 and a = λ/(λ+r).The differential equation ( 30) has two linearly independent solutions denoted by M (a, b, z) and U (a, b, z) [33].Hence, the most general solution of Eq. ( 29) can be expressed as where c 1 and c 2 are two arbitrary constants.The solution must be finite at all t, in particular at t = 0.The small argument asymptotics of the two solutions are given by [33] M (a, b, z) → 1 as z → 0 (32) Consequently, the divergence as t → 0 in Eq. ( 31) is avoided by setting c 2 = 0, and our solution simply reads where we have set c 1 = 1 by absorbing this constant into the amplitude of the eigenmode.Importantly, using the identity M (0, b, z) = 1, we recover f r (t) = 1 for λ = 0.This means that the solution ϕ(x) with λ = 0 in the eigenvalue equation (26), i.e., the Gibbs-Boltzmann distribution ϕ(x) = a 0 e − 1 D U (x) , is also a stationary state when r > 0. To obtain the asymptotic behaviour of f r (t) in time for λ > 0, it is useful to recall here some basic properties of the confluent hypergeometric function M (a, b, z) ≡ 1 F 1 (a, b, z).It has a simple power series expansion whereas for large negative arguments, we have [33] M Consequently, f r (t) in Eq. ( 34) behaves for large t as a power-law, Note that this result does not hold for r = 0, as the prefactor in Eq. (37) vanishes in this case.Rather, the exponential decay f r (t) = e −λt is recovered by setting r = 0 directly in Eq. ( 34), since Finally, combining the time-dependent part f r (t) and the space-dependent part that involves the set of all eigenvalues of the Schrödinger operator, we can express the full solution of P r (x, t) as the following exact spectral decomposition As in the r = 0 case before, we will now separate out the stationary mode corresponding to the ground state n = 0 with λ 0 = 0 from the transient part and rewrite Eq. (38) as At late times, using the asymptotic behaviour of f r (t) in Eq. (37), it follows that the n-th mode in the second term on the rhs decays at late times as a power law ∼ t − λn r+λn for λ n > 0. Consequently, the sum in Eq. (39) represents the transient part P tr r (x, t), while the first term that is independent of t represents the stationary solution.
Hence, we come to the two main results announced in the Introduction.First, the stationary solution is independent of r and remains with the Gibbs-Boltzmann form as in the r = 0 case.This solution corresponds to an equilibrium state, since local detailed balance is fulfilled even with r > 0. This can be understood by considering two small disjoint regions [x, x + δ x ] and [y, y + δ y ], which have been occupied until time t for an total amount of time τ x (t) and τ y (t), respectively.The probability that the particle takes a move from the region x to y by resetting at time t is δ x P st (x)w (r) x→y , where the transition rate is w (r) x→y = rτ y (t)/t due to the preferential revisit rule.Since δ x P st (x) ≃ τ x (t)/t at large t, the above transition probability reads rτ y (t)τ x (t)/t 2 , which is symmetric with respect to x and y, implying detailed balance.This property is a specificity of the linear reinforcement rule.We note that it also follows from Eq. ( 24): if P r (x, t) becomes independent of time for all x at large t, then the 4th term in the rhs becomes rP st (x) and cancels with the 3rd term.Therefore the equation for the stationary state is the same as with r = 0.
Secondly, the leading behaviour of the transient part decays as where λ 1 = min n∈N + {λ n } represents the gap between the ground state and the first excited state of the underlying Schrödinger operator of the r = 0 case.Thus, unlike the stationary solution that remains unaffected by resetting, the transient part gets modified drastically if r > 0. It now decays very slowly as a power law at late times with a non-universal exponent θ 1 = λ 1 /(r + λ 1 ) comprised in the interval (0, 1) and that decreases with increasing r, indicating that the dynamics get slower and slower due to resetting.This result is completely general, valid for any shape of the confining potential.

Unbounded case U (x) = 0
Before proceeding to examples of confining potentials, we show how the free particle results, in particular the logarithmic growth of the mean square displacement [1,2], can be recovered from the analysis above.If the particle diffuses in an unbounded space, or U (x) = 0 for x ∈ (−∞, ∞), the Schrödinger equation ( 12) reads ψ ′′ (x) = −(λ/D)ψ(x).
The non-diverging solutions are of the form ψ(x) = cos(kx + φ) with λ = λ k = Dk 2 and k real.Since the potential does not have a finite Z, a 0 = 0 and Eq.(39) becomes a continuous sum over the modes k with a k and φ k determined from the initial condition.Eq. ( 42) is equivalent to the solution presented in [2] under a slightly different form.At late times, with the asymptotics (36) the solution simplifies to where r + λ k has been replaced by r in the last integral due to the fact that at large ln(rt), the small k regime dominates the integration.Expression (44) is the inverse Fourier transform of e −σ(t)k 2 /2 , therefore the asymptotic form of P r (x, t) is Gaussian with variance σ(t) = (2D/r) ln(rt).

V-shaped potential
For the case U (x) = F |x|, it turns out that the transient part of the density behaves as As discussed in Section 2, the first excited state λ 1 = F 2 /(4D) is well separated from the ground state E 0 = 0, but there is a continuum of excited states above E 1 = λ 1 /(2D) whose density vanishes as √ E − E 1 as E → E 1 , as shown below.This leads to additional logarithmic corrections to the leading power law decay in Eq. ( 41).
Let us now derive this result.We recall that the excited eigenfunctions with energy E > E 1 of the Schrödinger equation with the potential given by Eq. ( 21) are spatially extended and are of two types, symmetric and anti-symmetric, or even and odd in x, respectively [22]: and λ = 2DE is taken as a continuous parameter.These functions are δ-normalised, or , therefore it is convenient to use the index k and expand the transient part of the density in this basis.Eq. (39) becomes where the dispersion relation is given from Eq. ( 49) by The coefficients α k follow from the initial condition P r (x, t) = δ(x − x 0 ), where the Kummer's function M in Eq. ( 51) is simply unity.By projecting, one obtains Therefore, We next expand the integrand of Eq. ( 54) at small k, as it dominates in the large t limit.For a fixed x, we have from Eqs, (47)-(48), Inserting these expressions into Eq.(54) and using the large time expressions in Eqs. ( 36)-(37) for the Kummer's function, one obtains the leading behaviour in k of the integrand.One notices that when ln t is large, the Fourier integral is actually controlled by the small k regime, This expansion cannot be used in the limit x → ∞, as we have considered kx ≪ 1 or |x| ≪ (ln t) 1/2 .Performing the Gaussian integral and using Eq.(50) leads to the final expression which is of the form of Eq. ( 46) at large ln t.From Eq. (58), we see that the number of states in [k, k + dk] is k 2 dk at small wave-numbers, which, from Eq. ( 52), is equivalent to a density of state going as

Exact solution for a particle in an interval with reflecting boundary conditions
In this section, we present the exact time-dependent solution for the position distribution of the particle with memory diffusing in a box potential, i.e., an interval [0, L] with reflecting boundary conditions.This is the case of a singular potential U (x) in Eq. ( 24).The general mapping to the quantum problem will still work, with V Q (x) representing a box potential which is zero in [0, L] and infinite outside.However, we can avoid using the quantum mapping and solve directly the Fokker-Planck equation instead.In this case, the Fokker-Planck equation ( 24) reduces to for 0 ≤ x ≤ L. The initial condition is P r (x, 0) = δ(x − x 0 ) and P r (x, t) satisfies the reflecting boundary conditions at x = 0 and x = L, We decompose P r (x, t) into its Fourier modes and write The reflecting boundary condition at x = 0 is automatically satisfied since we have only kept the cosine modes.The boundary condition at x = L imposes sin(kL) = 0, implying that k takes the values k n = nπ/L where n = 0, 1, 2, . ... Substituting (62) into the linear equation (60), one finds, as expected, that the different modes decouple and Pr (k, t) satisfies the evolution equation for any This equation is the same as Eq. ( 28) with λ = Dk 2 .Hence, we can directly write down the full exact solution for Pr (k, t) by replacing λ by Dk 2 in Eq. (34).Substituting these solutions back into the Fourier decomposition in Eq. (62) gives The coefficients a n 's can be determined from the initial condition P r (x, 0) = δ(x − x 0 ).At t = 0, the M -functions on the rhs in Eq. ( 64) are exactly 1.Hence, The orthogonality property of the cosine modes, where δ m,n is the Kronecker delta function, and m ̸ = 0, n ̸ = 0, gives from Eq. ( 65) while a 0 = 1 L .The full exact solution thus reads The first term 1/L represents the equilibrium steady state P st r (x) where the particle gets uniformly distributed over the full interval, independently of the resetting rate r.The second term, the sum over the excited modes, represents the transient part P tr r (x, t).To leading order for large t, the dominant contribution comes from the mode n = 1 and again decays algebraically as This result is perfectly consistent with the general result in Eq. (41) since λ 1 = Dπ 2 /L 2 indeed represents the gap between the first excited state and the ground state of a quantum particle in a box [0, L].

Numerical simulations
To verify and illustrate these theoretical results in the different cases, we have performed numerical simulations by incorporating resetting with preferential memory in standard Brownian dynamics algorithms [34].At each discrete time t = n∆t, where n a positive integer and ∆t the simulation time-step (∆t ≪ 1), with probability r∆t the particle is reset and its next position x n+1 is given by where n ′ is a random integer uniformly distributed in {0, 1, 2, . . ., n}.With the complementary probability 1 − r∆t, the position evolves through the discrete Langevin equation with ξ n+1 a Gaussian random variable of zero mean and unit variance.with reflecting boundary conditions.The mean square position slowly relaxes toward its asymptotic value ⟨x 2 ⟩(t → ∞) = L 2 /3 in a power-law way (symbols), while the memory-less case r = 0 follows an exponential behaviour (purple solid line).We have set the resetting rate to 2λ 1 = 2Dπ 2 /L 2 .The scaling law t −1/3 predicted by Eq. ( 69) is displayed as a guide to the eye.We have set L = 1, D = 1, x 0 = 0.25, ∆t = 10 −4 and the averages are performed over 10 4 independent trajectories.
Instead of extracting the full distribution P r (x, t) from the simulations, we have computed the simpler mean square position of the particle ⟨x 2 ⟩(t).If the position density obeys the generic form given by Eq. ( 4) and if no cancellation occurs, the second moment also relaxes toward its stationary value as Furthermore, the single power-law decay at late time in Eq. ( 72) is easier to observe if the separation between θ 1 and the exponent of the second excited mode, θ 2 , is as large as possible.In such a case, the condition t −θ 1 ≫ t −θ 2 is fulfilled from shorter times and the contribution of the sub-leading terms in Eq. ( 39) is actually negligible.For instance, for the particle confined in the box [0, L], one has θ n = D(nπ/L) 2 r+D(nπ/L) 2 from Eq. ( 68) and the gap between the first two exponent reads with λ 1 = Dπ 2 /L 2 .Notably, fixing λ 1 and varying the resetting rate from 0 to ∞, the above quantity reaches a maximum at a value r * given by t Relaxation of the particle with memory in a harmonic potential U (x) = µx 2 /2.The mean square position relaxes slowly (symbols) at finite resetting rate (here r = √ 2λ 1 = √ 2µ).The memory-less case is shown with the purple solid line.In this case, the leading power-law decay is t −θ2 with θ 2 = 2 − √ 2 (blue straight line as a guide to the eye).Here, ⟨x 2 ⟩(t → ∞) = D/µ and the stiffness is set to µ = 1, while D = 1, ∆t = 10 −3 , x 0 = 0.5 and the averages are preformed over 10 6 trajectories.
while the corresponding exponents are θ 1 (r * ) = 1/3 and θ 2 (r * ) = 2/3.Simulation results are shown in Figure 1, which represents the absolute value of the difference in Eq. (72) as a function of t, where ⟨x 2 ⟩(t → ∞) = L 2 /3 from the uniform stationary distribution.With a resetting rate set equal to the value r * in Eq. (74), one indeed observes a scaling law in excellent agreement with t −1/3 over nearly five decades, thus confirming the validity of the general expression (5).In the case r = 0, the much faster exponential decay toward the stationary state is recovered.
For the particle confined by the harmonic potential U (x) = µx 2 /2, the asymptotic second moment is ⟨x 2 ⟩(t → ∞) = D/µ, while θ n = nµ r+nµ .This case is a bit subtle since the eigenfunction of the first excited state is odd in x, given by Consequently, the contribution of ψ λ 1 (x) to ⟨x 2 ⟩(t) vanishes and the second excited state provides the leading correction for that quantity: The simulation results for r = √ 2µ, corresponding to θ 2 (r) = 2 − √ 2 = 0.5857..., are shown in Fig. 2 and exhibit a good agreement with theory.The slight mismatch at late times is attributed to a lack of precision due to the rather large ∆t.
Finally, in the case of the V-shaped potential U (x) = F |x|, Eq. (59) gives and the second moment of the equilibrium distribution is ⟨x 2 ⟩(t → ∞) = 2(D/F ) 2 .Figure 3 displays the mean square difference in a case where r = λ 1 , i.e., where the exponent of the leading power-law is 1/2 according to Eq. ( 5).As expected, the simulation curve exhibits small deviations from a t −1/2 law, in contrast to Fig. 1 of the box case that agrees perfectly with the predicted power-law.The logarithmic corrections in Eq. ( 77) are difficult to detect numerically, though; we rather measure an effective exponent of roughly 0.6 in the last two decades.

Summary and Conclusion
We have studied the relaxation of a Brownian particle with long range memory in confining potentials.The memory effects are introduced by stochastic resetting events, where the particle relocates to positions visited in the past, before it continues to diffuse in the potential from there.At large time, the particle's position tends toward a steady state given by the Gibbs-Boltzmann distribution.Contrary to confined systems subject to resetting to a fixed position [23,24,25,26,27,28,29,30,31], this stationary state is in equilibrium since local detailed balance is fulfilled.Importantly, we find that the relaxation toward the steady state is not exponential in time like in standard Brownian motion but follows an inverse power-law, for any non-zero value of the resetting rate.This sluggish decay is a consequence of the preferential character of memory in the model, where more frequently visited sites are more likely to be visited again in the future.Consequently, the particle tends to stay close to its starting point and takes a very long time to visit the accessible positions.
The exponent that controls the slow relaxation can be obtained exactly and lies in the interval (0, 1).It depends on the resetting rate and on the longest relaxation time 1/λ 1 of the memory-less particle through a simple general formula given by Eq. ( 5).In other words, knowing the asymptotic exponential relaxation of the standard Brownian particle in a potential U (x) allows us to predict the anomalous decay in the presence of memory for any value of the resetting rate.We have applied our theory to several analytically tractable potentials (box, harmonic, V-shaped), finding very good agreement with simulations.To sum up, the scale-free dynamics at large t is generic and can be compared to a self-organised critical state, independent of the initial position and characterised by a non-universal exponent.
Other models of anomalous diffusion such as continuous time random walks also exhibit a power-law relaxation toward Gibbs-Boltzmann equilibria at late times, as shown by a fractional Fokker-Planck approach [35].The mechanism for subdiffusive motion in those processes is of very different nature and is due to a broad distribution of waiting times with infinite mean.The exponent θ n of the power-law decay of the n-th eigenmode is also contained in the interval (0, 1) but turns out to be independent of n and equal to the exponent associated with the waiting time distribution, or the growth of the mean square displacement.Hence, for these systems under confinement the first excited state does not dominate the relaxation and all the modes need to be considered a priori.
Our model represents a simple description of memory-based movements of an animal in its home range, i.e., a limited region of space that can be produced by an effective confining potential [18], but where in addition to random motion, the living organism chooses from time to time to revisit familiar places.Empirically, space use by animals is known to be highly heterogeneous even within limited areas [36,1], a feature qualitatively captured by our finding that the exploration of a small accessible region can be very slow.It would be interesting to extend this model to situations where the external potential has a local minimum and to study the crossing of the particle over a barrier.

Figure 1 .
Figure1.Simulation results of the particle with memory confined in an interval [0, L] with reflecting boundary conditions.The mean square position slowly relaxes toward its asymptotic value ⟨x 2 ⟩(t → ∞) = L 2 /3 in a power-law way (symbols), while the memory-less case r = 0 follows an exponential behaviour (purple solid line).We have set the resetting rate to 2λ 1 = 2Dπ 2 /L 2 .The scaling law t −1/3 predicted by Eq. (69) is displayed as a guide to the eye.We have set L = 1, D = 1, x 0 = 0.25, ∆t = 10 −4 and the averages are performed over 10 4 independent trajectories.

1 Figure 3 .
Figure 3. Relaxation of the particle with memory in the potential U (x) = F |x|.Symbols represent the relaxation of the mean square position for a resetting rate of r = λ 1 = F 2 /(4D).The case r = 0 is represented with the purple solid line.The power-law t −1/2 is shown as a guide to the eye.Here, ⟨x 2 ⟩(t → ∞) = 2(D/F ) 2 and the force is set to F = 1, while D = 1, ∆t = 10 −3 , x 0 = 0.