Anomalous diffusion and run-and-tumble motion of a chemotactic particle in low dimensions

TL;DR

This paper investigates the anomalous diffusion behaviors of a self-chemotactic particle in low dimensions, revealing superdiffusive, logarithmic, and normal diffusion regimes depending on the interaction type and spatial dimension.

Contribution

It introduces a non-linear, non-Markovian model for chemotactic particles and maps the 1D case to a run-and-tumble dynamics, highlighting new diffusion regimes.

Findings

Superdiffusive MSD with exponent 4/3 in 1D repulsive chemotaxis

Logarithmic growth of MSD in attractive chemotaxis

Diffusion reverts to normal in 3D with a renormalized coefficient

Abstract

We study the stochastic dynamics of a symmetric self-chemotactic particle and determine the long-time behavior of its mean squared displacement (MSD). The attractive or repulsive interaction of the particle with the chemical field that it generates induces a non-linear, non-Markovian effective dynamics, which results into anomalous diffusion for spatial dimensions $d \leq 2$. In one spatial dimension, we map the case of repulsive chemotaxis onto a run-and-tumble-like dynamics, leading to an MSD which, as a function of the elapsed time $t$, grows superdiffusively with exponent $4/3$. In the presence of attractive chemotaxis, instead, the particle exhibits a slowdown, with the MSD growing logarithmically with time. In $d=2$, we find logarithmic aging of the diffusion coefficient, while in $d=3$ the motion reverts standard diffusive behavior with a renormalized diffusion coefficient.

Tables

Table 1. Table 1: Long-time behavior of the MSD ⟨ 𝒓 2 ​ ( t ) ⟩ \langle\bm{r}^{2}(t)\rangle depending on the dimensionality d d and on the sign of the interaction α \alpha with the field. D ​ ( ⋯ ) D(\cdots) is an unknown function.

chemotaxis	MSD
	$d=1$	$d=2$	$d=3$
repulsive ($\alpha >0$)	$t^{4/3}$	$D\left(\ln t\right)t$	$t$
attractive ()	$\ln t$	$t/{\left(\ln t\right)}^2$

Equations

