Experimental Realization of Diffusion with Stochastic Resetting

Stochastic resetting is prevalent in natural and man-made systems, giving rise to a long series of nonequilibrium phenomena. Diffusion with stochastic resetting serves as a paradigmatic model to study these phenomena, but the lack of a well-controlled platform by which this process can be studied experimentally has been a major impediment to research in the field. Here, we report the experimental realization of colloidal particle diffusion and resetting via holographic optical tweezers. We provide the first experimental corroboration of central theoretical results and go on to measure the energetic cost of resetting in steady-state and first-passage scenarios. In both cases, we show that this cost cannot be made arbitrarily small because of fundamental constraints on realistic resetting protocols. The methods developed herein open the door to future experimental study of resetting phenomena beyond diffusion.

Our experimental setup is based on a home built holographic optical tweezers (HOTs) system. It uses a spatial light modulator (Hamamtsu, X10468-04) to imprint a computer generated phase pattern on an expanded laser beam (Coherent, Verdi λ = 532nm). The beam is then projected on the back aperture of a 100x objective (oil immersion, NA = 1.4) mounted on an Olympus microscope (IX71). Samples consist of a dilute colloidal suspension of spherical silica particles with a diameter of d = 1.5±0.02µm and a refractive index of n p = 1.46 (Kisker Biotech, lot# GK0611140 02) in double distilled water sealed between a glass slide and a coverslip with a sample thickness of approximately 20µm. The dense silica particle sediments to the bottom of the sample chamber and diffuses there in quasi two dimensions. The motion of the particle is recorded by a CMOS camera (Grasshopper 3, Point Gray) at a rate of 20 fps. Particle position is extracted using conventional video microscopy algorithms [S1] with an accuracy of approximately 30nm. A laser power of 1W was used to ensure sufficient trapping of the particle. We utilize in-house developed code for hardware control and data analysis.
B. Extraction of an instantaneous return trajectory from a constant radial-return-velocity trajectory.
We start by identifying the return sections (Fig. S1a, red) and wait sections (Fig. S1a,green) in the constant return velocity trajectory. We then cut out, digitally, those sections from the original trajectory. As a result we retain only the diffusive part of the trajectory, which is a trajectory with instantaneous returns (Fig. S1b). Note the reassignment of time in the extracted instantaneous return trajectory.
FIG. S1. a) A typical trajectory of a colloidal particle diffusing under a resetting process with constant radial return velocity. The return sections are marked in red and wait sections are marked in green. b) The trajectory with instantaneous returns obtained by cutting the return and wait sections of the trajectory in (a).

C. Estimation of steady state distributions
For the first section of the paper, three different sets of resetting experiments with r = 0.05s −1 were preformed and captured in 20 FPS: • Constant high velocity (7.5µm/s) returns -29 movies of 15 resetting events and 1 movie of 10 resetting events were combined to a total of 445 resetting events, all preformed on the same particle. The combined movie consists of 193958 frames, which after removing the 1 sec wait time at the end of each resetting resulted in 185058 frames.
• Constant slow velocity (0.8µm/s) returns -30 movies of 15 resetting events were combined to a total of 450 resetting events, all preformed on the same particle. The combined movie consists of 222618 frames, which after removing the 1 sec wait time at the end of each resetting resulted in 213618 frames.
• Constant time (3.79s) returns -20 movies of 15 resetting events and 1 movie of 5 resetting events were combined to a total of 305 resetting events, all preformed on the same particle. The combined movie consists of 149783 frames, which after removing the 1s wait time at the end of each resetting resulted in 143683 frames.
In order to estimate the steady state distribution along the x-axis (or the radial steady state distribution), we first allow a relaxation time T and only then obtain the position distribution by time averaging over the remaining trajectory. To this end, we measured the transient distribution of the particle at a time T , after the onset of the experiment (Fig. S2). We find that for each combined movie a steady state distribution was reached after approximately T = 5/r, where 1/r is the average resetting time. We therefore obtained the steady state trajectory by removing from the beginning of the combined movie a number of frames equal to T = 5/r, i.e., which on average amounts to five resetting events.
FIG. S2. Position distributions in 1D for different transient times since the beginning of the experiment which was conducted at a rate r = 0.2s −1 . It can be seen that after 25s, the probability distribution coincides with the theoretical prediction of the steady state given by Eq. (1) in the main text.

SII. STEADY STATE OF DIFFUSION WITH STOCHASTIC RESETTING
In this section, we derive the theoretical formulas for the steady state position distributions of a Brownian particle diffusing in two dimensions in the presence of stochastic resetting. Here, we assume that the particle starts from the origin and that its motion is reset at a rate r as described in the main text. In what follows, we consider three different return protocols: (i) instantaneous return, (ii) return at a constant velocity v, and (ii) return at a constant time τ 0 .

