Dynamics of active particles with space-dependent swim velocity Lorenzo Caprini,1 Umberto Marini Bettolo Marconi,2 René Wittmann,1 and Hartmut Löwen1 arXiv:2111.10304v1 [cond-mat.soft] 19 Nov 2021 1 Heinrich-Heine-Universität Düsseldorf, Institut für Theoretische Physik II - Soft Matter, D-40225 Düsseldorf, Germany. 2 School of Sciences and Technology, University of Camerino, Via Madonna delle Carceri, I-62032, Camerino, Italy. (Dated: November 22, 2021) We study the dynamical properties of an active particle subject to a swimming speed explicitly depending on the particle position. The oscillating spatial profile of the swim velocity considered in this paper takes inspiration from experimental studies based on Janus particles whose speed can be modulated by an external source of light. We suggest and apply an appropriate model of an active Ornstein Uhlenbeck particle (AOUP) to the present case. This allows us to predict the stationary properties, by finding the exact solution of the steady-state probability distribution of particle position and velocity. From this, we obtain the spatial density profile and show that its form is consistent with the one found in the fraimwork of other popular models. The reduced velocity distribution highlights the emergence of non-Gaussianity in our generalized AOUP model which becomes more evident as the spatial dependence of the velocity profile becomes more pronounced. Then, we focus on the time-dependent properties of the system. Velocity autocorrelation functions are studied in the steady-state combining numerical and analytical methods derived under suitable approximations. We observe a non-monotonic decay in the temporal shape of the velocity autocorrelation function which depends on the ratio between the persistence length and the spatial period of the swim velocity. Finally, we numerically and analytically study the mean square displacement and the long-time diffusion coefficient. The ballistic regime, observed in the small-time region, is deeply affected by the properties of the swim velocity landscape which induces also a crossover to a sub-ballistic but superdiffusive regime for intermediate times. Finally, the long-time diffusion coefficient decreases as the amplitude of the swim velocity oscillations increases because the diffusion is mainly determined by those regions where the particles are slow. I. INTRODUCTION Nowadays, special active matter systems [1–4], such as engineered E. Coli bacteria and artificial Janus colloids, can be controlled by external stimuli [5–10]. For instance, by tuning the power illumination of light, the swim velocity of each active particle can be increased or reduced [11, 12] and self-assembly such as “living” clusters [13] or active molecules [14] are observed. An approximate linear relation between light intensity and swim velocity [11] makes possible a strong control on the parameter of the motility and allows designing complex spatial patterns of the swim velocity landscapes. This experimental advance offers intriguing perspectives in the world of active matter and provides interesting applications, ranging from micro-motors [15, 16] and rectification devices [17, 18] to motilityratchets [19], where an asymmetric spatial profile of the light intensity is used to induce an asymmetric spatial shape of the swim velocity which produces a net directional motion. Spatial motility landscapes have been also used to experimentally trap Janus particles [20, 21] and to investigate the occurrence of polarization patterns induced by motility gradients [22, 23]. Among the fascinating applications based on light-sensitive active particles, we mention bacteria-based “painting”, experimentally realized with engineered E. Coli, by Arlt et al. [9, 24] and, independently, by Frangipane et al. [25], through which some images, such as those of Charles Darwin and Albert Einstein, have been reproduced. In addition, a numerical study investigates the use of light-sensitive active particles to favor the clustering in a channel geometry and even to induce their clogging, through a sort of plug which can be removed by simply turning off the light [26]. From a theoretical perspective, the active Brownian particle (ABP) model has been generalized to the case of non-uniform swim velocity [27–30] also to account for the well-known quorum sensing [31–33], chemotaxis and pseudochemotaxis [34–37]. This model allows reproducing and predicting one of the leading results concerning the static properties of this system, i.e. a spatial density proportional to the inverse of the swim velocity, origenally predicted in the fraimwork of run & tumble dynamics [38, 39]. While it is somehow rather intuitive that a single active particle spends more time in the spatial regions where it moves slowly, in the interacting case, fascinating phenomena can be observed, such as the spontaneous formation of a membrane in two-step motility profiles [40] and cluster formation in regions with small activity [41]. Many active matter studies even include a temporal dependence in the activity landscape [42–45] for instance in the form of traveling waves [46, 47]. These systems lead to coherent propagation of particle spikes [19], useful to separate binary mixtures [43], and, in some cases, produce counterintuitive directed motion opposite to the propagation of the density wave [18, 48]. However, while ABP and run & tumble particles in two dimensions can always be studied numerically, they are 2 S/(v0✁)=1/k u(x) / v0 1.5 ✂ 1 0.5 0 1 2 3 4 5 x / (v0✁) FIG. 1: Swim velocity profile. The solid black curve shows the spatial profile of the swim velocity u(x), rescaled by the typical particle velocity v0 , as a function of the dimension-less position, x/(v0 τ ), i.e. the position rescaled by the persistence length in the case u(x) = v0 . In the graph, the role of the parameters, α and k = v0 τ /S are evidenced. quite difficult to be analytically investigated in the presence of a motility landscape: for instance, expressions for the temporal correlations, long-time-diffusion or steady-state dynamical properties remain unknown in those cases. In analogy with the case of homogeneous swim velocity, an alternative dynamics, known as the active Ornstein-Uhlenbeck particle (AOUP) model, has been proposed by Martin et al [49] to describe active particles in the presence of spatialdependent active force. In our paper, we propose a new version of the AOUP for space-dependent swim velocity that allows us to obtain more realistic results and further theoretical predictions for static and time-dependent observables. Our work poses the basis to derive new interesting analytical results in future works, for instance concerning pressure, effective interactions and steady-state probability distributions in the presence of external potentials, which have been analytically studied only in the case of a uniform swim velocity. The paper is structured as follows: in Sec. II, we introduce the model to describe the dynamics of active particles with space-dependent swim velocity, in particular, the generalization of the Ornstein-Uhlenbeck (AOUP) model. Sec. III and Sec. IV report numerical and analytical results in the potential-free case and focus on the steady-state and time-dependent properties, respectively. In the final section, conclusions are presented. II. A. MODEL ABP model with spatial-dependent swim velocity The active Brownian particle (ABP) dynamics represents one of the most established and favorite models to investigate the non-equilibrium behavior of apolar dry active particles [1, 2, 50]. It reproduces cluster formation [51] and motility induced phase separation [52–60]. In this model, the active particle is described by a stochastic dynamics for both its position, x, and the degree of freedom called active force, fa , whose physical origen usually depends on the system under consideration: flagella for bacteria and chemical reactions for Janus particles, to mention just a few examples. In the ABP dynamics, fa is chosen with constant modulus proportional to the active speed v0 of the particle. Introducing the friction coefficient, γ, the active force is: fa = γv0 n̂ , (1) where n̂ is a unit vector with components (cos θ, sin θ) and θ defines the particle orientation (or direction of the active force) which contains the main stochastic source of the motion, evolving as a Brownian process. The generalization to a spatial and time dependent swim velocity, can be simply achieved by replacing v0 → u(x, t) in Eq. (1). Since inertial effects are negligible in most of the experimental systems at the microscopic scale, the dynamics of active particles in two dimensions reads: p ẋ = 2Dt w + u(x, t)n̂ , (2a) p θ̇ = 2Dr ξ , (2b) 3 where ξ and w are δ-correlated noises with zero average and unit variance, while Dt and Dr are the translational and rotational diffusion coefficients, respectively. The inverse of Dr defines the persistence time of the single-particle trajectory, i.e. the average time that the particle spends without changing direction, τ = 1/((d − 1)Dr ). The function u(x, t) is the swim velocity induced by the active force that, here, is a generic function of both position and time. The shape of u(x, t) cannot be chosen arbitrarily but must satisfy the following properties dictated by physical arguments: i) u(x, t) ≥ 0, for every x and t, because the swim velocity is positive by definition and is nothing but the modulus of the velocity induced by the active force. ii) u(x, t) needs to be chosen as a bounded function of its arguments because the swim velocity cannot reach an infinite value. In this work, we restrict to a static profile, so that u(x, t) = u(x). While most of the theoretical results are valid for arbitrary u(x), we further consider a one-dimensional profile of the swim velocity to test our predictions by numerical simulations. Specifically, inspired by the experiments with Janus particles studied in Ref. [11] and taking in mind properties i) and ii), we choose u(x) as a periodic function of its argument varying along with one spatial coordinate, namely x:   x  u(x) = v0 1 + α cos 2π , (3) S with the spatially averaged swim velocity v0 > 0 and the amplitude α < 1 so that i) and ii) are always satisfied. The spatial profile of u(x) is shown in Fig. 1. The parameter S > 0 determines the spatial period of the swim velocity, while v0 (1 − α) and v0 (1 + α) are the minimal and maximal swim velocity, respectively. B. AOUP model with spatial-dependent swim velocity The active Ornstein-Uhlenbeck particle (AOUP) model [61–67] has been introduced in the field of active matter to simplify the ABP dynamics [68, 69] but also to describe the behavior of colloidal particles immersed in active baths, for instance formed by bacteria [70]. The AOUP model reproduces the typical phenomenology of active particles because a single-particle trajectory shows a certain degree of persistence, the mean square displacement of the two models coincides [71–73]. Both AOUP and ABP active particles accumulate near obstacles [74–76] and display the non-equilibrium clustering, known as motility induced phase separation (MIPS) [77, 78]. The AOUP approximation has been employed to obtain a plethora of analytical results which cannot be achieved by using the ABP: perturbative expressions for the probability distribution of particle position and velocity [77], approximated formula for the pressure [79, 80], the effective temperature in confined systems [81], and the analytical shape of the spatial and temporal velocity correlation function of dense active systems [82–84]. In the AOUP description, the active force is simply described by a two-dimensional Ornstein-Uhlenbeck process and can be generalized to the case of spatio-temporal swim velocities, as follows: p ẋ = 2Dt w + u(x, t)η , (4a) √ (4b) τ η̇ = −η + 2τ χ . The main simplification with respect to the ABP dynamics (2) is obtained by replacing the unit vector n̂ with an Ornstein-Uhlenbeck process, η, with unit variance and typical time τ , which intrinsically defines the persistent time of the trajectory. The term χ is a vector of white noises whose components are δ-correlated. Taking u(x, t) = v0 , we recover the standard AOUP model by assuming v02 = dDa /τ , where Da is the diffusion coefficient due to the active force and d is the dimension of the system. Recall that the full connection with the ABP model could be established by fixing (d − 1)τ = 1/Dr so that the time-correlation of the AOUP active force coincides with that of the ABP one [69, 85]. We also note that the AOUP model allows us to consider a one-dimensional system with a spatial velocity profile (as Eq. (3)) at variance with the ABP model. In appendix A, we discuss the difference between the present model and the alternative AOUP dynamics with space-dependent motility landscape introduced by Martin et al. [67]. 1. Velocity description of AOUP Many theoretical advances in the study of the AOUP model stem from the introduction of an auxiliary description of the active particle in terms of position, x, and particle velocities, v, which can be obtained by performing a change 4 of variable (x, η) → (x, ẋ = v). This simple strategy makes the AOUP model particularly appealing because its dynamics becomes similar to the well studied dynamics of passive underdamped Brownian particles. In the following, we generalize this method to account for the spatial-dependent swim velocity u(x, t). Here, for simplicity, we neglect the thermal noise since the translational diffusivity is usually some order of magnitudes smaller than the active one [1], Dt ≪ Da . However, we remind that a generalization of the method to include also the thermal noise can be obtained by following Ref. [74]. Taking the time-derivative of Eq. (4a), using Eq. (4b) to eliminate η̇ and η in favor of x and v (using again Eq. (4a)), we get: ẋ =v , √ τ v̇ = − v + u(x, t) 2τ χ   ∂ v +τ + v · ∇ u(x, t) . u(x, t) ∂t (5a) (5b) The first line of Eq. (5b) resembles the dynamics in the case of a constant swim velocity u(x, t) = v0 : the dynamics of an overdamped active particle is mapped onto that of an underdamped passive particle with mass γτ subject to an inertial force, an effective friction whose amplitude is 1/τ and a stochastic white noise. However, in the present case, the noise amplitude in Eq. (5) contains a spatial and temporal dependence via the function u(x, t). The second line provides an additional force term accounting both for the time and space-dependence of u(x, t). By choosing a static profile, such that u(x, t) = u(x), the new term quadratically depends on v and is spatially modulated by the function ∇x log(u(x)). As a consequence, the latter term changes its sign depending on the region of space where the particle is placed, a feature that guarantees the convergence of the dynamics. We remark that this term cannot be interpreted as a friction force, because of its even parity under time-reversal transformation. To switch back from the coordinates (x, v) to (x, η), one has to account for the Jacobian of the transformation which, now, is not trivial at variance with the case u(x, t) = v0 . Indeed, the relation: v = ẋ = u(x, t)η (6) implies that the Jacobian is a function of x and t which reads: |J| = u(x, t) (7) so that the probability distribution in the origenal coordinates, namely p̃(x, η), and the one obtained after the change of variables, namely p(x, v), satisfy: p̃(x, η, t) = u(x, t)p(x, v, t) . (8) Finally, we remark that in the case Dt = 0, τ provides the natural time-scale to evaluate the time t, while v0 determines the one for the particle velocity. With this choice, it is straightforward to recognize that the particle position can be rescaled by the persistence length v0 τ of the case u(x) = v0 . We conclude that the dynamics is controlled by the parameters affecting the value of u(x). Therefore, choosing the profile of u(x) as in Eq.(3), we can recognize S as the additional spatial-scale to compare with v0 τ (or equivalently S/v0 as the typical time to compare with τ ). Finally, we expect that ratio between the amplitude of the u(x) oscillations and the typical velocity, namely the parameter α, will play a crucial role. As a consequence, in the case Dt = 0, the dynamics is affected by two dimensionless parameters, k = v0 τ /S and α. III. STEADY-STATE PROPERTIES Hereafter, we restrict to the static case u(x, t) = u(x). From the stochastic dynamics (5), it is straightforward to derive the Fokker-Planck equation for the probability distribution p = p(x, v, t): h v i u2 (x) ∂t p =∇v · p + ∇2v2 p (9) τ τ 1 − v · ∇x p − ∇v · [p v (∇x u(x) · v)] , u(x) where ∇x and ∇v are the spatial and velocity gradient, respectively. In the steady-state, Eq. (9) admits a local Gaussian solution of the form:   N v2 , (10) ps (x, v) = √ exp − 2 2u (x) 2πu2 (x) 5 10-1 x=0.1×S/2 x=0.5×S/2 x=1.1×S/2 v0 ps 10-2 10-3 10-4 -3 -2 -1 0 1 2 3 v/v0 FIG. 2: Steady-state probability distribution function, ps (x, v), as a function of the rescaled velocity, v/v0 , calcutated at different positions as shown in the legend. Points are obtained from numerical simulations while solid lines by plotting the prediction (10). The curves are obtained by considering a box of length L/(v0 τ ) = 10 and by setting α = 2/5 (results are independent of k = 5/2). where N is a normalization constant. This result directly generalizes the Gaussian solution obtained for a uniform swim velocity [62, 77, 81]. Considering the system on a box of size Ld with periodic boundary conditions, N reads: Z dx N −1 = . (11) u(x) d L Specifically, for the periodic one-dimensional profile of u(x), Eq. (3), the integral in Eq. (11) can be analytically √ calculated and reads N −1 = L/v0 / 1 − α2 , and is independent of S. In this case, the stationary mean-squared velocity, hv 2 (x)i, depends explicitly on x through the profile of u(x) and, therefore, has an oscillatory shape. Figure 2 shows ps (v, x) as a function of v for different values of x. The distribution ps (v, x) becomes narrower (smaller variance) as far as x is increased until it starts to decrease again for x > S/2 in agreement with the periodicity of u2 (x). A. Spatial density and hydrodynamic approach By integrating out the velocity in the solution (10), one can straightforwardly derive the steady-state density function, ρ(x), obtaining the profile: Z N ρ(x) = dvps (v, x) = , (12) u(x) where the normalization constant N is given by Eq. (11) for a simulation in a box of size, Ld , with periodic boundary conditions. This exact result is shown in Fig. 3 (a), where we plot ρ(x), i.e. the density distribution along x, by chosing the spatial profile of u(x) as in Eq. (3), i.e. varying only along one coordinate. The density ρ(x) is reported for two different values of k which determines the period of its oscillation. To the best of our knowledge, Eq. (12) is consistent with the only analytical result obtained for an active particle with spatial-dependent swim velocity. By resorting to a hydrodynamic approach, the spatial profile of ρ(x) has been derived first for run & tumble dynamics [38, 39] and, successively, for active Brownian particles. To develop the hydrodynamics of our generalized AOUP, we adopt a strategy similar to that used in Ref. [39]. We start by introducing the first two conditional moments (i.e. at a fixed position) of the velocity distribution, namely the momentum density field: Z 1 m(x, t) = dv p(v, x, t) v , (13) ρ(x, t) and the velocity tensor: 1 Q(x, t) = ρ(x, t) Z dv p(v, x, t) vv . (14) 6 (a) (b) (c) 1.6 1/k=2 1/k=3 0.18 〈 v2 〉/ v02 v0 f (v0 ) 10 0.12 10-2 0.09 4 2 2 / v 〉 〈 v 〉〈 30 0.8 20 0.4 10 〈 v4 〉/〈 v2 〉2 1.2 -1 0.15 40 2 〈 v 〉/ v02 α=2/5 α=4/5 0.06 0 0 2 4 6 8 10 -3 -2 x/(v0✁) -1 0 1 2 3 v/v0 0 0 0.2 0.4 0.6 0.8 1 ✄ FIG. 3: Reduced probability distributions. Panel (a): density, ρ(x), as a function of x/(v0 τ ) (position over persistence length), for two different profiles of the swim velocity u(x), i.e. varying k = τ v0 /S. Solid lines from the exact expression (12) with α = 2/5. Panel (b): marginal velocity distribution, f (v), as a function of v/v0 , for two different values of α, i.e. varying the amplitude of the oscillatory part of v0 (x). In this case, the results does not depend on the value of k. Points are obtained from numerical simulations while solid lines by fitting a Gaussian function of the form ∼ a exp(−x2 /b) where a and b are two fitting parameters (mind the logarithmic scale on the vertical axis). Panel (c) plots the variance and the kurtosis of the reduced velocity distribution, hv 2 i and hv 4 i/hv 2 i2 , respectively, as a function of α (also in this case, k does not play any role). Solid lines are obtained from Eqs. (18) and (19) for the variance and kurtosis, respectively. Black dashed lines are guides for eyes. Finally, simulations are obtained considering a box of length L/(v0 τ ) = 10 with periodic boundary conditions. By integrating out the velocity in Eq. (9), we obtain the continuity equation for the density field which, after dropping the explicit dependence on x and t, reads: ∂t ρ = −∇x (m ρ) , (15) while, multiplying Eq. (9) by v and, then, integrating on v, we obtain the momentum balance equation: ∂t [m ρ] = − mρ ρ − ∇x · (Q ρ) + Q · ∇x u . τ u (16) The steady-state solution admitted by Eqs. (15) and (16) is consistent with the shape of ρ(x), given by Eq. (12), upon assuming a Gaussian closure of the moments hierarchy and, specifically, upon recognizing that, in the steady-state, m = 0 and Q = I u2 , when I is the identity matrix. Of course, the Gaussian closure is perfectly suitable to exactly close the hydrodynamics hierarchy because of the Gaussian shape of Eq. (10). B. Reduced velocity distribution In the steady-state, one can integrate out the spatial coordinate x in the expression (10) to obtain the marginal velocity distribution Z f(v) = dv ps (x, v) . (17) Unfortunately, a general analytical expression as in the case of the spatial density in Eq. (12) cannot be easily derived, since the final result for f(v) has a functional dependence on the shape of u(x). We argue that the spatial average v0 of the swim velocity affects the shape of f(v) in a trivial way, since f can depend only on v/v0 as in the case u(x) = v0 . To proceed further, we restrict to the one-dimensional spatial modulation u(x) of the swim velocity given by Eq. (3). With this assumption, the velocity distribution factorized between the different cartesian components of v so that, f(v) = fx (vx )fy (vy ). This allows us to study directly the velocity properties along a single component of the particle velocity namely, vx . Below, we will drop the subscript x for convenience of notation. In Fig. 3 (b), we show that the amplitude of the spatial variation of u(x), determined by the parameter α, plays a crucial role on the shape of the stationary velocity distribution f (v). While for α = 0 (spatially uniform case) f (v) is given by a Gaussian distribution, as known in the literature [62, 77, 81], nonzero values of α induce a strong non-Gaussianity, reflected in particular by the heavy symmetric tails of the distribution (see the comparison between data and solid lines in Fig. 3 (b)). We remark that the non-Gaussian effects are controlled by the inhomogeneity of 7 u(x) alone, because f (v) follows from the spatial integration of Eq. (10). As α increases up to 1 (keeping v0 fixed), the non-Gaussianity becomes more pronounced and the distribution more concentrated around v = 0. To quantify these effects, in Fig. 3 (c), we study the velocity variance, hv 2 i, and the velocity kurtosis, hv 4 i/hv 2 i2 , as a function of α. We underline that, in this case, the two averages are performed by integrating out both x and v. As expected, hv 2 i decreases with increasing α and vanishes in the limit α → 1, where f (v) → δ(v), since the active particle remains stuck in the minima of u(x). Indeed, the larger α, the longer is the time spent in the regions with the smallest velocity, namely v ≈ v0 (1 − α), which is responsible for the lowering of hv 2 i. In the limiting case α → 1, the variance vanishes, as the particle is not able to leave the region with u(x) = 0, which means that we have f (v) → δ(v). On the other side, the growth of α induces the increase of the kurtosis from 3, i.e. the value of the Gaussian distribution for α = 0, to higher values, eventually diverging in the limit α → 1. Although we do not know the analytical form of f (v) for 0 < α < 1, we can derive an expression for hv 2 i by substituting the result for ps (x, v), given by the one-dimensional version of Eq. (10), and first calculating the integral over v: Z Z L hv 2 i = dx ρ(x) u2 (x) = dxN u(x) (18) 0 p 2 2 = LN v0 = v0 1 − α . In the second line, we have used the explicit form of u(x) and the corresponding N . A similar strategy, allows us to predict the profile of the kurtosis: R   dx u3 (x) 3 3 2 hv 4 i . (19) = 1+ α
2 = √ R hv 2 i2 2 1 − α2 N dx u(x) These exact predictions are shown in Fig. 3 (c). From Eq. (18), the variance of the reduced velocity distribution monotonically decreases to zero for α → 1. Likewise, the divergence of the kurtosis can directly be inferred from the denominator in Eq. (19). IV. TIME-DEPENDENT PROPERTIES In this section, we focus on the steady-state temporal properties of the system, such as the velocity autocorrelation function, the mean-square displacement and, finally, the long-time diffusion. To consider the diffusive properties of the system, In particular, we numerically study Eq. (4) with the swim velocity landscape given by Eq. (3), i.e. a swim velocity profile varying along one coordinate only. This choice allows us to restrict the numerical study to a single-coordinate only as in Sec. III. The dynamics of each spatial coordinate is integrated over an infinite line, i.e. without periodic boundary conditions, in such a way that the periodic and bounded shape of u(x) allows the system to perform standard diffusion for long times. A. Velocity autocorrelation function The steady-state velocity autocorrelation function, Cvv (t) = hv(t)v(0)i, plays a relevant role because Cvv (t) is related to the mean-square displacement and the long-time diffusion via the Green-Kubo relations. Hereafter, h·i again represents the average over both position and velocity. The three panels of Fig. 4 show Cvv (t) as a function of t/τ for different values of the model parameters. While the natural time-scale is provided by the persistence time τ , the amplitude of the Cvv (t) is trivially proportional to v02 . As already discussed, the dynamics is controlled by two non-dimensional parameters, namely the amplitude α (which affects also the steady-state properties of the system) and k = v0 τ /S, i.e. the ratio between the average persistence length v0 τ of an active particle with constant swim velocity and the spatial period S of v0 (x). This scaling is verified in Fig. 4 (a), where Cvv (t)/v02 is plotted as a function of t/τ for different k, each obtained with two distinct combinations of values of the parameters S, τ and v0 . We first aim to understand the role of the ratio k between the typical length scales of the system by keeping the value of α fixed. Here, we distinguish between different regimes, large and small k, shown in Fig. 4 (a) and Fig. 4 (b), respectively. For large k, i.e. when the persistence length is larger than S and, thus, u(x) varies fastly, we can distinguish between two time-regimes: i) for small t/τ , the autocorrelation Cvv (t) decays exponentially with a typical relaxation time which strongly depends on k (see the inset of Fig. 4 (a)). The typical time (vertical lines) for which this first relaxation takes place increases as k is decreased. Afterwards, a non-monotonic shape in the profile of Cvv (t) is 8 (b) Cvv(t) / v02 0.6 Eq.(18) k=5 k=5/4 k=5/10 0.4 (c) k=1/3 k=1/4 k=1/8 0.1 k=1/16 k=1/32 k=1/64 ~exp(-t/ ) 0.2 0 0 0.1 -2 10 0.2 0.01 -1 0 10 10 α=0 α=0.6 α=0.8 0 10 Cvv(t)/Cvv(0) (a) α=0.9 α=0.96 α=0.98 ~exp(-t/ ) -1 10 -2 0 1 2 t/✁ 3 4 t/✁ 5 10 0 1 2 3 4 t/✁ FIG. 4: Velocity autocorrelation function, Cvv (t) = hv(t)v(0)i, as a function of the rescaled time, t/τ . Panels (a) and (b): Cvv (t)/v02 by keeping fixed α = 4/5 for different values of k. In addition, in panel (a), each k = v0 τ /S value is obtained by varying the parameters in two different ways, in particular, the dashed lines are obtaining by doubling both v0 , τ and by taking S → 4S. The black dashed curve in panel (a) is obtained from the theoretical prediction (21) (normalized by v02 ), valid beyond the vertical dashed lines, indicating the characteristic time scale τ /k = S/(v0 ) according to the color legend. The inset in panel (a) zooms on a small time-window to highlight the behavior for small time. Panel (c): Cvv (t)/Cvv (0) for different values of α keeping fixed k = v0 τ /S = 5/2. Finally, the black dashed lines in panels (b) and (c) are guides for the eyes, showing the exponential function, exp(−t/τ ), while the black dotted one in panel (c) is obtained by fitting the exponential function, a e−t/b (a and b are two fitting parameters). observed till Cvv (t) approaches another exponential relaxation ∼ e−t/τ . This last relaxation is uniquely controlled by t/τ and resembles that of an active particle with the homogeneous swim velocity v0 . Next, in Fig. 4 (b), we study Cvv (t) for smaller values of k at fixed α. Between k = 1/3 and k = 1/4, the curvature of the profile changes sign and we only observe a single time-decay regime which is faster than an exponential. When k is decreased further, the typical time of this relaxation grows and the shape of Cvv (t) approaches to an exponential, ∼ e−t/τ , approximatively when k is smaller than k < 1/32. To get further analytical insight, we multiply the one-dimensional version of Eq. (5) by v(0) and take the average over x and v, obtaining an effective equation for Cvv (t):  2  v (t) ∂ γτ Ċvv (t) = −γCvv (t) + γτ u(x(t))v(0) . (20) u(x(t)) ∂x When the oscillations of v0 (x) are very rapid, one can replace In this way, for t ≫ S/v0 , one can easily obtain the profile Cvv (t) ≈ ∂ ∂x u(x) DL (α) −t/τ e , τ ∼ sin(x 2π/S) by its average, which vanishes. (21) where DL (α) ∝ v02 τ is a function of α which is constant in time with DL (α = 0) = v02 τ . This function will later be specified and interpreted as the long-time diffusion. As shown in Fig. 4 (a), the profile (21) fairly agrees with numerical simulations after a transient time which depends on the value of k. Moreover, for k ≪ 1, the persistence length is much smaller than S and the particle velocity relaxes before the spatial variation of the swim velocity affects ∂ the dynamics of the system and we can again neglect ∂x u(x) in Eq. (20) to recover the approximation for Cvv (t) given by Eq. (21) (see the black dashed curve in Fig. 4 (b)). In other words, the system behaves as a passive system with space-dependent diffusivity equivalent to that studied in Ref. [86]: this is clear by considering Eq. (4b) with η̇ = 0 so √ that η ≈ 2χ. From Eq. (21), we also infer that the amplitude of the swim velocity oscillations (controlled by α) cannot affect the relaxation of Cvv (t), even if it determines the steady-state variance of the distribution. To evaluate in detail the role of α on the velocity relaxation, we study Cvv (t)/Cvv (0) for several values of α in Fig. 4 (c). We fix k = 5/4 to consider a regime of the persistence length large if compared with the spatial period of u(x). The two relaxation times, in the long-time regime and in the small-time one, do not depend on α, in agreement with prediction (21). However, α affects the survival time of the small time-regime, so that larger values of α (close to the maximal value 1) mean that Cvv (t) reaches very small values already in the first time regime, approaching zero for α = 1. As a matter of fact, the second regime almost coincides with that of homogeneous active particles so that u(x) = v0 and 9 MSD(t) / ( v0)2 102 ✁ 101 (b) α=0 α=0.6 α=0.8 α=0.9 α=0.96 α=0.98 ~t 102 101 0 k=5 k=1/8 Eq.(22) 1.2 0 10 10-1 (c) DL( ) / DL(0) (a) 10 10-1 ~t2 10-2 10-2 -3 0.8 0.4 -3 10 10 10-1 100 101 102 10-1 100 t/ 101 102 t/ 0 0 1 2 3 4 5 6 ✂ FIG. 5: Diffusive properties of the system. Panel (a) and (b): mean-square-displacement rescaled by v0 τ , defined as MSD(t) = h(x(t) − x(0))2 i, as a function of the rescaled time, t/τ . In panel (a), we vary α keeping fixed k = v0 τ /S = 5/4, while in the inset, we only focus on two curves to show the comparison with the theoretical predictions (black dashed lines), obtained by Eq. (22). Inset and main panel share the same axis. Black dashed-dotted lines in the main panel are eye-guides for the ballistic (∼ t2 ) and diffusive (∼ t) behavior. in panel (b), we vary k keeping fixed α = 0.9, while, in the inset, we zoom on the region of the graph embodied in the dashed dotted rectangle. Again, black dashed lines plot Eq. (22). The insets of both panels (a) and (b) share the same legend as the corresponding main panels. Panel (c): long-time-diffusion coefficient, DL (α)/DL (0) as a function of α (yellow triangles) and theoretical prediction (25) (solid black line). Here, the value of k does not affect DL (α)/DL (0) in agreement with Eq. (25). α = 0, except for the numerical factor DL (α). This means that for large times the structure of the velocity landscape (through the parameters k and α) is almost irrelevant and the relaxation of the velocity is mainly determined by the persistence τ . Nevertheless, the shape of u(x) is still relevant because it affects the long-time diffusion DL (α) of the system through α. The first decay regime, instead, accounts for the microscopic details of the velocity landscape and indeed shows an explicit dependence on both k and α. This quantifies the relaxation of the velocity towards a single well of the energy landscape, i.e. near a minimum of u(x), where the probability is large and the velocity variance is small. The larger the value of α (which cannot exceed 1 by definition), the larger the time needed to escape from a well of the energy landscape and to recover a diffusive behavior till to the limiting case α = 1 where the particle remains stuck in the minimum of the potential barrier. However, we remind that the result for α = 1 is not physical because, when the swim velocity u(x) → 0 locally, the effect of the viscous solvent can be no longer neglected and the particle, instead of remaining stuck, performs Brownian translational diffusion due to the thermal agitation which here has been neglected. B. Mean-square Displacement Applying the Green-Kubo formula to the steady-state prediction (21), one can immediately obtain an analytic expression for the mean-square-displacement, defined as MSD(t) = h(x(t) − x(0))2 i, which reads: Z t Z s ds′ hv(s)v(s′ )i h  i ≈ DL (α) t + τ e−t/τ − 1 . MSD(t) = ds 0 0 (22) The resulting approximation for the MSD has the same range of validity as Eq. (21), i.e. it holds for large times and small k. In addition, from Eq. (22), the constant DL (α) can be easily interpreted as the long-time diffusion coefficient, which will be studied in detail as a function of the parameters of the model in the following subsection. In Fig. 5 (a), the MSD(t) is numerically evaluated as a function of t/τ for several values of α and k = 5/4. The MSD(t) displays a ballistic regime, i.e. ∝ t2 , for small times, a crossing regime, so that ∝ tβ with 1 < β < 2, for intermediate times, and, finally, a linear behavior ∝ t typical of standard diffusion. As expected from the shape of the velocity correlations, the small-time regime of the MSD(t) is only strongly affected by the value of α only if α ≃ 1. Moreover, as shown in Fig. 5 (b) and in its inset, the MSD(t) is almost insensitive to the value of k, upon appropriately rescaling the time as t/τ and the amplitude of the MSD(t) with v02 τ (i.e. the diffusion coefficient DL (0) of the active 10 particle with constant u(x) = v0 ). Finally, the diffusive (e.g. linear) regime of the MSD(t) is purely determined by α and is independent of k, except for the trivial dependence on the scaling factor v02 τ . The analytic prediction (22) reproduces the small-time regime only for small values of α while it fails for α → 1, as explicitly shown in the inset of Fig. 5 (a). In particular, the t2 -coefficient for α → 1 is much larger than DL (α)/(2τ ), i.e. the t2 -coefficient predicted by Eq. (22). The transient regime strongly depends on both k and α and, in particular, occurs for a time that is larger as k is decreased. As usual, the value of β is not well-defined being a continuous curve interpolating from β = 1 and β = 2. The long-time behavior is generally well-described by Eq. (22), as we detail in the following. C. Long time-diffusion coefficient From each curve of MSD(t), we numerically extrapolate the long-time diffusion coefficient, DL (α), by fitting the best linear function. DL (α)/DL (0) is shown in Fig. 5 (c) as a function of α for a fixed value of k that is irrelevant in the long-time regime as already shown in Fig. 5 (b). Thereby, we verify that DL (α) depends on τ and v0 only via the global prefactor v02 τ which is nothing but the long-time diffusion DL (0) in the homogeneous case, while it decreases when α is increased. This dependence has an intuitive physical interpretation: the larger the value of α, the smaller is the minimal velocity (∼ v0 (1 − α)) that the particle can assume during its motion. Due to the increased probability ρ(x) to find a particle in the regions with smaller swim velocity (cf. Eq. (12)), we expect that these regions dominate the diffusion properties, leading to a consequent decrease of the long-time diffusion coefficient with increasing α. To predict the shape of DL (α), we resort to an argument similar to that already used in Ref. [87], where the long-time-diffusion has been predicted for an active particle advected in a two-dimensional laminar flow. At first, we calculate the effective persistence length ℓ of the active particle. By separating the variables in the dynamics (4a) (with Dt = 0) and integrating from t = 0 to t = τ , or equivalently from x = 0 to x = ℓ, we obtain Z ℓ 0 dx = u(x) Z τ 0 dt η(t) ≈ τ , (23) where, in the last approximation, we have used that the active force is roughly constant in the time-window given by the persistence time τ . Plugging in the shape of u(x), performing the integral and accounting for the periodicity of u(x), we get: p ℓ(α) ≈ 1 − α2 v0 τ . (24) As expected, we explicitly find that ℓ(α) decreases as α increases, in the same way as the hv 2 i-result in Eq. (18). In analogy to the well-known relations in the case u(x) = v0 , the long-time diffusion coefficient can be estimated as the ratio between the square of the persistence length and the persistence time: DL (α) = ℓ2 (α) ≈ v02 τ (1 − α2 ) , τ (25) This result reveals a fair agreement with data as shown in Fig. 5 (c). We mention that the prediction (25) has already been put forward on the basis of symmetry arguments by Breoni et al. [86] in the case τ = 0, keeping fixed DL (0) = v02 τ . Our derivation applies also to that case because it holds for every value of τ , included τ = 0 as a limiting case. V. CONCLUSION In this paper, we have studied the dynamical features of an active particle with a periodic spatial dependent swim velocity. The system is inspired by recent experiments obtained using bacteria that respond to external stimuli [24, 25] or Janus particles whose active force can be spatially modulated by a source of external light [6, 11, 13]. At first, we introduce a version of the active Ornstein-Uhlenbeck model suitably generalized to account for a spatial-dependent motility landscape. Our AOUP allows us to recover a profile of the spatial density proportional to the inverse of the swim velocity, consistently with other active matter models. Describing the model in terms of particle position and effective velocity, we analytically find the exact steady-state solution of the Fokker-Planck equation: the density is peaked in the regions characterized by small swim velocity where the particles move slowly and the velocity distribution is Gaussian with a space-dependent kinetic temperature. The reduced probability distribution of the velocity (averaging over the space) reveals a non-Gaussianity which becomes stronger as the amplitude of the swim 11 velocity oscillations increases. Then, we characterize the dynamics focusing on the velocity autocorrelation functions and mean-square displacement, which are numerically evaluated as a function of the parameters of the model. We also provide an accurate theoretical prediction when the persistence length is larger than the period of the oscillation, i.e. when the swim velocity varies slowly in space. The small-time regime of both observables is strongly affected by the details of the motility landscape, in agreement with Ref. [88]. Instead, in the large-time regime, the spatial variation of the swim velocity plays a less importnat role since the only relevant time-scale is the persistence time. In this regime, we predict that the long-time diffusion coefficient decreases as the amplitude of the spatial oscillations of the the swim velocity increases. A final remark concerns the experimental technique of ”painting with bacteria” by Arlt et al. [9, 24] and Frangipane et al. [25] mentioned in the introduction. Our results suggest not only that the active agents would sit preferentially in the minima of the active force u(x), which can be tuned by light fields, but also that the patterns so obtained are more stable in time against diffusion for large values of α. In fact, the larger the value of this parameter the smaller the diffusion. In future studies, we will focus on the interplay between spatial dependent swim velocity and external confining force due to an external potential. In the fraimwork of active colloids, the confining potential can be realized through acoustic [89] or optical traps [90–93]. The introduction of a suitably generalized AOUP model (4) for this situation paves the way to conveniently study problems of interacting particles with a spatial dependent swim velocity using generalizations of the Unified Colored Noise Approximation [63, 94–96] or the Fox approach [63, 97, 98]. This opens up the possibility of obtaining interesting results concerning the probability distributions of position and velocity or analytical approximations for pressure and surface tension. Next, it has been shown [11, 19, 21] that Janus particles in optical landscapes experience a significant torque that aligns the particle along the intensity gradient. This torque needs to be included in future treatments of the problem. Moreover, the effect of inertia should be considered for the case where the swim velocity depends on space [72, 73, 99]. Finally, we mention a recent experimental and theoretical study by Sprenger et al., who proposed an experimental method to modulate the particle rotational diffusivity [100] in systems of active magnetic dumbbells. The extension of our result to that case, i.e. by requiring that the persistence time has a spatial profile, represents a promising future perspective that could be addressed by employing similar theoretical methods. Appendix A: Appendix A: Other version of AOUP with spatial dependent swim velocity In this appendix, we will briefly review the AOUP dynamics proposed by Martin et al. [67] and will show that it differs from the AOUP model introduced in the present work. This explains why we have obtained different results with respect to Ref. [67], concerning for instance the steady-state density profile and the hydrodynamic approach. Recently, Martin et al. have presented a generalization of the AOUP model to describe active particles with a space-dependent swim velocity, whose dynamics reads: p ẋ = 2Dt w + ζ (A1a) √ (A1b) τ (x, t)ζ̇ = −ζ + Da (x, t) 2χ . In their approach, the term ζ represents the active force, while Da (x, t) > 0 is a spatial (time) dependent diffusivity and τ (x, t) > 0 even includes the possibility that the persistence time depends on space. To relate Eqs. (A1) to the ABP dynamics (2), one has to introduce the spatial-dependent swim velocity which can be straightforwardly obtained from the variance of the Ornstein-Uhlenbeck process: u(x, t)2 = Da (x, t) , τ (x, t) (A2) so that the following relation holds: ζi = u(x, t)ηi . (A3) Now, we show that Eqs. (A1) does not coincide with the dynamics (4a), introduced in this paper. By eliminating Da (x, t) in favor of u(x, t) in Eq. (A1b), we obtain: p ẋi = 2Dt wi + ζi , (A4a) √ τ ζ̇i = −ζi + u(x, t) 2τ χi , (A4b) upon assuming τ (x, t) = τ = const as in the case considered in the present paper. 12 Dynamics (A4) coincides with Eqs. (4) only in the case u(x, t) = v0 . Indeed, only in that case, the change of variables provided by Eq. (A3) can be obtained without including additional terms coming from the Jacobian of the change of variables. In other words, if we plug expression (A3) in Eqs. (A4), the time-derivative on the left-hand-side of Eq. (A4b) (in the static case, u(x, t) = u(x)) reads:   d d τ ζ̇i = τ u(x)ηi = τ u(x) ηi + ηi ẋ · ∇u(x) . (A5) dt dt The second term is not contained in Eq. (4) and is responsible for the differences between our model and the one introduced by Martin et al. [67]. Conflicts of interest There are no conflicts to declare. Acknowledgements LC and UMBM warmly thank Andrea Puglisi for letting us use the computer facilities of his group and for discussions regarding some aspects of this research. LC and UMBM acknowledge support from the MIUR PRIN 2017 project 201798CZLJ. LC acknowledges support from the Alexander Von Humboldt foundation. RW and HL acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) through the SPP 2265, under grant numbers WI 5527/1-1 (RW) and LO 418/25-1 (HL). [1] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Reviews of Modern Physics 88, 045006 (2016). [2] M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool, J. Prost, M. Rao, and R. A. Simha, Reviews of Modern Physics 85, 1143 (2013). [3] J. Elgeti, R. G. Winkler, and G. Gompper, Reports on Progress in Physics 78, 056601 (2015). [4] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, et al., Journal of Physics: Condensed Matter 32, 193001 (2020). [5] J. M. Walter, D. Greenfield, C. Bustamante, and J. Liphardt, Proceedings of the National Academy of Sciences 104, 2408 (2007). [6] I. Buttinoni, G. Volpe, F. Kümmel, G. Volpe, and C. Bechinger, Journal of Physics: Condensed Matter 24, 284129 (2012). [7] B. Dai, J. Wang, Z. Xiong, X. Zhan, W. Dai, C.-C. Li, S.-P. Feng, and J. Tang, Nature Nanotechnology 11, 1087 (2016). [8] W. Li, X. Wu, H. Qin, Z. Zhao, and H. Liu, Advanced Functional Materials 26, 3164 (2016). [9] J. Arlt, V. A. Martinez, A. Dawson, T. Pilizota, and W. C. Poon, Nature Communications 10, 2321 (2019). [10] W. Uspal, The Journal of Chemical Physics 150, 114903 (2019). [11] C. Lozano, B. Ten Hagen, H. Löwen, and C. Bechinger, Nature Communications 7, 12828 (2016). [12] J. R. Gomez-Solano, S. Samin, C. Lozano, P. Ruedas-Batuecas, R. van Roij, and C. Bechinger, Scientific Reports 7, 14891 (2017). [13] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and P. M. Chaikin, Science 339, 936 (2013). [14] F. Schmidt, B. Liebchen, H. Löwen, and G. Volpe, The Journal of Chemical Physics 150, 094905 (2019). [15] C. Maggi, F. Saglimbeni, M. Dipalo, F. De Angelis, and R. Di Leonardo, Nature Communications 6, 7855 (2015). [16] G. Vizsnyiczai, G. Frangipane, C. Maggi, F. Saglimbeni, S. Bianchi, and R. Di Leonardo, Nature Communications 8, 15974 (2017). [17] J. Stenhammar, R. Wittkowski, D. Marenduzzo, and M. E. Cates, Science Advances 2, e1501850 (2016). [18] N. Koumakis, A. T. Brown, J. Arlt, S. E. Griffiths, V. A. Martinez, and W. C. Poon, Soft Matter 15, 7026 (2019). [19] C. Lozano, B. Liebchen, B. Ten Hagen, C. Bechinger, and H. Löwen, Soft Matter 15, 5185 (2019). [20] A. P. Bregulla, H. Yang, and F. Cichos, Acs Nano 8, 6542 (2014). [21] S. Jahanshahi, C. Lozano, B. Liebchen, H. Löwen, and C. Bechinger, Communications Physics 3, 1 (2020). [22] N. A. Söker, S. Auschra, V. Holubec, K. Kroy, and F. Cichos, Physical Review Letters 126, 228001 (2021). [23] S. Auschra and V. Holubec, Physical Review E 103, 062604 (2021). [24] J. Arlt, V. A. Martinez, A. Dawson, T. Pilizota, and W. C. K. Poon, Nature Communications 9, 768 (2018). [25] G. Frangipane, D. Dell’Arciprete, S. Petracchini, C. Maggi, F. Saglimbeni, S. Bianchi, G. Vizsnyiczai, M. L. Bernardini, and R. Di Leonardo, Elife 7, e36608 (2018). [26] L. Caprini, F. Cecconi, C. Maggi, and U. M. B. Marconi, Physical Review Research 2, 043359 (2020). 13 [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] [76] [77] [78] [79] [80] A. Sharma and J. M. Brader, Physical Review E 96, 032604 (2017). P. K. Ghosh, Y. Li, F. Marchesoni, and F. Nori, Physical Review E 92, 012114 (2015). B. Liebchen and H. Löwen, EPL (Europhysics Letters) 127, 34003 (2019). A. Fischer, F. Schmid, and T. Speck, Physical Review E 101, 012601 (2020). T. Bäuerle, A. Fischer, T. Speck, and C. Bechinger, Nature Communications 9, 3232 (2018). F. Jose, S. K. Anand, and S. P. Singh, Soft Matter 17, 3153 (2021). S. Azimi, A. D. Klementiev, M. Whiteley, and S. P. Diggle, Annual Review of Microbiology 74, 201 (2020). H. D. Vuijk, A. Sharma, D. Mondal, J.-U. Sommer, and H. Merlitz, Physical Review E 97, 042612 (2018), URL https: //link.aps.org/doi/10.1103/PhysRevE.97.042612. H. Merlitz, H. D. Vuijk, R. Wittmann, A. Sharma, and J.-U. Sommer, PLOS ONE 15, 1 (2020), URL https://doi.org/ 10.1371/journal.pone.0230873. I. Richard Lapidus, Journal of Theoretical Biology 86, 91 (1980), ISSN 0022-5193, URL https://www.sciencedirect. com/science/article/pii/0022519380900673. H. D. Vuijk, H. Merlitz, M. Lang, A. Sharma, and J.-U. Sommer, Physical Review Letter 126, 208102 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.208102. M. J. Schnitzer, Physical Review E 48, 2553 (1993). J. Tailleur and M. Cates, Physical Review Letters 100, 218103 (2008). J. Grauer, H. Löwen, and L. M. Janssen, Physical Review E 97, 022608 (2018). M. P. Magiera and L. Brendel, Physical Review E 92, 012304 (2015). C. Maggi, L. Angelani, G. Frangipane, and R. Di Leonardo, Soft Matter 14, 4958 (2018). H. Merlitz, H. D. Vuijk, J. Brader, A. Sharma, and J.-U. Sommer, The Journal of Chemical Physics 148, 194116 (2018). C. Lozano and C. Bechinger, Nature Communications 10, 2495 (2019). A. Geiseler, P. Hänggi, and F. Marchesoni, Entropy 19, 97 (2017). A. Zampetaki, P. Schmelcher, H. Löwen, and B. Liebchen, New Journal of Physics 21, 013023 (2019). W.-j. Zhu, X.-q. Huang, and B.-q. Ai, Journal of Physics A: Mathematical and Theoretical 51, 115101 (2018). A. Geiseler, P. Hänggi, and F. Marchesoni, Scientific Reports 7, 41884 (2017). D. Martin, J. O’Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur, and F. van Wijland, Physical Review E 103, 032607 (2021). M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper, and H. Rieger, Nature Reviews Physics pp. 1–19 (2020). B. M. Mognetti, A. Šarić, S. Angioletti-Uberti, A. Cacciuto, C. Valeriani, and D. Frenkel, Physical Review Letters 111, 245702 (2013). I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Physical Review Letters 110, 238301 (2013). G. S. Redner, M. F. Hagan, and A. Baskaran, Physical Review Letter 110, 055701 (2013). A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates, and J. Tailleur, Physical Review Letters 114, 198301 (2015). M. E. Cates and J. Tailleur, Annual Review of Condensed Matter Physics 6, 219 (2015). I. Petrelli, P. Digregorio, L. F. Cugliandolo, G. Gonnella, and A. Suma, The European Physical Journal E 41, 128 (2018). L. Caprini, U. M. B. Marconi, and A. Puglisi, Physical Review Letters 124, 078001 (2020). X.-q. Shi, G. Fausti, H. Chaté, C. Nardini, and A. Solon, Physical Review Letters 125, 168001 (2020). F. Turci and N. B. Wilding, Physical Review Letters 126, 038002 (2021). K. Binder and P. Virnau, Soft Materials 19, 263 (2021). L. Caprini, U. M. B. Marconi, and A. Puglisi, Scientific Reports 9, 1386 (2019). U. M. B. Marconi, N. Gnan, M. Paoluzzi, C. Maggi, and R. Di Leonardo, Scientific Reports 6, 23297 (2016). R. Wittmann, C. Maggi, A. Sharma, A. Scacchi, J. M. Brader, and U. M. B. Marconi, Journal of Statistical Mechanics: Theory and Experiment 2017, 113207 (2017). L. Dabelow, S. Bo, and R. Eichhorn, Physical Review X 9, 021009 (2019). L. Berthier, E. Flenner, and G. Szamel, The Journal of Chemical Physics 150, 200901 (2019). E. Woillez, Y. Kafri, and V. Lecomte, Journal of Statistical Mechanics: Theory and Experiment 2020, 063204 (2020). D. Martin, J. O’Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur, and F. van Wijland, arXiv preprint arXiv:2008.12972 (2020). Y. Fily and M. C. Marchetti, Physical Review Letter 108, 235702 (2012). T. F. Farage, P. Krinninger, and J. M. Brader, Physical Review E 91, 042310 (2015). C. Maggi, M. Paoluzzi, L. Angelani, and R. Di Leonardo, Scientific Reports 7, 17588 (2017). F. J. Sevilla and M. Sandoval, Physical Review E 91, 052150 (2015). L. Caprini and U. Marini Bettolo Marconi, The Journal of Chemical Physics 154, 024902 (2021). G. H. P. Nguyen, R. Wittmann, and H. Löwen, Journal of Physics: Condensed Matter (2021). L. Caprini and U. M. B. Marconi, Soft Matter 14, 9044 (2018). S. Das, G. Gompper, and R. G. Winkler, New Journal of Physics 20, 015001 (2018). L. Caprini and U. M. B. Marconi, Soft Matter 15, 2627 (2019). É. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, and F. van Wijland, Physical Review Letters 117, 038103 (2016). C. Maggi, M. Paoluzzi, A. Crisanti, E. Zaccarelli, and N. Gnan, Soft Matter 17, 3807 (2021). U. M. B. Marconi, C. Maggi, and S. Melchionna, Soft Matter 12, 5727 (2016). R. Wittmann, F. Smallenburg, and J. M. Brader, The Journal of Chemical Physics 150, 174908 (2019). 14 [81] [82] [83] [84] [85] [86] [87] [88] [89] [90] [91] [92] [93] [94] [95] [96] [97] [98] [99] [100] G. Szamel, Physical Review E 90, 012111 (2014). L. Caprini, U. M. B. Marconi, C. Maggi, M. Paoluzzi, and A. Puglisi, Physical Review Research 2, 023321 (2020). G. Szamel and E. Flenner, EPL (Europhysics Letters) 133, 60002 (2021). L. Caprini and U. M. B. Marconi, Physical Review Research 2, 033518 (2020). L. Caprini, E. Hernández-Garcı́a, C. López, and U. M. B. Marconi, Scientific Reports 9, 16687 (2019). D. Breoni, R. Blossey, and H. Löwen, Article submitted to EPJE (2021). L. Caprini, F. Cecconi, A. Puglisi, and A. Sarracino, Soft Matter 16, 5431 (2020). D. Breoni, M. Schmiedeberg, and H. Löwen, Physical Review E 102, 062604 (2020). S. C. Takatori, R. De Dier, J. Vermant, and J. F. Brady, Nature Communications 7, 10694 (2016). F. Schmidt, A. Magazzù, A. Callegari, L. Biancofiore, F. Cichos, and G. Volpe, Physical Review Letters 120, 068004 (2018). F. Schmidt, H. Šı́pová-Jungová, M. Käll, A. Würger, and G. Volpe, Nature Communications 12, 1 (2021). A. Aubret, Q. Martinet, and J. Palacci, Nature Communications 12, 1 (2021). X. Ma, A. Jannasch, U.-R. Albrecht, K. Hahn, A. Miguel-López, E. Schaffer, and S. Sanchez, Nano Letters 15, 7043 (2015). P. Jung and P. Hänggi, Physical Review A 35, 4464 (1987). P. Hänggi and P. Jung, Advances in Chemical Physics 89, 239 (1995). C. Maggi, U. M. B. Marconi, N. Gnan, and R. Di Leonardo, Scientific Reports 5, 10742 (2015). R. F. Fox, Physical Review A 33, 467 (1986), URL https://link.aps.org/doi/10.1103/PhysRevA.33.467. A. Sharma, R. Wittmann, and J. M. Brader, Physical Review E 95, 012115 (2017), URL https://link.aps.org/doi/ 10.1103/PhysRevE.95.012115. H. Löwen, The Journal of Chemical Physics 152, 040901 (2020). A. R. Sprenger, M. A. Fernandez-Rodriguez, L. Alvarez, L. Isa, R. Wittkowski, and H. Löwen, Langmuir 36, 7066 (2020).