$$\{\begin{array}{l}\gamma \dot{\mathbf{r}}(t)=-\alpha \nabla {\varphi }_{t_R}(\mathbf{r}(t),t)+\mathbf{\zeta }(t),\\ {\partial }_t\varphi (\mathbf{x},t)=\mathcal{D}{\nabla }^2\varphi (\mathbf{x},t)+{\delta }^d(\mathbf{x}-\mathbf{r}(t)),\end{array}$$

$${\varphi }_{t_R}(\mathbf{x},t)={\int }_0^t\mathrm{d}{t}^{\prime }{G}_d(\mathbf{x}-\mathbf{r}(t^{\prime }),t+t_R-t^{\prime }),$$

$$g(t_R)=\{\begin{array}{ll}(\ln t_R)/(8\pi {\mathcal{D}}^2)& \text{for}d=2,\\ -1/(12{\pi }^{3/2}{\mathcal{D}}^2\sqrt{\mathcal{D}{t}_R})& \text{for}d=3.\end{array}$$

$$\{\begin{array}{l}\gamma b^{1-z}\dot{\mathbf{r}}=-b^{d_{\alpha }+d_{\varphi }-1}[\alpha \nabla {\varphi }_{t_R}-{\beta }_{\gamma }\dot{\mathbf{r}}\ln b]+b^{-z/2}\mathbf{\zeta }(t)\\ b^{d_{\varphi }-z}{\partial }_t\varphi (\mathbf{x},t)=b^{d_{\varphi }-2}\mathcal{D}{\nabla }^2\varphi (\mathbf{x},t)+b^{-d}{\delta }^d(\mathbf{x}-\mathbf{r}(t)).\end{array}$$

$$\left[-2+\alpha d_{\alpha }{\partial }_{\alpha }+zt{\partial }_t+{\beta }_{\gamma }(t_R){\partial }_{\gamma }\right]⟨{\mathbf{r}}^2⟩=0.$$

$${\gamma }_E(t)=\{\begin{array}{ll}\gamma -\frac{\alpha }{8\pi {\mathcal{D}}^2}\ln \left(t/t_R\right)& \text{ for }d=2,\\ \gamma -\frac{\alpha }{12{\pi }^{3/2}{\mathcal{D}}^2\sqrt{\mathcal{D}{t}_R}}& \text{ for }d=3.\end{array}$$

$$\dot{r}(t)-\frac{\alpha }{2\sqrt{\pi }}{I}_t[\dot{r}]=\zeta (t),$$

$$⟨r^2(t)⟩\simeq \frac{8\pi T}{{\alpha }^2{a}^2(0)}\ln t+c,$$

$$\mathcal{L}[\delta r](s)=\frac{\alpha \mathcal{L}\left[{\ddot{r}}_a/v_a^3\right](s)}{\left[s+\frac{\alpha }{2}{v}_a-\frac{\alpha }{4}\left(\frac{2v_a^2+4s}{\sqrt{v_a^2+4s}}\right)\right]},$$

$${\lambda }^{ϵ}(f_o)\simeq \frac{\alpha (f_c+ϵf_o)B}{4\sqrt{2\pi TA}}\exp \left[-8A\frac{(f_c+ϵf_o{)}^2}{{\alpha }^{2}T}\right].$$

$$f_o(x,t)=2{ϵ}_1{G}_1(x,t)+4\sum _{i=1}{ϵ}_i{G}_1(x-x_i,t-t_i).$$

$$\{\begin{array}{l}{\partial }_{t}P(r,t)=\frac{{\alpha }^2}{8{\lambda }_0}{\partial }_r^{2}P(r,t)+\frac{\alpha {\lambda }_1}{4{\lambda }_0}{\partial }_r[\bar{f}(r,t)P(r,t)],\\ \bar{f}(r,t)=-\alpha {\int }_0^{t}d{t}^{\prime }\int d{r}^{\prime }P(r^{\prime },t^{\prime }){\partial }_r{G}_1(r-r^{\prime },t-t^{\prime }).\end{array}$$

$$\delta (t_1-t_2)={\int }_{-\infty }^{+\infty }\frac{du}{2\pi }{t}_1^{iu-1/2}{t}_2^{-iu-1/2},$$

$$⟨\dot{r}(t_1)\dot{r}(t_2)⟩=2T{\int }_{-\infty }^{+\infty }\frac{du}{2\pi }\{1-\tilde{\alpha }{I}_{t_1}{\}}^{-1}[t_1^{iu-\frac{1}{2}}]\{1-\tilde{\alpha }{I}_{t_2}{\}}^{-1}[t_2^{-iu-\frac{1}{2}}].$$

$$a(p)=\sqrt{\pi }\frac{\Gamma (p)}{\Gamma (p+1/2)}-\frac{1}{p},$$

$$⟨\dot{r}(t_1)\dot{r}(t_2)⟩=\frac{8\pi }{{\alpha }^2}{\int }_{-\infty }^{+\infty }\frac{du}{2\pi }\frac{t_1^{iu-1}{t}_2^{-iu-1}}{|a(iu){|}^2}.$$

$$⟨r^2(t)⟩=\frac{8\pi }{{\alpha }^2}{\int }_{-\infty }^{+\infty }\frac{du}{2\pi }\frac{4{\sin }^2\left((u/2)\ln \left(t/t_c\right)\right)}{u^2|a(iu){|}^2}+c,$$

$$-\alpha \nabla \varphi (r_s(t),t)=\frac{\alpha }{4\sqrt{\pi }}{\int }_{t_{lt}}^t\frac{r_s(t)-r_s(t^{\prime })}{(t-t^{\prime }{)}^{3/2}}\exp \left[-\frac{r_s(t)-r_s(t^{\prime })}{4(t-t^{\prime })}\right]+f_o(r_s(t),t),$$

$$\frac{\alpha }{4\sqrt{\pi }}{\int }_{t_{lt}}^t\frac{r_s(t)-r_s(t^{\prime })}{(t-t^{\prime }{)}^{3/2}}\exp \left[-\frac{r_s(t)-r_s(t^{\prime })}{4(t-t^{\prime })}\right]=$$

$$\frac{\alpha }{4\sqrt{\pi }}{\int }_{t_{lt}}^t\frac{r_a(t)-r_a(t^{\prime })}{(t-t^{\prime }{)}^{3/2}}\exp \left[-\frac{r_a(t)-r_a(t^{\prime })}{4(t-t^{\prime })}\right]+\frac{\alpha }{4\sqrt{\pi }}{\int }_{t_{lt}}^t\frac{\delta r(t)-\delta r(t^{\prime })}{(t-t^{\prime }{)}^{3/2}}\left[1-\frac{r_a(t)-r_a(t^{\prime })}{2(t-t^{\prime })}\right]\exp \left[-\frac{r_a(t)-r_a(t^{\prime })}{4(t-t^{\prime })}\right].$$

$$\dot{\delta r}(t)={\int }_0^{t}d{t}^{\prime }\gamma (t-t^{\prime })\dot{\delta r}(t^{\prime })+\alpha \frac{{\ddot{r}}_a}{v_a^3}(t),$$

$$\dot{\eta }(t)=-\frac{1}{\tau }\eta (t)+\sqrt{2T}\zeta (t).$$

$$\bar{\tau }=\frac{\sqrt{2\pi T\tau }}{(f_c+ϵf_o)}\exp \left[\frac{(f_c+ϵf_o{)}^2\tau }{2T}\right].$$

$$-v{\partial }_x{\varphi }_{+}(x)-{\partial }_x^2{\varphi }_{+}(x)=\delta (x).$$

$${\partial }_x{\varphi }_{+}(x)=\{\begin{array}{ll}0& \text{for}x<0,\\ -e^{-vx}& \text{for}x>0.\end{array}$$

$${\partial }_t{P}_{\pm }=\pm v{\partial }_r{P}_{\pm }+{\lambda }_0(P_{\mp }-P_{\pm })\pm {\lambda }_1(P_{+}+P_{-})\bar{f}.$$

$$\{\begin{array}{l}{\partial }_{t}P=-v{\partial }_x\bar{s},\\ {\partial }_t\bar{s}=-v{\partial }_{r}P-2{\lambda }_0\bar{s}-{\lambda }_1\bar{f}P.\end{array}$$

$$\bar{f}(r,t)=-\alpha t^{1-2z}{\int }_0^{1}du{u}^z{U}^{\prime }\left(\frac{r}{u{t}^z}\right).$$

$$⟨t_{\mathrm{tbl}}⟩=(1-e^{-{\lambda }_0\bar{t}})/{\lambda }_0.$$

Full text

††preprint: APS/123-QED

Introduction.

Chemotaxis — the ability to move in response to chemical gradients— is a fundamental mechanism of motion in biological and synthetic systems, governing processes such as bacterial and cellular migration [1, 2], the motion of catalytic colloids and swimming droplets [3, 4, 5, 6, 7, 8, 9, 10], and even the activity of individual enzymes [11, 12, 13]. In these systems, the environment is continuously reshaped by the particles themselves, as they emit the chemical species that drive their own motion and create an evolving chemical landscape. This self-chemotactic behavior represents a paradigmatic example of active matter [14], and it has been extensively studied for its ability to induce self-propulsion when the particle maintains a chemical gradient along its surface [15, 3]. Janus colloids, i.e., particles engineered with asymmetric surface activity, represent a standard route to generate such gradients [15, 16, 17], but this built-in asymmetry is not essential. In fact, spontaneous symmetry breaking can occur when the interaction between the colloid and the released chemical is repulsive: in this case, the chemical accumulation around the particle destabilizes its stationary configuration, resulting in autonomous propulsion. In the steady state and in the absence of thermal fluctuations, the particle moves with a constant speed that can be predicted analytically [18, 19, 20, 21]. However, because of their microscopic size, chemotactic particles are affected significantly by stochastic forces. Accordingly, in order to understand their transport properties, it is essential to characterize their mean squared displacement (MSD), defined as the mean $\langle\bm{r}^{2}(t)\rangle$ of the particle position $\bm{r}(t)$, as a function of time $t$. If the chemical produced by the colloid decays in time, the temporal correlations of the velocity of the particle vanish rapidly and the the particle behaves diffusively at large times with modified diffusion coefficient [22, 23]. However, most chemotactic particles emit chemical species with negligible or no decay on the timescale of their dynamics. In this case, modifications of the chemical landscape by the particle motion induce persistent, non-linear memory effects on its dynamics. Persistent memory typically leads to anomalous diffusion, as demonstrated in various low-dimensional systems, including trail-interacting particles [24], self-interacting random walks [25, 26, 27], and passive tracers in active baths [28]. Despite this, existing studies on low-dimensional chemotactic dynamics mostly focus on determining the effective diffusion coefficients [29, 30, 31, 32, 33], leaving the precise functional form of the MSD at long times largely unexplored. In this work, we aim at filling this gap by analyzing a paradigmatic model of a self-chemotactic particle and systematically studying the influence of the spatial dimensionality $d$ and the sign of the (attractive or repulsive) particle-field interaction $\alpha$ on their long-time dynamics. Using a scaling analysis, we identify a critical spatial dimension $d_{c}=2$ separating distinct dynamical regimes. For $d>d_{c}$, the particle exhibits normal diffusion at long times with renormalized diffusion coefficient. For $d\leq d_{c}$, instead, we find anomalous diffusion: superdiffusion for repulsive interactions and subdiffusion for attractive interactions. To further characterize these behaviors, we map the dynamics onto simpler effective models – a generalized Langevin equation with memory for the attractive case and a biased run-and-tumble model for the repulsive one. These mappings allow us to determine analytically the asymptotic behavior of the MSD, yielding logarithmic growth for the attractive case and superdiffusive growth scaling as $t^{4/3}$ in $d=1$. These predictions are summarized in Tab. 1.

The model.

The dynamics of the system is determined by the evolution of the position $\bm{r}$ of the particle and of the concentration field $\phi(x,t)$ of the produced chemical. The particle, with overdamped dynamics, is subject to a phoretic force proportional to the gradient of the chemical field at its current position, to a viscous friction with coefficient $\gamma$, as well as to a stochastic force due to its interaction with a bath in equilibrium at temperature $T$. The field, instead, diffuses in space with a point source located at the colloid position, due to the chemical emission from the particle. The resulting equations are:

[TABLE]

where $\mathcal{D}$ is the diffusivity of the chemical field, $\alpha$ the phoretic constant, and the noise $\bm{\zeta}(t)$ has zero mean and variance $\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=2T\delta(t-t^{\prime})\delta_{ij}$. A similar set of equations appears in diffusion through random media, studied for both passive [34, 35] and active [32, 33] particles. Here, however, chemical production and the absence of decay leads to anomalous diffusion at all parameter values. Note that, in the presence of a point source, $\nabla\phi$ actually diverges at the particle position for $d>1$. Accordingly, one needs to regularize the divergence: this introduces a timescale $t_{R}$, indicated as a subscript of $\phi$ in Eq. (1), which is assumed to be smaller than any other scale of the system. Different regularization schemes introduced in the literature include smoothing the source term in Eq. (1) with a Gaussian [34, 31] (accounting for the finite size of the colloid) or computing the field produced by the particle only up to time $t^{\prime}<t$ [29] (accounting for the time delay between the production of chemical and the phoretic response). Changing regularization scheme does not alter the dynamics (up to numerical factors), provided that the scheme itself does not break the rotational symmetry of the problem. We regularize Eq. (1) as follows: being $G_{d}(\bm{x},t)=(4\pi\mathcal{D}t)^{-d/2}e^{-x^{2}/(4\mathcal{D}t)}$ the solution of the diffusion equation, the field generated at position $\bm{x}$ by the colloid with trajectory $\{\bm{r}(t^{\prime}),0\leq t^{\prime}\leq t\}$ is $\phi(\bm{x},t)=\int_{0}^{t}{\rm d}t^{\prime}G_{d}(\bm{x}-\bm{r}(t^{\prime}),t-t^{\prime})$. Then we define:

[TABLE]

i.e., the field $\phi$ at position $\bm{x}$ but at a later time $t+t_{R}$, assuming that no chemical is released from $t$ to $t+t_{R}$. Considering $\phi_{t_{R}}(\bm{x},t)$ for $t^{\prime}\simeq t$, it can be shown that $\nabla\phi_{t_{R}}(\bm{x},t)$ diverges at $\bm{x}=\bm{r}(t)$ as $t_{R}\rightarrow 0$ for $d\geq 2$. In particular, this divergence is proportional to the particle velocity: $\nabla\phi_{t_{R}}(\bm{r},t)\rightarrow g(t_{R})\dot{\bm{r}}(t)+\mbox{finite}$, where the finite part is independent of $t_{R}$, while $g$ is given by:

[TABLE]

Scaling analysis.

To study the long-time behavior of the system we perform a scaling analysis of Eq. (1) by applying the transformations $\bm{r}\rightarrow b\bm{r}$ and $\bm{x}\rightarrow b\bm{x}$ to the space variables with a scale parameter $b\simeq 1$. Time, field, and the phoretic constant are correspondingly scaled as $t\rightarrow b^{z}t$, $\phi\rightarrow b^{d_{\phi}}\phi$, $\alpha\rightarrow b^{d_{\alpha}}\alpha$, where $z$, $d_{\phi}$, and $d_{\alpha}$ are determind such that Eq. (1) is invariant under the scale transformation. Note that the scale transformation reduces the cutoff $t_{R}$ to $t_{R}/b^{z}$. We use Eq. (3) to restore the original cutoff, by noticing that for $b\simeq 1$ one has $\nabla\phi_{t_{R}/b^{z}}=\nabla\phi_{t_{R}}-zt_{R}\frac{\partial g}{\partial t_{R}}\dot{\bm{r}}(t)\ln b$. The rescaled equations of motion are:

[TABLE]

with $\beta_{\gamma}(t_{R})=\alpha z\frac{\partial g}{\partial\ln t_{R}}$ for $d=3$ or 2 and $\beta_{\gamma}=0$ for $d=1$. Imposing that Eq. (4) is the same as Eq. (1) requires $z=2$, $d_{\phi}=2-d$, $d_{\alpha}=-d_{\phi}$, but this is not sufficient. We also transform the viscosity as $\gamma\rightarrow\gamma+\beta_{\gamma}(t_{R})\ln b$ to account for the short-scale correction to the gradient due to the change in cutoff. Since Eq. (1) is unchanged under the aforementioned transformations, also $\langle\bm{r}^{2}\rangle$ has to be unchanged. Imposing its invariance for $b\simeq 1$ gives the Callan-Symanzik equation:

[TABLE]

Using the method of characteristics, we solve Eq. (5) to express the MSD as $\langle\bm{r}^{2}(t)\rangle=D_{E}(\alpha_{E}(t),\gamma_{E}(t),t_{R})\,t$, where the effective diffusivity $D_{E}$ depends on the effective coupling $\alpha_{E}(t)=\alpha\left(t/t_{R}\right)^{-d_{\alpha}/z}$ and viscosity:

[TABLE]

For $d<2$, $|\alpha_{E}|$ grows with time and thus the MSD depends on the asymptotic form of $D_{E}$ at large $\alpha$, which is unknown from scaling alone. This case will be discussed separately below. For $d\geq 2$ and $\alpha<0$, Eq. (6) shows how the chemical cloud surrounding the particle generally slows down its dynamics by increasing the effective viscosity $\gamma_{E}$(t) of the medium compared to its bare value $\gamma$. Moreover, as for $d>2$ $\alpha_{E}$ vanishes at long times, $D_{E}$ is determined by the Stokes-Einstein relation $D_{E}(0,\gamma_{E},t_{R})=2dT/\gamma_{E}^{2}$ (for $d=2$ this is an approximation which requires small $\alpha$). In particular, for $d=3$, only the chemical produced at times close to $t$ affect the dynamics of the colloid and thus the effective viscosity approaches a constant value which is function of the bare coupling $\alpha$. This is shown in Fig. 1(a). By contrast, in $d=2$, the dynamics is affected by the chemical produced at early times, causing a logarithmic growth of the effective viscosity $\gamma_{E}$ and thus ageing of the system. The corresponding reduction in the effective diffusivity upon increasing time is shown in Fig. 1(b). Finally, at $\alpha>0$ the Stokes-Einstein relation does not correctly capture the dependence of the diffusivity on $\alpha$, as $\gamma_{E}$ vanishes at $\alpha=\alpha_{c}=12\pi^{3/2}\mathcal{D}^{2}\sqrt{\mathcal{D}t_{R}}\gamma$ (in $d=3$) or $\alpha>0$ (in $d=2$), causing a divergence of $D_{E}$.

This change of sign corresponds to the onset of the self-propelling regime of the chemotactic particle, similarly to other active systems [37, 38, 39, 28]. This implies that for $\alpha>0$ and $d=2$, the precise dependence of $D_{E}$ on $\ln t$ remains undetermined. The numerical investigation of this case is also difficult, due to the long timescales required to resolve a logarithmic growth, and it is left for future works. For $d>2$ and $\alpha>0$, instead, the scaling analysis still predicts a diffusive behavior at long times, although with an undetermined diffusion coefficient.

Dynamics in one dimension.

In $d=1$ the chemical interaction generates non-Markovian effects which affect more strongly the particle dynamics. Note that, now, we can set $t_{R}=0$ without encountering divergences. Without loss of generality, we also set $\mathcal{D}=1$ and $\gamma=1$ The dependence on $\mathcal{D}$ and $\gamma$ can be recovered by applying the transformations: $\alpha\rightarrow\alpha/(\gamma\mathcal{D}^{2})$, $T\rightarrow T/(\gamma^{2}\mathcal{D}^{2})$, $t\rightarrow t/\mathcal{D}$, $\phi\rightarrow\phi/\mathcal{D}$ . Since for $d<2$ the scaling analysis is inconclusive, a separate study of the attractive and repulsive interactions is required. Consider, first, the attractive case $\alpha<0$: as a subdiffusive growth of the MSD is expected, the displacement $[r(t)-r(t^{\prime})]^{2}$ grows slower than the elapsed time $t-t^{\prime}$ if the latter is large and/or $T$ is low. Accordingly, we can approximate: $\nabla G_{1}(x,t)\simeq-\frac{1}{4\sqrt{\pi}}\frac{x}{t^{3/2}}$ in the expression of $\nabla\phi_{t_{R}=0}$ following from Eq. (2), to linearize the memory kernel induced by the self chemotaxis. This leads to the generalized Langevin equation for the colloid:

[TABLE]

where $I_{t}$ is an effective retarded friction due to chemotaxis, given by $I_{t}[f]=\int_{0}^{t}dt^{\prime}\left[(t-t^{\prime})^{-1/2}\!-\!t^{-1/2}\right]f(t^{\prime})$. The solution of Eq. (7) (see Ref. [36] for details) implies, at long times,

[TABLE]

where $a(0)=1.39\ldots$ while $c$ is an undetermined constant. Numerical simulations confirm this prediction [36]. In the repulsive case $\alpha>0$ both nonlinearity and memory play a crucial role in the dynamics. In the noiseless limit $T\to 0$, the particle exhibits self-propulsion: in fact, by inserting the ansatz $r(t)=\pm vt$ into Eq. (1), we find two steady-state solutions with speed $v=\alpha/2$, introducing a natural timescale $\tau_{\alpha}=4/v^{2}=16/\alpha^{2}$. At small but finite temperature $T$, the noise induces tumbling events between these two states, inverting the direction of motion. When a tumble occurs, the particle visits again previously explored regions and interacts with the residual chemical field it left behind. This self-generated chemical memory modifies the dynamics and plays a key role in determining its long-time behavior. To investigate this, one needs to understand both the mechanism of tumbling and how it is influenced by the chemical memory. Note that these tumblings — defined as changes in the sign of the velocity — cannot be straightforwardly identified from $\dot{r}(t)$, which is non-differentiable due to the white noise. To single out the slow dynamics, we project the trajectory onto low-frequency modes by convolving Eq. (1) with the exponential kernel $g(t)=\tau^{-1}e^{-t/\tau}$. In particular, we define the coarse-grained position as $r_{s}(t)=(g*r)(t)$, where $(g*f)(t)=\int_{-\infty}^{t}g(t-t^{\prime})f(t^{\prime})$. We choose the convolution timescale to be $\tau\propto\tau_{\alpha}$. We assume that fast fluctuations average out and we approximate the convoluted chemotactic force as $g*\nabla\phi(r(t),t)\simeq\nabla\phi(r_{s}(t),t)$. Under this approximation, $r_{s}(t)$ evolves according to the same dynamics as $r(t)$ in Eq. (1), but with the white noise $\zeta(t)$ replaced by the smoothened noise $\zeta_{s}(t)=g*\zeta(t)$.This filtering allows one to identify the tumbling events from the slow variable $r_{s}(t)$. To understand how the tublings occurr, we split the chemotactic force $-\alpha\nabla\phi(r_{s}(t),t)$ in Eq. (2) into two contributions: one arising from the chemical produced since the last tumble and the other from those released before that time, which we denote by $f_{o}(r_{s}(t),t)$. Note that if $T$ is sufficiently low, tumblings are rare and the chemical field has time to relax in the meanwhile. This causes $f_{o}$ to be a slowly varying function, consistent with the requirement that $r_{s}$ is also slow. For infinitely slow $f_{o}$ and $\zeta_{s}$, the systems admits an adiabatic solution with velocity $\dot{r}_{a}(t)=\epsilon\,\alpha/2+f_{o}(r_{s}(t),t)+\zeta_{s}(t)$, where $\epsilon=\pm 1$ indicates the direction of the self-propelling motion. Expanding $r_{s}(t)=r_{a}(t)+\delta r(t)$ and Laplace-transforming the equation for $\delta r$ we arrive to [36]:

[TABLE]

where $\mathcal{L}$ denotes the Laplace transform and $v_{a}=|\dot{r}_{a}(t)|$. The tumbling instability happens when at least one pole of $\mathcal{L}[\delta r](s)$ has a positive real part, leading to exponential growth of $\delta r(t)$. According to Eq. (9) this requires that the particle is slowed down to $v_{a}<(\sqrt{5}-2)\alpha$ and therefore that $-\epsilon[f_{o}(r_{s},t)+\zeta_{s}(t)]$ exceeds the critical force $f_{c}=(\sqrt{5}-5/2)\alpha$. The rate $\lambda^{\epsilon}(f_{o})$ with which fluctuations in $\zeta_{s}$ cause $f_{o}+\zeta_{s}$ to satisfy this tumbling condition are given by [36]:

[TABLE]

Here, $A$ and $B$ are two fitting parameters: $A$ optimizes the scale separation between fast and slow dynamics (thus fixing the coarse-grained timescale $\tau=A\tau_{\alpha}$), while $B$ accounts for sub-leading corrections of order $O(T)$ in $f_{c}$, $v_{a}$, and $\tau_{\alpha}$. Their values are determined by measuring the rate $\lambda_{0}=\lambda^{\pm}(0)$ of the first tumbling from simulations of Eq. (1) at various $\alpha$ and $T$ and fit Eq. (10). To conclude the derivation of the effective run-and-tumble dynamics we express $f_{o}(x,t)$ in terms of the positions at which previous tumblings occurred. As at low $T$ the separation between two consecutive tumblings is large compared to the typical length scale $l=\mathcal{D}/v$ of the dynamics, we can study $f_{o}$ in the limit $l\rightarrow 0$. We start by considering the gradient $\partial_{x}\phi_{\pm}$ of the chemical field $\phi_{\pm}$ produced by a colloid moving with constant velocity $\pm v$ up to time $t$. As $\phi_{\pm}$ changes significantly only close to the particle position $r(t)$, at large scales one has $\partial_{x}\phi_{\pm}\simeq\mp l^{-1}\delta(x-r(t))$, see Ref. [36]. The chemical field generated by a run-and-tumble trajectory, with tumblings occurring at the space-time points $(x_{i},t_{i})$, can then be decomposed as illustrated in Fig. 2(a). For each tumbling, one adds (blue) and subtracts (red) the field $\phi_{-}$ or $\phi_{+}$ due to a fictitius particle reaching $x_{i}$ from $|x|=\infty$ by moving with a velocity opposite to that of the trajectory right before the tumbling. By pairing these contributions, one finds that the $i$-th tumbling generates a field $\epsilon_{i}(\phi_{+}-\phi_{-})$, where $\epsilon_{i}$ indicates the sign of the velocity along the trajectory immediately before the tumbling. This field acts as a point source of chemical gradient, causing a chemotactic force field $f_{i}(x,t)=4\epsilon_{i}G_{1}(x-x_{i},t-t_{i})$.

$f_{o}$ is then given by the sum of the contributions due to the various tumblings, plus the one of the unpaired trajectory from the initial condition, i.e.,

[TABLE]

Equations (10) and (11) define an effective run-and-tumble dynamics influenced by the chemotactic memory, which we simulate and compare with the original dynamics, as shown in Fig. 3. In particular, by fitting the measured $\lambda_{0}$ with Eq. (10), we obtain $A=0.42$ and $B=8.05$. The quality of the resulting fit is shown in Fig. 3(a) at varying $T$ and $\alpha$. The MSD as a function of $t$ and the probability distribution (PDF) of $r/\sigma(t)$ ($\sigma^{2}(t)=\langle r^{2}(t)\rangle$) at two (large) $t$, resulting from the effective run-and-tumble dynamics and from the chemotactic Eq. (1) are compared in Figs. 3(b) and 3(c), respectively, showing good agreement. In particular, the MSD for both models show the same scaling exponent and PDF at long times, which we now derive for the latter. Let $P(r,t)$ denote the probability of finding the chemotactic particle at position $r$ and time $t$ and $\bar{f}(r,t)=\langle f_{o}(r,t)\rangle$ the force acting on the particle, averaged over all possible past trajectories. Within the mean-field approximation and for large $t$ they satisfy [36]:

[TABLE]

with \lambda_{1}=\frac{\partial\lambda^{+}(f)}{\partial f}\big{|}_{f=0}. In the first equation, the diffusion term stems from the run-and-tumble motion and competes with the advective term due to $f_{o}$. Using the ansatz $P(r,t)=t^{-z}U(r/t^{z})$ in Eq. (12) with $z>1/2$, one finds $z=2/3$, such that $\langle r^{2}(t)\rangle\propto t^{4/3}$ [36]. This prediction is confirmed by the numerical simulations of both the chemotactic model in Eq. (1) and the effective run-and-tumble dynamics given by Eqs. (10) and (11), see Fig. 3(b). Remarkably, this exponent coincides with that one of the “true” self-avoiding random walk (TSAW) [26, 40, 25], a random walk in which the probability of visiting a site decreases with the number of previous visits to it. In fact, since the motion of the particle is superdiffusive, at long times the spreading of the chemical by diffusion is negligible compared to the typical particle displacement: the chemical is effectively frozen at its emission points and this fact turns the run-and-tumble into a TSAW. This is confirmed by the agreement between the PDF obtained from both models and the analytical solution of the TSAW [41, 42], see Fig. 3(c).

Conclusions.

We have shown that non-Markovian, chemically mediated interactions cause anomalous diffusion of self-chemotactic particles for $d\leq 2$. This observation has profound implications for the collective behavior of chemically active matter [43, 44]. While most theoretical descriptions of chemotactic particles rely on the assumption of fast chemical diffusion $\mathcal{D}\rightarrow\infty$ [45, 46], this assumption breaks down for $\alpha>0$ and $d\leq 2$, as particles become superdiffusive. In this case, the relevant scaling limit is $\mathcal{D}\rightarrow 0$, suggesting that collective chemotaxis in low dimensions may belong to a fundamentally different universality class than its higher-dimensional counterpart. Understanding these emergent behaviors remains an open and compelling challenge. An interesting generalization of our work involves replacing the production of chemical in Eq. (1) with an interaction which conserves its overall mass, encompassing the case of a probe non-reciprocally coupled to a bath [47]. Finally, another direction is to consider statistical properties beyond transport, for example examining how memory influences first-passage properties, which are central to search and sensing strategies in active systems.

S.I Derivation of the MSD and simulations for the attractive case in $d=1$

In the attractive case, the self-chemotactic dynamics (1) can be approximated by the generalized langevin equation (7), allowing for analytical progress. In order to derive the expression of the mean-squared displacement (MSD) in Eq. (8) we use the following representation of the delta-function:

[TABLE]

with $t_{1}$, $t_{2}>2$. This expression is obtained by doing a change of variable, which allows one tho write $\delta(t_{1}-t_{2})=t_{1}^{-1/2}t_{2}^{-1/2}\delta(\ln t_{1}-\ln t_{2})$. Then, representing $\delta(\ln t_{1}-\ln t_{2})$ via its Fourier transform, i.e., $\delta(x)=1/(2\pi)\int\!du\,e^{iux}$, one arrives at Eq. (S1).

We now derive the two-time velocity correlation function $\langle\dot{r}(t_{1})\dot{r}(t_{2})\rangle$ by substituting the formal solution $\dot{r}(t)=\{1-\tilde{\alpha}I_{t}\}^{-1}[\zeta]$ of Eq. (7) in the correlator, where we introduced $\tilde{\alpha}=\alpha/(2\sqrt{\pi})$ and the inverse is understood in the sense of operators acting on functions of time. Using Eq. (S1), the linearity of $I_{t}$, and averaging over the white noise we find:

[TABLE]

Although the action of $I_{t}$ [see its definition after Eq. (7)] on a generic function cannot be given an explicit expression, this can be done on monomials. Specifically, one finds $I_{t}[t^{p}]=a(p+1)t^{p+1/2}$, where

[TABLE]

and $\Gamma$ denotes the Gamma function. Conversely, the action of $I_{t}^{-1}$ on monomials is given by: $I_{t}^{-1}[t^{p}]=\/a(p+1/2)t^{p-1/2}$. Moreover, for large $t\gg\alpha^{-2}$ we have $\tilde{\alpha}I_{t}[t^{p}]\gg t^{p}$, and therefore we can approximate $[1-\tilde{\alpha}I_{t}]^{-1}\simeq\/\tilde{\alpha}I_{t}^{-1}$. Using the latter in Eq. (S2), we find:

[TABLE]

To evaluate the aysmptotic behaviour of the MSD we integrate both $t_{1}$ and $t_{2}$ in Eq. (S4) from $t_{c}$ to $t$, where $t_{c}$ is a cutoff time of order $\alpha^{-2}$, before which the large-$t$ approximation is no longer valid. This leads to:

[TABLE]

where $c$ is an undetermined constant containing the short-time contributions. Finally we note that, for large $t$, the integrand in Eq. (S5) is a fast-oscillating function, and thus it accumulates around its maximum at $u\simeq 0$. Approximating $a(iu)$ by its value at $u=0$ and integrating yield Eq. (8).

To confirm our findings, we simulate both Eqs. (1) and (7) for $\alpha<0$ at low and moderate temperatures, yielding the results shown in Fig. S1.

As expected, the MSD for the two dynamics coincide at low $T$. Upon increasing $T$, the short-time behaviors differ, but both models converge to the same asymptotics in the limit $t\rightarrow\infty$ of large times, which in all cases is well described by Eq. (8). The probability distribution (PDF) of $r$ at large $t$ is Gaussian due to the linearity of Eq. (7).

S.II Equation for $\delta r$ in reciprocal space

In order to derive Eq. (9), we decompose the chemotactic force as

[TABLE]

where $t_{lt}$ denotes the time of the last tumbling before $t$ and $f_{o}$ is due to the chemical field which has been produced before $t_{lt}$. We now substitute $r_{s}=r_{a}+\delta r$ in the first term on the r.h.s. of Eq. (S6) and expand in $\delta r$:

[TABLE]

Since $\dot{r}_{a}$ is a slow function of $t$ we can evaluate Eq. (S7) by expanding $r_{a}(t^{\prime})$ in Taylor series in $t^{\prime}$ around $t$. In the first term on the r.h.s. we expand up to the second order, i.e., $r_{a}(t^{\prime})\simeq r_{a}(t)-\dot{r}_{a}(t)(t-t^{\prime})+\frac{1}{2}\ddot{r}_{a}(t)(t-t^{\prime})^{2}$, ignoring contributions of order higher than 2 in the time derivative. In the second term, instead, the expansion is considered only up to the first order, i.e., $r_{a}(t^{\prime})\simeq r_{a}(t)-\dot{r}_{a}(t)(t-t^{\prime})$ since $\delta r$ is already vanishingly small in the adiabatic limit. Moreover, at sufficiently low $T$, $t_{lt}$ in the integrals can be replaced by $=-\infty$ because tumblings become really rare. After some algebraic manipulations, we find the equation for $\delta r$ to be:

[TABLE]

with $\gamma(t)=\frac{\alpha v_{a}}{4\sqrt{\pi}}\left[\frac{1}{2}\Gamma(-1/2,v_{a}^{2}t/4)-\Gamma(1/2,v_{a}^{2}t/4)\right]$, where $v_{a}=|\dot{r}_{a}(t)|$, while $\Gamma(a,b)$ is the incomplete Gamma function. By taking the Laplace transform of both sides of Eq. (S8) and by using the fact that $v_{a}$ can be treated as a constant while calculating the Laplace transform of $\gamma$, one obtains Eq. (9) of the main text.

S.III Derivation of the tumbling rate from the tumbling condition

As discussed in the main text, the tumbling rate $\lambda^{\epsilon}(f_{o})$ is the rate at which the stochastic force $\zeta_{s}(t)$ overcomes (in absolute value) $f_{c}+\epsilon f_{o}$ while having a sign opposite to that of the self-propulsion speed of the particle. In order to compute $\lambda^{\epsilon}(f_{o})$, we note that the slow noise $\zeta_{s}(t)$ — defined in the main text as $\zeta_{s}(t)=\tau^{-1}\int_{-\infty}^{t}dt^{\prime}\,e^{-(t-t^{\prime})/\tau}\zeta(t^{\prime})$, where $\tau$ is a coarse-graining time — can be alternatively written as $\zeta_{s}(t)=\tau^{-1}\eta(t)$, where $\eta(t)$ is an Ornstein–Uhlenbeck (OU) process [48] with

[TABLE]

The tumbling rate $\lambda^{\epsilon}(f_{o})$ is thus the inverse of the mean-first passage time $\bar{\tau}$ of $\eta$ to surpass a barrier at $-\epsilon\tau(f_{c}+\epsilon f_{o})$. At low $T$ this is given by [49]

[TABLE]

As explained in the main text, we choose $\tau$ to be proportional to the intrinsic timescale of the dynamics $\tau_{\alpha}$: $\tau=A\tau_{\alpha}\equiv 16A/\alpha^{2}$, where $A$ is a fitting parameter. By using this definition and by introducing an additional fitting parameter $B$ for the overall amplitude, we find Eq. (10).

S.IV Chemical field emitted by a particle moving with a constant velocity

We derive here the expression of the chemical field $\phi_{+}$ emitted by a particle moving at a constant velocity $v>0$. For convenience, we work in a reference frame comoving with the particle and centered such that the particle is at $x=0$. The chemical field satisfies the equation:

[TABLE]

For $x<0$ boundedness of $\phi_{+}$ requires $\partial_{x}\phi_{+}=0$, while for $x>0$ one has $\partial_{x}\phi_{+}=Ce^{-vx}$. To fix the value of $C$ we integrate Eq. (S11) from $-\delta$ to $\delta$ for $\delta\rightarrow 0$. This give the condition: $\partial_{x}\phi_{+}(-\delta)-\partial_{x}\phi_{+}(\delta)=1$, which fixes $C=-1$ and gives:

[TABLE]

The case of $v<0$ can be obtained analogously, finding $\partial_{x}\phi_{-}(x)=-\partial_{x}\phi_{+}(-x)$. For $v$ which is small compared with the typical distance $l$ traveled by the particle, Eq. (S12) is strongly peaked around $x=0$ and, accordingly, $\partial_{x}\phi_{+}(x)$ can be approximated by a delta function, as reported in the main text.

S.V Mean-field equations and scaling for the run and tumble model

In this section we derive Eq. (12) for the run-and-tumble model. According to Eq. (11), the maximum magnitude of the contribution of a tumbling event to $f_{o}$ scales as $\Delta t^{-1/2}$, where $\Delta t$ is the time elapsed since the tumbling event. This means that the chemotactic force from the past trajectory is (at most) of order $f_{o}\simeq\sqrt{\lambda_{0}}$ and, at low $T$, remains small compared to the critical force $f_{c}\propto\tau_{\alpha}^{-1/2}$, allowing us to linearize the tumbling rates $\lambda^{\pm}(f_{o})$ as $\lambda^{\pm}(f_{o})=\lambda_{0}\pm\lambda_{1}f_{o}$. We then perform a mean-field approximation by replacing the field $f_{o}$ experienced by each particle with its ensemble average $f_{o}(r,t)\simeq\bar{f}(r,t)$. At time $t$, $\bar{f}$ is given by the second of Eq. (12). Note that, while $f_{o}$ is the force due to the chemical field produced up to the last tumbling, we use the fact that, at sufficiently long time $t$, the time interval between two consecutive tumblings is negligible compared to $t$. This allows us to extend the time integral in the second of Eq. (12) to $t$. Consider now the probability $P_{\pm}(r,t)$ of having a particle at position $r$ and time $t$ traveling with velocity $\pm v$. $P_{\pm}(r,t)$ evolves according to the master equation

[TABLE]

We now rewrite this equation in terms of $P=P_{+}+P_{-}$ and $\bar{s}=P_{+}-P_{-}$, obtaining:

[TABLE]

As we are interested in scaling solutions, t large $t$ $\bar{s}$ varies slowly as a function of $t$ and $\partial_{t}\bar{s}\ll 2\lambda_{0}s$. Accordingly, the l.h.s. of the second equation in Eq. (S14) vanishes and $\bar{s}$ is enslaved to $P$ by the resulting equation $\bar{s}=(-v\partial_{r}P-\lambda_{1}\bar{f}P)/(2\lambda_{0})$.

This equation, combined with the first of Eq. (S14) renders the first of Eq. (12). In order to determine the scaling exponent of $\bar{f}$ we then notice that this involves a convolution of $P$ with the Green’s function $G$, which is a function of $r/\sqrt{t}$: As we are searching for solutions of Eq. (12) with $z>1/2$, $P$ is significantly wider than $G$ for large $t$. This means that for $t\rightarrow\infty$, in the spatial integral which determines $\bar{f}$ according to Eq. (12), $G$ can be approximated by a delta function. Using the scaling ansatz for $P(r,t)$ reported in the main text and involving the scaling function $U$, one can integrate in space the expression for $\bar{f}$ and, after changing the integration variable as $t^{\prime}=ut$, one finds

[TABLE]

The remaining terms in Eq. (12) scale with time as $\partial_{t}P\sim t^{-z-1}$, $\partial_{r}^{2}P\sim t^{-3z}$. Since $z>1/2$, $\partial_{r}^{2}P$ is negligible compared to $\partial_{t}P$ for large $t$. Balancing the latter with the advective term yields $z=2/3$.

S.VI Numerical methods

The simulations on the dynamics induced by Eq. (1) presented in the main text are obtained by numerical integration with an Euler scheme. At each timestep, the gradient of the chemical field at position $r(t)$ is computed from the integral representation given by the gradient of Eq. (2). To discretize this integral, we split it into contributions from all previously recorded particle positions $r_{i}$ at discrete timesteps $i$.

Specifically, we approximate the integral as follows: we linearly interpolate the particle trajectory between two consecutive positions, $r_{i}$ and $r_{i+1}$, yielding a piecewise-linear representation and compute the contribution of each linear region separately. For numerical stability and accuracy, the integral contributions are computed differently for older and more recent time steps.

For older contributions (those occurring from timestep 0 up to timestep $i-m$), we employ a simple trapezoidal integration scheme. For the last $m$ contributions, i.e., those closer in time to the current timestep, we utilize an adaptive Gaussian quadrature integration method. In the simulations reported in the main text, we set $m=100$.

In the numerical simulations reported in Fig. 1(a) we set $dt=0.1$, and run $n=500$ instances of the discretized dynamics for $N=10^{5}$ steps for all but the last two points at $\alpha=0.913$ and $\alpha=1.513$. For these two, we reduce $dt$ by a factor $1/32$ to improve convergence close to the critical point. In Fig. 1(b) instead, we start from $dt=10^{-10}$ and run $n=1000$ instances of the dynamics. Then, we increase $dt$ by a factor of $1000$ and coarsen each trajectory by taking the position every $1000$ steps. Each of these trajectories is then used as an initial condition for the coarser dynamics at larger $dt$. We repeat this scaling procedure $6$ time to obtain the data for Fig. 1(b). For Fig. 3(a), we run our simulation at $dt=0.01$, $t_{R}=0$ for $N=10^{5}$ timesteps and $n=1000$ instances for each value of $\alpha$ and $T$. We acquire the first tumbling time $t_{\rm{tbl}}$, defined as the first instance when the velocity, smoothed via a running average over a window $\tau$, changes sign. Trajectories that do not tumble during the simulation are assigned a tumbling time $t_{\rm{tbl}}$ equal to the simulation time $\bar{t}=N*dt$, i.e., $t_{\rm{tbl}}=\bar{t}$. Assuming that tumbling events follow a Poisson distribution with rate $\lambda_{0}$, the mean of $t_{\rm{tbl}}$ is given by

[TABLE]

We compute $\langle t_{\rm{tbl}}\rangle$ from simulation and use Eq. (S16) to estimate $\lambda_{0}$ in Fig. 3. The simulation of the dicretized chemoticatic dynamics for panels (b) and (c) of Fig. 3 are performed with $dt=0.025$, $N=24000$, and $n=1000$. The effective run-and-tumble model is obtained as follows: at each timestep, the force $f_{o}(r(t),t)$ is calculated using Eq. (11). A slow noise term $\zeta_{s}(t)$ is drawn at each step and the particle velocities are inverted whenever the noise exceeds the threshold $-\epsilon(f_{o}(r(t),t)+f_{c})$, where $\epsilon=\pm 1$ is the sign of the current velocity.

Finally, in Fig. S1, both the original dynamics [see Eq. (1)] and its linearized approximation [ in Eq. (7)] are simulated using a timestep $dt=0.5$, with trajectories of $N=1000$ steps each and averaging over $n=15000$ instances. The linearized dynamics is obtained by approximating $\nabla G_{1}(r(t)-r(t^{\prime}),t-t^{\prime})\simeq-\frac{1}{4\sqrt{\pi}}\frac{r(t)-r(t^{\prime})}{(t-t^{\prime})^{3/2}}$ for the gradient in Eq. (2).