S5
A. Instantaneous returns: derivation of the theoretical result in Fig. 2b Here, we will consider the motion of a Brownian particle with resetting and instantaneous returns. The steady state distribution for the position of such a particle which is diffusing in d-dimensions in known and given by [S2] where ν = 1 − d 2 , α 0 = r/D and K ν (z) is the modified Bessel function of the second order. Here, we will be interested only in the radial position distribution ρ inst (R = | x|). To derive it, we need to know the Jacobian which is given by Using the Jacobian we find In particular, by substituting d = 2 into Eq. (S3) we obtain the steady state radial density for diffusion with instantaneous resetting in two dimensions This result is corroborated experimentally in Fig. 2b in the main text.

B. Constant velocity returns: derivation of Eq. (1) in the main text
We now consider 2d Brownian motion with resetting and constant radial return velocity Here, v x and v y are the velocity components in the x and y directions; and note that while v is kept constant v x and v y vary depending on the angle v makes with the x-axis. As before, we will be interested in the steady state radial density ρ(R). This will be obtained by advancing on a methodology that was recently developed in [S3, S4] as described below. Diffusion with resetting and non-instantaneous returns has two phases of motion: diffusion and return. We will denote the conditional density of the particle in the diffusive and return phases as ρ diff (R) and ρ ret (R) respectively. The steady state density can then be written as where p c.v. D is the steady state probability to find the particle in the diffusive phase [S4]. This probability is given by average time spent in the diffusive phase average time spent in the diffusive and return phases where τ (R) is the mean time spent in the return phase [S4].
To proceed, we recall that invariance properties of stochastic motion with resetting assert that the conditional density describing the diffusive phase at the steady state, ρ diff (R), must be identical to ρ inst (R) from Eq. (S4) as generally proven in [S4]. We thus have In addition, we have

S6
with τ (R ) = R /v standing for the return time of the particle from a point whose distance from the origin is R . Plugging this result into Eq. (S7) and utilizing Eq. (S8), we find (S10) Finally, it was shown in [S4] that which is our case gives Plugging in the results for ρ diff (R) and ρ ret (R) into Eq. (S6), we arrive at the following expression for the steady state radial density which is Eq. (1) in the main text. Eq. (S13) is in very good agreement with experimental data as shown in Fig. S3.
In Fig. S4, we compare the result in Eq. (S13) to the result in Eq. (S4). We now consider 2d Brownian motion with resetting and constant return time τ (R) = τ 0 . In this case, the steady state probability to find the particle in the diffusive phase is given by and the complementary probability to be in the return phase is given by p c.t. R = 1 − p c.t. D = rτ 0 p c.t. D . As before, due to invariance properties [S4], we can immediately write down the conditional steady state density in the diffusive phase Now, since the particle returns to the origin at a fixed time (regardless of its distance from the origin at beginning of the return phase), the radial return velocity is not constant but rather given by v(R) = R/τ 0 . The conditional probability density to find the particle in the return phase is then given by [S4] ρ ret (R) = rp c.t. (S13)] (blue). Both densities are plotted vs. the radial distance R and the radial velocity v. It can be seen that the blue and red manifolds closely overlap when v is large. This is expected since the large v limit essentially corresponds to a situation where return is almost instantaneous. On the other extreme, for small return velocities, clear difference are seen between the red and blue manifolds. Here: r = 0.05s −1 and D = 0.18µm 2 s −1 .
Substituting Eq. (S4) in above we find where L L L n is modified Struve function of integer order n [S5]. Thus, the total density is simply given by which is Eq.
(2) in the main text. Eq. (S18) is in very good agreement with experimental data as shown in Fig. S5.
We have plotted the steady state radial density Eq. (S18) in a generic phase space in Fig. S6.

A. Derivation of Eq. (4) in the main text
To derive Eq. (4) in the main text we start from Eq. (3) which is a general equation that provides the probability density function of the return time (S18)] (blue). Both densities are plotted vs. the radial distance R and the return time τ0. It can be seen that the blue and red manifolds closely overlap when τ0 is small. This is expected since the small τ0 limit essentially corresponds to a situation where return is almost instantaneous. On the other extreme, for large return times, clear difference are seen between the red and blue manifolds.
where R is the d-dimensional position vector of the particle measured from the origin, τ ( R) is the return time and G 0 ( R, t) is the propagator of the underlying stochastic dynamics. The integration of R is performed over D, the domain of our interest. This formula can be readily applied to our setup and this gives the probability density function of the return time in two-dimensions where R is simply the Euclidean distance of the particle from the origin, τ (R) the return time, G 0 (R, t) = 1 4πDt e −R 2 /4Dt the diffusion propagator in polar coordinates, and f (t) = re −rt the exponential density of the resetting time. Now, since the return time only depends only on the radial distance, we integrate over the angular coordinate θ. In addition, we observe that for a constant radial return velocity v, we have τ (R) = R/v. Integrating over R and t then gives with α 0 = r/D. To get Eq. (4) in the text we simply observe that E(t) = Pt. A simple calculation then yields the probability density for the energy spent per resetting event with E 0 = α −1 0 v −1 P. Eq. (S22) allows us to compute all the moments of E. In particular, the mean energy spent per resetting event is found to be E = πE 0 /2.

SIV. FIRST PASSAGE OF DIFFUSION WITH STOCHASTIC RESETTING
To study the first passage time distribution we performed 5 sets of experiments at different resetting rates and a constant return time of τ 0 = 3.79s the movies were captured at 20 FPS: • r = 0.05s −1 -20 movies of 15 resetting events and 1 movie of 5 resetting events were combined to a total of 305 resetting events, all preformed on the same particle.
• r = 0.0667s −1 -17 movies of 10 resetting events, 2 movies of 15 resetting events and 6 movies of 20 resetting events were combined to a total of 320 resetting events, all preformed on the same particle.
• r = 0.125s −1 -45 movies of 15 resetting events and 7 movies of 5 resetting events were combined to a total of 710 resetting events, all preformed on the same particle.
• r = 0.5s −1 -51 movies of 20 resetting events were combined to a total of 1020 resetting events, all preformed on the same particle.
• r = 1s −1 -44 movies of 30 resetting events, 1 movie of 20 resetting events and 1 movie of 10 resetting events were combined to a total of 1350 resetting events, all preformed on the same particle.
To extract the series of FPT events, we first project the motion of the particle onto the x-axis (Fig. S7a). Then, we identify segments of the trajectory which start at the origin and end at the virtual absorbing wall that is marked in a solid line in Fig. S7a. In Fig. S7b, these segments are marked in cyan and each such trajectory segment stands for a first passage event. First passage times are obtained by taking the time lengths of the first passage events.
FIG. S7. a) Projection of a particle's trajectory onto the x-axis. The position of the virtual absorbing wall at L = 1µm is marked with a solid line. Here, the resetting rate was set to r = 0.05s −1 . b) The same trajectory from a), but with first passage segments colored in cyan. These segments of the trajectory start at the origin and end when the particle first reaches the absorbing wall.

A. Derivation of Eq. (6) in the main text
In this subsection, we provide a derivation of Eq. (6) for the mean first passage time of diffusion with stochastic resetting and a constant return time. To derive this equation, we take advantage of a recently developed framework for home range search [S6] and observe that the mean first passage time in our case can be written as Substituting these expressions in Eq. (S23), we find that the mean first time to a target positioned at L, starting from the origin, is given by which is Eq. (6) in the main text.

SV. ENERGY COST PER FIRST PASSAGE EVENT
In this section, we derive analytical expressions for the energy cost associated with a single first passage event. To this end, we once again require the probability density of the return time. In a first passage time scenario, we can write the probability density function of the return time, following Eq. (7) in the main text, as

S10
where p is the probability that a resetting event occurs before a first passage event, R is a d-dimensional position vector, τ ( R) is the return time from R, and G abs ( R, t) is the reset-free propagator in the presence of the absorbing target. The integration is performed over the domain D. We also need the probability distribution for the number of return events (N r ) per first passage event. This is simply given by a geometric distribution P r(N r = k) = p k (1 − p) , for k = 0, 1, 2, 3, · · · . (S26) Therefore the mean number of return events per first passage is given by N r = p 1−p . We are now left with the task of computing p which is essentially the probability that reset occurs before first passage. This probability reads where we recall that f (t) is the resetting time density.
We are now ready to specialize the general formulas given above to the cases discussed in the main text. Initially, we considered the constant return time protocol. In this case, we have τ ( R) = τ 0 for all R, which is thus also the mean return time. The mean energy spent per first passage event E FP is the mean energy spent per return event (Pτ 0 ) multiplied with the mean number of return events ( N r [S7]. We then considered a constant velocity return protocol with v < v max . In this case, we have τ ( R) = |x|/v, and the mean return time can be computed using Eq. (S25) to give where again we have assumed that resetting times are taken from an exponential distribution with mean 1/r. Thus the mean energy spent per first passage event is given by which was used in the main text. Note that when r → 0 (equivalently α 0 → 0), we find E FP = PL/v. Thus while E FP ≡ 0 for r = 0, when r > 0 one has E FP > PL/v which is essentially the minimum energy required to drag a particle at constant velocity v from the origin to a target located at L.
A. Physical realizability of the constant return time protocol and details of Fig. 5d Observe that there exists a critical value of r * below which the average return velocity in the constant return time protocol exceeds any prefixed value v. To see this, define the mean return velocity for the constant time protocol. This is given by the average return distance divided by τ 0 . The former is given by S11 which in our case reads Thus, the mean return velocity in the constant return time protocol is given byv(r) = |x| /τ 0 , which is a monotonically decreasing function of the resetting rate r. In particular, note thatv(r → 0) = ∞ whilev(r → ∞) = 0.
Comparing this with the velocity v (in the constant return velocity protocol), we conclude that there exists a critical resetting rate for whichv In particular, note that for r > r * , we havev < v. Otherwise for r < r * , we havev > v. The point beyond which the blue line in Fig. 5d was dashed was found by numerically solving Eq. (S32). Also, as the above analysis holds for any v, setting v = v max yields a critical resetting rate r * below which the constant return time protocol is not physically realizable.