# Recoil of a driven tracer in a correlated medium

Marcin Piotr Pruszczyk, Davide Venturelli, Andrea Gambassi

► To cite this version:

Marcin Piotr Pruszczyk, Davide Venturelli, Andrea Gambassi. Recoil of a driven tracer in a correlated medium. Journal of Statistical Mechanics: Theory and Experiment, 2025, 2025 (11), pp.113202. 10.1088/1742-5468/ae0d7b hal-05402091

HAL Id: hal-05402091

https://hal.sorbonne-universite.fr/hal-05402091v1

Submitted on 6 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Recoil of a driven tracer in a correlated medium

To cite this article: Marcin Piotr Pruszczyk et al J. Stat. Mech. (2025) 113202

View the [article online](#) for updates and enhancements.

You may also like

- Self-interacting processes via Doob conditioning

Francesco Coghi and Juan P Garrahan

- Quench dynamics of entanglement entropy under projective charge measurements: the free fermion case

Riccardo Travaglino, Colin Rylands and Pasquale Calabrese

- Phase diagram of two-component mean-field Bose mixtures

Oskar Stachowiak and Pawel Jakubczyk
PAPER: *Classical statistical mechanics, equilibrium and non-equilibrium*

# Recoil of a driven tracer in a correlated medium

Marcin Piotr Pruszczyk<sup>1,2,*</sup> ID, Davide Venturelli<sup>3,4</sup> ID
and Andrea Gambassi<sup>1,2</sup> ID

E-mail: mpruszcz@sissa.it

Received 17 June 2025

Accepted for publication 24 September 2025

Published 19 November 2025

Online at stacks.iop.org/JSTAT/2025/113202

https://doi.org/10.1088/1742-5468/ae0d7b

<figure><img src="image_2.png" /></figure>

CrossMark

**Abstract.** We study the stochastic dynamics of a Brownian particle after it is suddenly released from a harmonic trap moving with constant velocity through a fluctuating correlated medium, described by a scalar Gaussian field with relaxational dynamics and in contact with a thermal bath. We show that, after the release, the particle exhibits recoil, i.e. it moves in the direction opposite to the drag. As expected, this effect vanishes if the field equilibrates instantaneously. The final value of the average position of the particle is reached algebraically in time in the case of conserved dynamics of the field or for non-conserved dynamics at the critical point. Our predictions are expected to be relevant, at least qualitatively, to driven colloidal particles in liquid media close to critical points.

**Keywords:** colloids, bio-colloids and nano-colloids, fluctuation phenomena, stochastic particle dynamics, viscoelasticity

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ SISSA—International School for Advanced Studies, via Bonomea 265, 34136 Trieste, Italy</footnote>

<footnote n="2">² INFN, via Bonomea 265, 34136 Trieste, Italy</footnote>

<footnote n="3">³ Laboratoire de Physique Théorique de la Matière Condensée, CNRS/Sorbonne Université, 4 Place Jussieu, 75005 Paris, France</footnote>

<footnote n="4">⁴ Laboratoire Jean Perrin, CNRS/Sorbonne Université, 4 Place Jussieu, 75005 Paris, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
# Contents

<table><tbody><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. The model</td><td>5</td></tr><tr><td>2.1. Dynamics</td><td>6</td></tr><tr><td>3. Position of the particle</td><td>8</td></tr><tr><td>3.1. Perturbative expansion (weak-coupling approximation)</td><td>9</td></tr><tr><td>3.2. Steady-state average particle position</td><td>10</td></tr><tr><td>3.3. Recoil</td><td>11</td></tr><tr><td>3.4. Long-time behavior</td><td>15</td></tr><tr><td>4. Non-equilibrium steady-state configuration of the field</td><td>19</td></tr><tr><td>4.1. The noiseless limit</td><td>20</td></tr><tr><td>4.2. The steady-state field configuration</td><td>21</td></tr><tr><td>5. Conclusions</td><td>26</td></tr><tr><td>Acknowledgments</td><td>28</td></tr><tr><td>Appendix A. Perturbative calculation of the correction to the particle's position</td><td>28</td></tr><tr><td>A.1. Driven particle</td><td>28</td></tr><tr><td>A.2. Recoil after the particle is released</td><td>31</td></tr><tr><td>Appendix B. Analysis of the recoil</td><td>33</td></tr><tr><td>B.1. Small-<i>v</i> behavior of <i>F</i><sup>A/B</sup></td><td>33</td></tr><tr><td>B.2. Long-<i>t</i> and large-<i>ξ</i> behavior of <i>F</i><sup>A/B</sup></td><td>35</td></tr><tr><td>B.3. Exact expressions for the critical recoil with model <i>A</i> dynamics</td><td>37</td></tr><tr><td>Appendix C. Field configuration with a driven particle</td><td>38</td></tr><tr><td>C.1. Model <i>A</i></td><td>38</td></tr><tr><td>C.2. Model <i>B</i></td><td>39</td></tr><tr><td>References</td><td>40</td></tr></tbody></table>

## 1. Introduction

Ever since the pioneering works of Einstein and Langevin [1, 2], the dynamics of a mesoscopic particle immersed in a simple medium is commonly described in terms of the Langevin equation for the particle position $X(t)$, featuring an instantaneous friction and a random force in the form of a Gaussian white noise. The validity of such a description hinges on the assumption of timescale separation, i.e. that the evolution of
the medium occurs over timescales much shorter than those characterizing the motion of
the particle. In contrast, complex media with macroscopic relaxation times manifestly
violate this assumption [3]. Viscoelastic fluids [4] may serve as an example: indeed,
storage and dissipation of energy within their complex microstructure typically result in
macroscopically long stress-relaxation times. Consequently, dragging a colloidal particle
through such a fluid drives the medium out of equilibrium, which in turn affects the
statistics of the particle position [3], as observed in active microrheology experiments
[5–10].

The overdamped generalized Langevin equation (GLE [11, 12]) offers a description of such a system on a coarse-grained scale by modeling the interaction between the medium and the particle via a friction $\int_{-\infty}^{t} dt' \Gamma(t-t') \dot{\mathbf{X}}(t')$—characterized by the memory kernel $\Gamma(t)$—together with a colored Gaussian noise whose temporal correlations ensure detailed balance. Such a description, though computationally convenient due to the linearity of the evolution equation, has its limitations. First, in general, the functional form of $\Gamma(t)$ is not known *a priori*: in practice, it is often assumed to be of the form $\Gamma(t) = \sum_i a_i e^{-t/\tau_i}$, where the coefficients $\{a_i, \tau_i\}$ are obtained by fitting experimental or numerical data. For experiments involving viscoelastic media, such an analysis was performed for example in [13, 14]. This simplified model—in which each term of the sum may be heuristically visualized as due to a fictitious *bath particle* harmonically coupled to the tracer—is actually able to capture the effective particle dynamics within the linear response regime. Conversely, as the system is driven further out of equilibrium, it turns out that an increasing number of terms has to be included in the sum to accurately fit the data [15, 16], possibly with $a_i < 0$. Second, in constructing the GLE from first principles (see e.g. [12]), the linear projection of the degrees of freedom of the medium that are being coarse grained may lead to a memory kernel that depends not only on the interaction between the particle and the medium, but also on other details of the system—such as the strength of external (confining) forces acting solely on the particle [17, 18].

These considerations suggest the need for more sophisticated models which are able
to describe the particle dynamics in complex media beyond the linear-response regime,
including the peculiar phenomena exhibited by tracer particles in this context. One of
them is the *recoil* displayed by an optically confined particle dragged at constant velocity
through a viscoelastic medium [8, 13, 19–21]: after switching off the confining force, the
particle turns out to move backwards, i.e. in the direction opposite to the previous
drag, before coming to rest. This somehow unexpected fact can actually be understood
heuristically as follows. While the particle is driven through the medium, it locally
perturbs its structure; accordingly, the average position of the particle—which reaches
a steady-state value in the frame of reference co-moving with the trap—is determined
by the mechanical balance between the confining force due to the optical trap, the
friction, and the overall force exerted on the particle by the nonhomogeneous structure
acquired by the medium. This is schematically represented in figure 1 (in which the
field $\phi$ describes the relevant thermally fluctuating field within the medium, as clarified
further below). Generically, the latter force tends to increase the total friction, pulling
the particle in the direction opposite to the drag. When the trap is suddenly switched
off, this force *does not* vanish instantaneously: consequently, the particle is temporarily
pulled backward, until the medium has time to rearrange and relax to equilibrium.
Recoil of a driven tracer in a correlated medium

<figure><img src="image_4.png" /></figure>

**Figure 1.** Schematic representation of a Brownian particle (black sphere), driven by a moving optical trap (providing the potential represented by the dark yellow line) and interacting with a fluctuating correlated field $\phi$ (dark blue line), which is preferably enhanced in the vicinity of the particle (see equation (4) with $\lambda > 0$). (a) While the particle is driven, i.e. for $t < 0$, the field attains an average configuration (blue line) which is stationary and spatially modulated in the reference frame co-moving with the trap. At time $t = 0^{-}$, the particle is subject to the harmonic force $\mathbf{F}_{\kappa}$ exerted by the moving trap, to the field-induced force $\mathbf{f}$ defined in equation (8) which pulls the particle along the gradient of the field, and to the friction force $\mathbf{F}_{\nu}$. (b) At time $t = 0^{+}$, the particle is released from its confinement and thus it starts moving backwards (recoil), due to the fact that $\mathbf{f}$ remains non-vanishing due to the slow relaxation of the field, and that the friction changes both its sign and amplitude $\mathbf{F}_{\nu} \mapsto \tilde{\mathbf{F}}_{\nu} = -\mathbf{f}$.

The recoil described above essentially stems from two key ingredients: (i) a spatial structure acquired by the complex medium, and (ii) a long relaxation time, which renders the effective particle dynamics non-Markovian. A simple model that incorporates both these features can be formulated in terms of the coupled overdamped dynamics of a particle and a scalar Gaussian field $\phi(\mathbf{x}, t)$, characterized by a macroscopic correlation length $\xi$ and a finite relaxation time [18, 22–34, 52]. Modeling the interaction between the field and the particle via a coupling that is linear in $\phi$ allows one to integrate out the latter exactly, yielding the effective (non-linear, non-Markovian) dynamics of the particle. This procedure connects the emergent memory terms that describe the motion of the particle with the underlying properties of the medium. Originally introduced to describe tracer diffusion in a thermally fluctuating environment [25], this minimal model has been more recently interpreted as a simplified description of colloidal particles immersed in a medium close to criticality, such as a binary liquid mixture in the vicinity of the critical point of its demixing transition [35–39]—with a particular interest in the emergence of fluctuation-induced critical Casimir forces [40–43]. In such media, characterized by extended spatial and temporal correlations, $\phi$ should be regarded as the order parameter describing a second-order phase transition [44] (e.g. the relative concentration of the two components in the case of a binary liquid mixture). More broadly, such spatio-temporal correlations also affect the dynamics of inclusions
in lipid membranes [45–47], microemulsions [48–50], or defects in ferromagnetic systems [22–27].

A recent theoretical study [31] investigated the motion of an overdamped particle driven through a correlated medium, within the framework described above. Notably, it was shown that the particle may exhibit oscillatory relaxation modes (generally forbidden in Markovian overdamped systems), analogous to those experimentally observed with a colloidal particle driven through a viscoelastic medium [15]. In order to further investigate the phenomenological similarity between viscoelastic and correlated media, in this work we examine the possibility of observing a correlation-induced analogue of the elastic recoil observed in experiments [8, 19–21]. To do so, in section 2 we consider a particle subject to a harmonic potential moving with a constant velocity $v$, coupled linearly to a scalar Gaussian field which follows relaxational dynamics [51]. After the system attains a non-equilibrium steady state at long times in the reference frame co-moving with the trap, the harmonic trap is abruptly switched off, releasing the particle from its original confinement. The effective, nonlinear dynamics of the particle obtained by integrating out the field degrees of freedom is then solved perturbatively in the particle-field coupling in section 3. The subsequent average trajectory of the particle turns out to exhibit a recoil, which we characterize analytically and verify via numerical simulations of the system. Both the cases of locally conserved and non-conserved dynamics of the field are considered, also in the critical limit in which the corresponding correlation length and relaxation time of the field diverge. We complement our study by the analysis of the steady-state field configuration presented in section 4, and finally in section 5 we summarize our results and present our outlook.

## 2. The model

We consider a system consisting of a tracer, modeling a colloidal particle, and a correlated medium, with the effective Hamiltonian [52]

$$ \mathcal{H}[\phi, \mathbf{X}] = \mathcal{H}_{\phi}[\phi] + \mathcal{H}_{\text{int}}[\phi, \mathbf{X}] + \mathcal{U}(\mathbf{X}, t). \quad (1) $$

Here $\mathbf{X} = (X_1, X_2, \dots, X_d)$ is the position of the particle in $d$ spatial dimensions, which is subject to a time-dependent potential

$$ \mathcal{U}(\mathbf{X}, t) = \frac{\kappa}{2}(\mathbf{X} - \mathbf{v}t)^2 \vartheta(-t), \quad (2) $$

describing a harmonic trap with stiffness $\kappa$. The trap has a center that moves with constant velocity $\mathbf{v}$, and it is switched off abruptly at $t=0$, as encoded by the Heaviside step function $\vartheta(-t)$. The correlated medium is modeled by a fluctuating scalar field $\phi(\mathbf{x}, t) \in \mathbb{R}$ characterized by the quadratic Hamiltonian [44]

$$ \mathcal{H}_{\phi}[\phi] = \int d^d\mathbf{x} \left[ \frac{1}{2}(\nabla\phi)^2 + \frac{1}{2}r\phi^2 \right]. \quad (3) $$

The parameter $r > 0$ determines the spatial extension of correlations of the field, and is related to the correlation length $\xi$ via $\xi = r^{-1/2}$. In this respect, the thermally
fluctuating medium undergoes a second-order phase transition at its critical point $r = 0$, characterized by a diverging correlation length $\xi$ and relaxation time. Additionally, the particle is coupled linearly to the field in equation (1) by

$$ \mathcal{H}_{\text{int}}[\phi, \mathbf{X}] = -\lambda \int d^d\mathbf{x} \, \phi(\mathbf{x}) V(\mathbf{x} - \mathbf{X}). \qquad (4) $$

In passing, we note that such a coupling is invariant under translations in space. The interaction kernel $V(\mathbf{x})$ actually models the shape of the particle and is normalized by requiring that its integral over the whole space is one. Accordingly, the strength of the interaction between the field and the tracer is determined solely by the value of the parameter $\lambda$. Note that the coupling (4) breaks the $\phi \mapsto -\phi$ symmetry of $\mathcal{H}_{\phi}[\phi]$ and thus, for positive $\lambda V(\mathbf{x})$, it favors configurations of the field in which $\phi$ is enhanced in the vicinity of the particle$^5$. Additionally, $V(\mathbf{x})$ is assumed to be spherically symmetric and characterized by a single lengthscale, i.e. $V(\mathbf{x}) = \hat{V}(|\mathbf{x}|/R)$, where $R$ plays the role of the effective radius of the particle. The Gaussian kernel

$$ V(\mathbf{x}) = \frac{1}{(2\pi R^2)^{d/2}} \exp\left(-\frac{\mathbf{x}^2}{2R^2}\right) \qquad (5) $$

may serve as an example. In this effective description of the interaction between the particle and the field introduced in equation (4), the ‘particle’ is identified with the interaction kernel $V(\mathbf{x})$, and the field has non-zero values in its interior. Hence, the model introduced here characterizes the coupling of the particle to the field in terms of two effective parameters, namely the strength $\lambda$ and range $R$ of the interaction.

## 2.1. Dynamics

The dynamics of the particle is described by the overdamped Langevin equation

$$ \dot{\mathbf{X}} = -\nu \nabla_{\mathbf{X}} \mathcal{H} + \xi, \qquad (6) $$

where $\nu$ is the mobility of the particle, and the stochastic noise $\xi$ is a Gaussian variable with vanishing mean and variance given by (for simplicity the Boltzmann constant is set to unity here and henceforth, i.e. $k_B \equiv 1$)

$$ \langle \xi_i(t) \xi_j(t') \rangle = 2\nu T \delta_{ij} \delta(t-t'), \qquad (7) $$

where $T$ is the temperature of the equilibrium heat bath within which the particle moves. Note that, assuming $\lambda V(\mathbf{x}) > 0$, the force $\mathbf{f}[\phi, \mathbf{X}]$ exerted on the particle due to the coupling (4), given by

$$ \mathbf{f}[\phi, \mathbf{X}] = -\nabla_{\mathbf{X}} \mathcal{H}_{\text{int}}[\phi, \mathbf{X}] = \lambda \int d^d\mathbf{x} V(\mathbf{x} - \mathbf{X}) \nabla_{\mathbf{x}} \phi(\mathbf{x}), \qquad (8) $$

<footnote>$^5$ Recall that $\phi$ can be interpreted, for instance, as the relative concentration of the two components of a binary liquid mixture close to its demixing point. In this case, the enhancement of the field around the particle may indicate preferential adsorption of one of the two components by the surface of the particle.</footnote>
pulls the particle along the gradient of the field. For future convenience, we note that equation (8), and consequently the dynamics of the particle, exhibit invariance under {$\phi, \lambda$} $\mapsto \{-\phi, -\lambda\}$. The evolution of the field $\phi$ is governed by the relaxational dynamics [51]

$$ \partial_t \phi(\mathbf{x}, t) = -D_\alpha (i\nabla)^\alpha \frac{\delta\mathcal{H}(\phi, \mathbf{X})}{\delta\phi(\mathbf{x}, t)} + \eta(\mathbf{x}, t), \qquad (9) $$

where $D_\alpha$ is the mobility of the field, and $\alpha \in \{0, 2\}$. Here, for $\alpha = 0$, the dynamics of the field does not obey any local conservation, while for $\alpha = 2$, the field is locally conserved by the dynamics. In fact, equation (9) can then be cast in the form $\partial_t \phi = -\nabla \cdot \mathbf{J}(\mathbf{x}, t)$, where $\mathbf{J}$ is the corresponding current. We note that the cases $\alpha = 0$ and $\alpha = 2$ correspond to the so-called model A and model B, respectively, according to the classification of [44]. In the current analysis, however, we neglect the self-interaction term $\propto \phi^4$ present in the original formulation of these models, and study them within their Gaussian approximation. Note that the physical dimension [$D_\alpha$] of the field mobility $D_\alpha$ depends on the choice of the dynamics. In particular, [$D_\alpha$] = $\mathcal{L}^z \mathcal{T}^{-1}$, where $z = 2 + \alpha$ is the dynamical critical exponent of the field, $\mathcal{L}$ is a unit of length, and $\mathcal{T}$ is a unit of time.

Both the field and the particle are assumed to be in contact with the same heat reservoir characterized by the temperature $T$. Consequently, the noise $\eta(\mathbf{x}, t)$ acting on the field $\phi(\mathbf{x}, t)$ is a Gaussian variable with vanishing mean and characterized by the variance

$$ \langle \eta(\mathbf{x}, t) \eta(\mathbf{x}', t') \rangle = 2D_\alpha T (i\nabla)^\alpha \delta^{(d)}(\mathbf{x} - \mathbf{x}') \delta(t - t'). \qquad (10) $$

The correlations of the noise given by equations (7) and (10) satisfy the fluctuation-dissipation relation. Accordingly, in the absence of external driving, i.e. for $\mathbf{v} = 0$, the joint dynamics of the particle and the field is invariant under time-reversal, and the joint probability distribution function of the position of the particle and the field configuration in the equilibrium steady state is given by the Boltzmann distribution $\propto \exp(-\mathcal{H}[\phi, \mathbf{X}]/T)$, where $\mathcal{H}$ is the Hamiltonian in equation (1).

In order to solve for the dynamics of the field as a function of the particle's position, it will prove useful to express the field in terms of its Fourier transform as

$$ \phi(\mathbf{x}, t) = \int \frac{d^d\mathbf{q}}{(2\pi)^d} e^{i\mathbf{q}\cdot\mathbf{x}} \phi_{\mathbf{q}}(t), \quad \phi_{\mathbf{q}}(t) = \int d^d\mathbf{x} \ e^{-i\mathbf{q}\cdot\mathbf{x}} \phi(\mathbf{x}, t), \qquad (11) $$

where the wavevector $\mathbf{q}$ labels each mode of the field. In these terms, the dynamics of the field imposed by equations (1), (3), (4), (9), and (10) takes the form

$$ [\partial_t + D_\alpha q^\alpha (q^2 + r)] \phi_q(t) = \lambda D_\alpha q^\alpha V_q e^{-i\mathbf{q}\cdot\mathbf{X}(t)} + \eta_q(t), \qquad (12) $$

$$ \langle \eta_q(t) \eta_{q'}(t') \rangle = 2D_\alpha q^\alpha T \delta(t-t') (2\pi)^d \delta^{(d)}(\mathbf{q}+\mathbf{q}'). \qquad (13) $$
Here, $V_q$ is the Fourier transform (see equation (11)) of the interaction kernel introduced in equation (4). Similarly, the dynamics of the particle position can be expressed as<fnref n="6" />

$$ \dot{\mathbf{X}}(t) = -\nu\kappa(\mathbf{X}(t) - \mathbf{v}t)\vartheta(-t) + i\lambda\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} e^{i\mathbf{q}\cdot\mathbf{X}(t)} \phi_{\mathbf{q}}(t) + \boldsymbol{\xi}(t). \quad (14) $$

Note that for $\lambda=0$ and $t<0$, according to equations (12) and (14), each component of the particle position $\mathbf{X}$ and each mode of the field $\phi_q$ are described by non-interacting Ornstein–Uhlenbeck processes with relaxation times given, respectively, by

$$ \tau_{\bar{X}}^{-1} = \nu\kappa, \quad (15) $$

$$ \tau_{\phi}^{-1}(q) = \alpha_q \equiv D_{\alpha} q^{\alpha} (q^2 + r). \quad (16) $$

Accordingly, the field $\phi(\mathbf{x},t)$ exhibits not only long-range spatial correlations across the distance set by $\xi$, but also correlations over long times $\propto \xi^z$. In particular, for model A, the relaxation time of the long-wavelength (small-wavevector) modes diverges as the system reaches criticality ($r \to 0$). Such a divergence is present in model B both in the critical and off-critical case (i.e. finite $r$) due to the local conservation of the order parameter [51]. Consequently, the coupling between the tracer and the slow-relaxing field yields a non-Markovian effective dynamics of the particle, due to a generic lack of separation of timescales.

Note that, due to the simplicity of its formulation, this minimal model may not only describe a near-critical medium, but also serve as a simple approximation of various complex media exhibiting both spatial and temporal correlations.

### 3. Position of the particle

In order to derive the effective dynamics of the particle, we notice that the equation of motion (12) of the Gaussian field is linear in $\phi$, enabling its formal integration. This yields

$$ \phi_q(t) = \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t')} \left[ \lambda D_{\alpha} q^{\alpha} V_q e^{-i\mathbf{q}\cdot\mathbf{X}(t')} + \eta_q(t') \right], \quad (17) $$

where we assumed that the initial conditions $\phi_q(t_0)$ for the field configuration were imposed in the distant past, i.e. at $t_0 \to -\infty$, deeming them inconsequential. In order to circumvent possible issues of equilibration of critical systems, we assume that the system is taken adiabatically to criticality after it has reached a stationary state, such that the limit $\xi \to \infty$ is taken after $t_0 \to -\infty$. Notably, the field $\phi_q(t)$ depends dynamically only on the position of the particle $\mathbf{X}(t)$ and the configuration of the noise $\eta_q(t)$. Accordingly,

<footnote n="6">⁶ Note that, in the case of model A ($\alpha=0$), the average field mode $\phi_q=0$ may take arbitrarily large absolute values as criticality is approached for $r \to 0$. However, this fact is inconsequential for the particle dynamics: indeed, the mode with $\mathbf{q}=0$ does not contribute to equation (14). In order to counteract this growth, one would need to add a suitable chemical potential—e.g. by replacing $\mathcal{H}_{\phi} \mapsto \mathcal{H}_{\phi} + \lambda V_{q=0} \int d\mathbf{x} \phi(\mathbf{x})$ in equation (3).</footnote>
the expression (17) for the dynamics of the field can be substituted into equation (14), describing the position of the particle, rendering the effective dynamics

$$ \dot{\mathbf{X}}(t) + \nu\kappa(\mathbf{X}(t) - \mathbf{v}t)\vartheta(-t) = \mathbf{F}[\{\mathbf{X}(t)\}] + \boldsymbol{\xi}(t), \quad (18) $$

where

$$ \mathbf{F}[\{\mathbf{X}(t)\}] = i\lambda\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} e^{i\mathbf{q}\cdot\mathbf{X}(t)} \left[ \phi_{\mathbf{q}}^{(0)}(t) + \lambda D q^{\alpha} V_{\mathbf{q}} \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t')} e^{-i\mathbf{q}\cdot\mathbf{X}(t')} \right], \quad (19) $$

and where

$$ \phi_{\mathbf{q}}^{(0)}(t) = \int_{-\infty}^{t} dt' e^{-\alpha_{\mathbf{q}}(t-t')} \eta_{\mathbf{q}}(t') \quad (20) $$

is a colored Gaussian noise. Note that the effective dynamics is not only non-Markovian, but also non-linear, and thus cannot be solved analytically. To overcome this difficulty, we will resort to the perturbative expansion introduced in section 3.1.

### 3.1. Perturbative expansion (weak-coupling approximation)

In order to obtain analytical predictions for the effects of the field-particle interaction on the particle dynamics, following [18, 29, 32, 52], we assume that $\lambda$ is a small parameter in terms of which one is allowed to perform a perturbative expansion:

$$ \mathbf{X}(t) = \mathbf{X}^{(0)}(t) + \lambda \mathbf{X}^{(1)}(t) + \lambda^2 \mathbf{X}^{(2)}(t) + \mathcal{O}(\lambda^3), \quad (21) $$

$$ \phi_{\mathbf{q}}(t) = \phi_{\mathbf{q}}^{(0)}(t) + \lambda \phi_{\mathbf{q}}^{(1)}(t) + \lambda^2 \phi_{\mathbf{q}}^{(2)}(t) + \mathcal{O}(\lambda^3). \quad (22) $$

In the studies cited above, the expansion at the second order in $\lambda$ already proved sufficient to capture non-trivial effects on the particle motion. The results obtained within this weak-coupling approximation will be compared with those of numerical simulations of the system in cf figure 2.

When the field and the particle are decoupled (i.e. for $\lambda=0$), the position of the latter is given by (see equations (12) and (14))

$$ \mathbf{X}^{(0)}(t) = \begin{cases} \mathbf{v}(t - \frac{1}{\nu\kappa}) + \int_{-\infty}^{t} dt' e^{-\nu\kappa(t-t')} \boldsymbol{\xi}(t'), & \text{for } t \le 0, \\ \mathbf{X}^{(0)}(0) + \int_{0}^{t} dt' \boldsymbol{\xi}(t'), & \text{for } t > 0, \end{cases} \quad (23) $$

where we chose the coordinate system such that at $t=0$ the minimum of the trap is located at the origin. In the expression above, for $t \le 0$, one identifies the time integral of the noise as the Ornstein-Uhlenbeck process $\mathbf{O}(t)$ with the relaxation time $\tau_X$ reported in equation (15), and for $t > 0$, as the Wiener process $\mathbf{W}(t)$ with diffusion coefficient $\nu T$. The symmetry of the particle dynamics reported below equation (8) implies that $\langle \mathbf{X}(t) \rangle$ is actually an even function of $\lambda$. Accordingly, the second-order correction to the
decoupled dynamics is the lowest-order one with a non-vanishing mean, and it is given by

$$ \dot{\mathbf{X}}^{(2)}(t) + \vartheta(-t) \nu \kappa \mathbf{X}^{(2)}(t) = i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} e^{i\mathbf{q}\cdot\mathbf{X}^{(0)}(t)} V_{-\mathbf{q}} \left[ \phi_{\mathbf{q}}^{(1)}(t) + i\mathbf{q} \cdot \mathbf{X}^{(1)}(t) \phi_{\mathbf{q}}^{(0)}(t) \right]. \quad (24) $$

Here, we recognize the solution $\phi_{\mathbf{q}}^{(0)}(t)$ of the field dynamics with $\lambda=0$ as the colored noise introduced in equation (20).

## 3.2. Steady-state average particle position

Using the expression for $\langle e^{i\mathbf{q}\cdot(\mathbf{O}(t)-\mathbf{O}(t'))} \rangle$ reported in equation (A.7), where $\mathbf{O}(t)$ is the stationary Ornstein–Uhlenbeck process at time $t$ (see e.g. [18]), we calculate the lowest-order correction to the average position of the particle while it is subject to the harmonic confinement. The details of the computation are provided in appendix A, while we report here only the final expression

$$ \langle \lambda^2 \mathbf{X}^{(2)} \rangle = -\hat{\mathbf{e}}_1 \lambda^2 \frac{v}{\kappa} \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + r} |V_{\mathbf{q}}|^2 \int_0^{\infty} ds e^{-(\alpha_q - i q_1 v) s - (T/\kappa) q^2 (1 - e^{-\kappa v s})}, \quad (25) $$

where, for the sake of concreteness, we assumed $\mathbf{v} = \hat{\mathbf{e}}_1 v$, with $\hat{\mathbf{e}}_1$ being the unit vector in the direction of $x_1$. This time-independent average shift of the position of the particle from the center of the moving trap caused by the interaction with the field was derived previously in [32]. Note that this coupling actually induces an additional friction, since the displacement in equation (25) is negative, i.e. it occurs in the direction opposite to the velocity of the trap, as we prove in appendix A.1. An analogous phenomenon was reported in [53] while discussing a tracer driven in a colloidal bath. This friction can be understood as a result of the force in equation (8) pulling the particle along the gradient of the steady-state field configuration forming behind it. A detailed description of this steady-state field profile is presented in section 4 below.

In order to better understand the qualitative behavior of the average particle displacement in equation (25), we first consider the noiseless limit $T \to 0$ in which the fluctuations due to the thermal bath are negligible, and the right-hand side of the equation takes the simpler form

$$ \langle \lambda^2 \mathbf{X}^{(2)} \rangle_{T=0} = -\hat{\mathbf{e}}_1 \lambda^2 \frac{v}{\kappa} \int \frac{d^d \mathbf{q}}{(2\pi)^d} D_{\alpha} q^{\alpha} |V_{\mathbf{q}}|^2 \frac{q_1^2}{\alpha_q^2 + q_1^2 v^2}. \quad (26) $$

Accordingly, as a function of the velocity $v$, one observes first a linear growth of $\langle \lambda^2 \mathbf{X}^{(2)} \rangle_{T=0}$ followed by a decay to zero $\propto v^{-1}$ as $v$ is further increased$^7$. Accordingly, the shift vanishes both as $v \to 0$ and as $v \to \infty$. The former limit is easy to understand, because when the trap is not moving, the dynamics of the system is invariant under the inversion $x_1 \mapsto -x_1$, and thus the shift has to vanish. In the latter limit, the particle

<footnote>$^7$ In fact, provided that the $\mathbf{q}$-integral in equation (26) converges, the right-hand side of this equation behaves like $v/(1 + Av^2)$, where $A$ is a constant—for further discussion, see the supplemental material of [32].</footnote>
moves through the medium too quickly for the field to react to its presence, and thus the dynamics of the field effectively decouples from that of the particle. A detailed analysis of the case of dissipative (model A) dynamics in spatial dimension $d=1$ is presented in the supplemental material of [32].

### 3.3. Recoil

The main focus of this work is on the motion of the particle after it has been released from the harmonic trap. It turns out that the particle exhibits recoil, i.e. it moves in the direction opposite to the drag, in a manner similar to what is observed experimentally in viscoelastic media [8, 19–21]. We note that this backwards motion arises as a result of the field-particle interaction, and vanishes for $\lambda=0$: indeed, due to the symmetry of the particle’s dynamics reported under equation (8), one can verify that

$$ \langle \mathbf{X}(t) \rangle - \langle \mathbf{X}(0) \rangle = \langle \lambda^2 \Delta \mathbf{X}^{(2)}(t) \rangle + \mathcal{O}(\lambda^4), \quad (27) $$

where we denoted the time-dependent recoil by $\lambda^2 \Delta \mathbf{X}^{(2)}(t)$. We postpone the derivation of the recoil to appendix A and present here the final expression for its average:

$$ \begin{aligned} \langle \lambda^2 \Delta \mathbf{X}^{(2)}(t) \rangle = & -\hat{\mathbf{e}}_1 \lambda^2 \nu v \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + r} |V_q|^2 \frac{1 - e^{-(\alpha_q + q^2 \nu T)t}}{\alpha_q + q^2 \nu T} \\ & \times \int_0^\infty ds \, e^{-(\alpha_q - i q_1 v)s - (T/\kappa)q^2(1 - e^{-\kappa \nu s})}. \end{aligned} \quad (28) $$

Note that equation (28), as a function of time $t$, possesses the mathematical structure of a weighted exponential decay towards its final value. However, when the relaxation time $\tau_\phi(q)$ of the small-wavevector modes $\phi_{q \simeq 0}$ diverges (see equation (16)), i.e. generically for model B and for model A at criticality, algebraic decays are expected. Note also that a non-vanishing temperature $T$ causes a renormalization of the (inverse) relaxation times $\alpha_q$ according to

$$ \alpha_q \mapsto \alpha'_q = \alpha_q + q^2 \nu T, \quad (29) $$

which is equivalent to $D_0 \mapsto D_0 + \nu T$ for model A, and $r \mapsto r + \nu T/D_2$ for model B. Consequently, a finite temperature $T$ suppresses critical fluctuations in the case of model B, as it effectively drives the system away from the critical point $r=0$. The short-time behavior of the recoil can be easily obtained from equation (28) by approximating $1 - e^{-\alpha'_t t} \simeq \alpha'_q t$, which renders

$$ \langle \lambda^2 \Delta \mathbf{X}^{(2)}(t) \rangle \simeq -\hat{\mathbf{e}}_1 \lambda^2 \nu v t \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + r} |V_q|^2 \int_0^\infty ds \, e^{-(\alpha_q - i q_1 v)s - (T/\kappa)q^2(1 - e^{-\kappa \nu s})}. \quad (30) $$

In passing, we note that this expression coincides with equation (25) for the steady-state shift $\langle \lambda^2 \mathbf{X}^{(2)}(t) \rangle$ multiplied by the dimensionless time $t/\tau_X$, where $\tau_X$ is the relaxation time of the particle in the trap, introduced in equation (15). This is because, before the
particle is released, the field-induced friction experienced by the particle is counterac-
ted by the force $-\mathbf{f}_0 = -\kappa \langle \lambda^2 \mathbf{X}^{(2)}(t < 0) \rangle$ generated by the trap as a consequence of the
additional shift $\langle \lambda^2 \mathbf{X}^{(2)}(t < 0) \rangle$ from its center, maintaining mechanical equilibrium in
the co-moving reference frame in which the trap is at rest. When the particle is released
from confinement, the force $\mathbf{f}_0$ is exerted on the particle by the field in its initial con-
figuration at $t = 0$. This force, due to the overdamped nature of the particle dynamics
characterized by the mobility $\nu$, gives an initial velocity $\mathbf{v}_0 = \nu \mathbf{f}_0$ to the particle, caus-
ing a recoil which initially grows as $\mathbf{v}_0 t = \nu \kappa t \langle \lambda^2 \mathbf{X}^{(2)}(t < 0) \rangle = (t/\tau_X) \langle \lambda^2 \mathbf{X}^{(2)}(t < 0) \rangle$; see
equation (15), as noted above. In the noiseless limit $T \to 0$, the dynamics is deterministic
and equation (28) becomes

$$
\langle \lambda^2 \Delta \mathbf{X}^{(2)}(t) \rangle = -\hat{\mathbf{e}}_1 \lambda^2 \nu v \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + r} |V_{\mathbf{q}}|^2 \frac{1 - e^{-\alpha_q t}}{\alpha_q^2 + q_1^2 v^2}, \quad (31)
$$

which does not depend on the trap strength $\kappa$ (which controls the spatial extent of the
thermal fluctuations of the particle position while it is driven by the trap). The corres-
ponding recoil is plotted as a function of time in figure 2 for the non-critical model A
in $d=1$ and with the Gaussian interaction kernel given in equation (5), along with the
results of a numerical simulation of the system. For the data presented, the relative
difference between the analytic prediction and the results of numerical simulations is
of the order of 4%. This difference increases upon increasing the coupling strength $\lambda$,
reaching $\approx 6\%$ for $\lambda = 0.5$.

It is convenient to express the recoil in the noiseless limit, see equation (31), in terms
of a scaling function $\mathcal{F}^{A/B}$, i.e.

$$
\langle \lambda^2 \Delta \mathbf{X}^{(2)}(t) \rangle = -\hat{\mathbf{e}}_1 g \frac{v \tau_R^2}{\tau_X} \mathcal{F}^{A/B} \left( \frac{R}{\xi}, \frac{t}{\tau_R}, \frac{v \tau_R}{R} \right), \quad (32)
$$

where

$$
\mathcal{F}^{A/B} \left( 1/\tilde{\xi}, \tilde{t}, \tilde{v} \right) = \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + 1/\tilde{\xi}^2} |\hat{V}_{\mathbf{q}}|^2 \frac{1 - e^{-\tilde{t}q^\alpha (q^2 + 1/\tilde{\xi}^2)}}{q^{2\alpha} \left(q^2 + 1/\tilde{\xi}^2\right)^2 + \tilde{v}^2 q_1^2}, \quad (33)
$$

and where $\hat{V}_{\mathbf{q}}$ is the Fourier transform (see equation (11)) of the interaction kernel
rescaled by the particle's effective radius introduced after equation (4). The superscript
A or B to $\mathcal{F}$ and to similar scaling functions introduced further below indicates that
they refer to the case of model A ($\alpha=0$) or model B ($\alpha=2$). The prefactor $g\nu\tau_R^2/\tau_X$ in
equation (32) has the physical dimensions of a length, and it is expressed in terms of a
dimensionless coupling constant

$$
g = \frac{\lambda^2}{\kappa R^d}, \tag{34}
$$

of the timescale $\tau_X$ defined in equation (15), and of

$$
\tau_R = R^z / D_\alpha, \tag{35}
$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Dependence of the recoil amplitude $|\Delta X|$ on time, for the non-critical model A in spatial dimension $d=1$, in the absence of noise ($T=0$) and for various values of the correlation length $\xi$ of the field. The interaction of the particle with the field is assumed to be Gaussian as in equation (5) with $R=1$, while the other parameters of the model are $\kappa=0.5$, $D=1$, $\nu=1$, $\lambda=0.1$, and $v=1$. The solid lines correspond to the predictions of the perturbative calculation reported in equation (31), while the dashed lines are the result of a numerical simulation in which the particle was dragged through the medium for $t \in [0, 250]$, before being released from its confinement at time $t=250$.</figcaption></figure>

which characterizes the relaxation time of the critical field over a distance of the order of the particle size $R$ (its definition follows from equation (16) with $r=0$ and $q \sim 1/R$). In the convention introduced in equation (32), distances are measured in terms of $R$, and times in terms of $\tau_R$, allowing for the introduction of the dimensionless correlation length and the dimensionless time according to

$$ \tilde{\xi} = \xi/R \quad \text{and} \quad \tilde{t} = t/\tau_R. \qquad (36) $$

We note that the last argument of the scaling function in equation (32), i.e. the dimensionless velocity

$$ \tilde{v} \equiv v\tau_R/R = 2\text{Wi} = (R/l_{v,\alpha})^{\alpha+1}, \qquad (37) $$

can be expressed in terms of the Weissenberg number Wi used in the context of micro-rheological experiments in viscoelastic media [15, 54], which compares the typical times describing the medium relaxation and its deformation. Indeed, the dragged particle deforms the medium over the length of its diameter in a time $2R/v$. Note that the length scale

$$ l_{v,\alpha} = (D_{\alpha}/v)^{1/(\alpha+1)} \qquad (38) $$
<table><caption>**Table 1.** Behavior of the scaling function $\mathcal{F}^{A/B}(1/\tilde{\xi}, \tilde{t}, \tilde{v})$ defined in equation (32) for $\tilde{v} \to 0$. The prefactors $f_{0,2}(d; \tilde{t})$, $g_{0,2}(d)$, and $h(d; \tilde{\xi})$, which depend on the dimensionless scaling variables $\tilde{t}$ and $\tilde{v}$ defined in equations (36) and (37), are determined by means of the Mellin transform in equations (B.5) and (B.6). Physically, the generic non-analytic dependence of the scaling function on $\tilde{v}$ is due to the fact that the relaxation time of the field, whose inverse is given by equation (16), diverges. The dependence on $d$ of the power of $\tilde{v}$ appearing in the table originates from the infrared divergence of the q-integral in equation (33), which is regularized for sufficiently large $d$.</caption><thead><tr><th></th><th>$\mathcal{F}^{A/B}(1/\tilde{\xi}, \tilde{t}, \tilde{v} \to 0)$</th><th>$\tilde{\xi} < \infty$</th><th>$\tilde{\xi} \to \infty$</th></tr></thead><tbody><tr><td rowspan="2">$\alpha = 0$</td><td>$\tilde{t} < +\infty$</td><td>$\sim 1$</td><td>$\sim 1 + f_0(d; \tilde{t}) \tilde{v}^{d-2}$</td></tr><tr><td>$\tilde{t} \to \infty$</td><td>$\sim 1$</td><td>$\sim 1 + g_0(d) \tilde{v}^{d-4}$</td></tr><tr><td rowspan="2">$\alpha = 2$</td><td>$\tilde{t} < +\infty$</td><td>$\sim 1$</td><td>$\sim 1 + f_2(d; \tilde{t}) \tilde{v}^{(d-4)/3}$</td></tr><tr><td>$\tilde{t} \to \infty$</td><td>$\sim 1 + h(d; \tilde{\xi}) \tilde{v}^{d-2}$</td><td>$\sim 1 + g_2(d) \tilde{v}^{(d-8)/3}$</td></tr></tbody></table>

arises naturally while describing the steady-state field configuration at criticality, as discussed in cf section 4.2.

We note that in the limit of instantaneous equilibration of the field, corresponding (for non-critical model A) to $D_0 \to \infty$ or, equivalently, $\tau_R \to 0$, the expression for the recoil in equations (32) and (33) generically vanishes; in this *adiabatic* limit, the timescale separation is recovered<fnref n="8" />. We therefore emphasize that the presence of the recoil stems from the fact that the dynamics of the particle and of the field occur over comparable timescales.

The convention adopted here allows us to investigate the following relevant limits of equation (32): (i) the limit of small velocities $v$; (ii) the long-time limit which we denote by

$$ \mathcal{R}^{A/B} (1/\tilde{\xi}, \tilde{v}) = \lim_{\tilde{t} \to \infty} \mathcal{F}^{A/B} (1/\tilde{\xi}, \tilde{t}, \tilde{v}), \qquad (39) $$

where we will henceforth refer to the quantity $\tilde{v} \times \mathcal{R}^{A/B}(1/\tilde{\xi}, \tilde{v})$ as the *dimensionless recoil range* because it describes the final particle position (recall that $\langle \lambda^2 \Delta \mathbf{X}^{(2)} \rangle \propto v \mathcal{F}^{A/B}$ in equation (32)), and (iii) the critical limit $\xi \to \infty$. Concerning (i), in the case of finite $\xi$ and $t$, at small velocities, the recoil is characterized by a linear growth upon increasing $v$, as indicated by the prefactor on the right-hand side of equation (32). However, when at least one of the relevant tunable scales in equation (32) diverges (i.e. either $t$ or $\xi$), the scaling function $\mathcal{F}^{A/B}$ acquires a non-trivial behavior at small $v$, which is analyzed by means of the Mellin transform in appendix B.1 and is reported in table 1.

<footnote n="8">⁸ More precisely, timescale separation requires that the relaxation time of the slowest mode of the field (i.e. the one corresponding to the smallest wavevector, see equation (16)) vanishes. Hence, in the case of model B and critical model A, the presence of zero modes with diverging relaxation time precludes a bona-fide *adiabatic* limit (see [30, 52] for further discussions).</footnote>
Depending on the type of dynamics and on the space dimensionality $d$, this might even result in a diverging recoil range as $v \to 0$: this occurs for critical model A in $d < 3$ and for critical model B in $d < 5$ (recall the presence of $v$ in the prefactor in equation (32)). For critical model A in $d=3$ and model B in $d=1$ and $d=5$ with $r > 0$ and $r=0$, respectively, the limit $v \to 0$ yields a finite recoil range. The physical interpretation of this result is that, in a system characterized by diverging relaxation time, any small deformation rate induced by the driving velocity of the particle is still much larger than the relaxation rate of the medium. Hence, for any infinitesimal, yet nonzero $\tilde{v}$, the recoil range may attain a finite value or even diverge. A similar behavior for $\langle \lambda^2 X^{(2)}(t < 0) \rangle$ given by equation (25) was reported in the supplemental material of [32].

In order to discuss in more detail the asymptotic long-time behavior of the recoil, which is the object of the next section, the kernel $V(x)$ describing the field-particle interaction is henceforth assumed to have the Gaussian form reported in equation (5), unless stated otherwise.

### 3.4. Long-time behavior

In this section we analyze the behavior of the scaling function of the recoil defined in equation (33), starting with the long-time asymptotics of the critical recoil (i.e. the average particle recoil when $\tilde{\xi} \to \infty$), and following with the description of $\mathcal{R}^{A/B}$ defined in equation (39) for a large but finite correlation length. The corresponding calculations are outlined in appendix B.2. Moreover, explicit expressions for the critical recoil in the case of model A dynamics with $d \in \{1, 2, 3\}$ are reported in appendix B.3. Here we summarize the key results of our analysis.

For $d > 2$, at long times the critical recoil satisfies

$$ \mathcal{F}^{A/B}(0, \tilde{t}, \tilde{v}) \sim \mathcal{R}^{A/B}(0, \tilde{v}) + \frac{\tilde{t}^{(2-d)/z}}{\tilde{v}^2} |\hat{V}_0|^2 \mathcal{C}_\alpha(d), \quad \text{for } \tilde{t} \gg 1, \qquad (40) $$

where the dynamical critical exponent $z$ of the field was introduced after equation (9), and a finite dimensionless velocity $\tilde{v}$ is assumed. Here, $\mathcal{C}_\alpha(d)$ is a constant determined solely by the spatial dimensionality and by the dynamics of the field, whose expression is reported in equation (B.9). Interestingly, $\mathcal{R}^{A/B}(0, \tilde{v})$ does not exist for $d \le 2$ due to infrared divergences: physically, this corresponds to the particle moving backwards indefinitely due to diverging relaxation times. The presence of such a diverging recoil range for $d \le 2$ can also be understood as a consequence of a peculiarity of the scalar Gaussian field theory at the critical point, discussed in more detail in section 4.2 (see cf equation (60)). Notably, for $d > 2$, $\mathcal{C}_\alpha(d)$ is negative and the recoil converges to its final value $\mathcal{R}^{A/B}$ from below, with a power-law dependence on time with exponent $(2-d)/z$. However (as shown in appendix B.3, where the full expressions for the scaling function in the case of critical dissipative dynamics with $d \in \{1, 2, 3\}$ are provided), in $d=1$, the long-time asymptotic behavior of the recoil is governed by the second term on the right-hand side of equation (40), i.e. it grows linearly as a function of time with the corresponding prefactor $\propto \mathcal{C}_0(1)$ which is positive. In $d=2$, instead, both $\mathcal{C}_0(d)$ and $\mathcal{C}_2(d)$ have simple poles as functions of $d$, and thus the recoil increases logarithmically
<table><caption>**Table 2.** Long-time behavior of the scaling function $\mathcal{F}^{A/B}(1/\tilde{\xi}, \tilde{t}, \tilde{v})$ in equation (33). The divergence of the recoil range for $\tilde{t} \to \infty$, occurring if $\tilde{\xi} \to \infty$ and $d \le 2$, indicates an endless backward motion of the particle after its release from the trap (see the discussion in sections 4.2 and 3.4). Note that $C_{0,2}(d)$ is positive for $d < 2$, and negative for $d > 2$.</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>$\mathcal{F}^{A/B}(1/\tilde{\xi}, \tilde{t} \gg 1, \tilde{v})$</td><td>$\tilde{\xi} < \infty$</td><td>$\tilde{\xi} \to \infty$</td><td></td></tr><tr><td>$\alpha = 0$</td><td>$d \le 2$<br>$d > 2$</td><td>$\sim \mathcal{R}^A(1/\tilde{\xi}, \tilde{v}) - \frac{|\hat{V}_0|^2}{2(4\pi)^{d/2}} \frac{\tilde{\xi}^6 e^{-i/\tilde{\xi}^2}}{\tilde{t}^{d/2+1}}$</td><td>$\sim \frac{\tilde{t}^{(2-d)/2} |\hat{V}_0|^2}{\tilde{v}^2} C_0(d)$<br>$\sim \mathcal{R}^A(0, \tilde{v}) + \frac{|\hat{V}_0|^2}{\tilde{t}^{(d-2)/2} \tilde{v}^2} C_0(d)$</td></tr><tr><td>$\alpha = 2$</td><td>$d \le 2$<br>$d > 2$</td><td>$\sim \mathcal{R}^B(1/\tilde{\xi}, \tilde{v}) - \frac{|\hat{V}_0|^2}{(4\pi)^{d/2}} \frac{\tilde{\xi}^{d+2}}{\tilde{v}^2 \tilde{t}^{d/2}}$</td><td>$\sim \frac{\tilde{t}^{(2-d)/4} |\hat{V}_0|^2}{\tilde{v}^2} C_2(d)$<br>$\sim \mathcal{R}^B(0, \tilde{v}) + \frac{|\hat{V}_0|^2}{\tilde{t}^{(d-2)/4} \tilde{v}^2} C_2(d)$</td></tr></tbody></table>

upon increasing time. We emphasize that the second term on the right-hand side of equation (40) is largely universal, as it depends only on the behavior of the interaction kernel $\hat{V}_q$ at wavevector $q=0$, rendering it insensitive to the details of the field-particle interaction (as long as it is linear, translationally and rotationally invariant).

Let us now turn to non-critical systems, where $\mathcal{R}^{A/B}$ turns out to be always finite. In particular, in the case of dissipative dynamics, the particle relaxes exponentially to its final position according to

$$ \mathcal{F}^A(1/\tilde{\xi}, \tilde{t}, \tilde{v}) \sim \mathcal{R}^A(1/\tilde{\xi}, \tilde{v}) - \frac{|\hat{V}_0|^2}{2(4\pi)^{d/2}} \frac{\tilde{\xi}^6 e^{-\tilde{t}/\tilde{\xi}^2}}{\tilde{t}^{d/2+1}}, \quad \text{for } \tilde{t} \gg 1. \quad (41) $$

In the case of model B, instead, the particle recoil attains its maximal value following an algebraic approach according to

$$ \mathcal{F}^B(1/\tilde{\xi}, \tilde{t}, \tilde{v}) \sim \mathcal{R}^B(1/\tilde{\xi}, \tilde{v}) - \frac{|\hat{V}_0|^2}{(4\pi)^{d/2}} \frac{\tilde{\xi}^{d+2}}{\tilde{v}^2 \tilde{t}^{d/2}}, \quad \text{for } \tilde{t} \gg 1. \quad (42) $$

A summary of the long-time behavior of the scaling function $\mathcal{F}^{A/B}$ is presented in table 2. We point out that when the dynamics is characterized by diverging relaxation times, and thus the approach to the maximal recoil occurs algebraically in time, the difference $\mathcal{F}^{A/B} - \mathcal{R}^{A/B}$ between the corresponding scaling functions scales with the dimensionless velocity $\tilde{v}$ like $\sim \tilde{v}^{-2}$, see equations (40) and (42). In the case of non-critical model A (which is characterized by a finite relaxation time), instead, the corresponding difference in equation (41) does not depend on the velocity. These two dependences on the velocity coincide with the behavior of $\mathcal{F}^{A/B}$ in the regime of $\tilde{v} \gg 1$ (far from equilibrium) and $\tilde{v} \ll 1$ (close to equilibrium), respectively. $\mathcal{F}^{A/B}$ for $d=3$ as a function of time are plotted in figures 3(a) and 4(a) for model A and B, respectively. The asymptotic curves introduced in equations (40) and (42) are reported as dashed lines.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Scaling function $\mathcal{F}^A$ of the recoil for model A in spatial dimension $d=3$, see equation (33). Panel (a) shows $\mathcal{F}^A$ as a function of the dimensionless time $\tilde{t}$ with $\tilde{v}=1$, and various values of the dimensionless correlation length $\xi/R$. The dashed line corresponds to the long-time behavior at the critical point reported in equation (40), where $z=2$. Panel (b) shows the dependence of the long-time limit $\mathcal{R}^A$ of the scaling function $\mathcal{F}^A$ of the recoil (see equation (39)) as a function of $\xi/R$, for various values of the scaled driving velocity $\tilde{v}$. The behaviors of these curves for large values of $\xi/R$, given by equation (43), are reported as dashed lines.</figcaption></figure>

For $d<2$, it is instructive to investigate the asymptotic behavior of $\mathcal{R}^{A/B}$ for large but finite values of the correlation length $\xi$. It turns out that the expressions corresponding to model A and B actually coincide, and that they are given by

$$ \mathcal{R}^{A/B} (1/\tilde{\xi}, \tilde{v}) \sim |\hat{V}_0|^2 \frac{\Gamma(1-d/2) \tilde{\xi}^{2-d}}{(4\pi)^{d/2}} \tilde{v}^2, \quad \text{for} \quad \tilde{\xi} \gg 1. \qquad (43) $$

Recall that the definition of the dimensionless velocity $\tilde{v}$ in equation (37) differs for models A and B. Note that the exponent that describes the power-law divergence of the recoil upon increasing the correlation length $\xi$ in equation (43) agrees with that reported in equation (40), which describes the long-time asymptotics of the critical recoil, taking into account the scaling relation $t \sim \xi^z$, where $z$ is the dynamical critical exponent. It turns out that, for $d=3$, equation (43) describes the asymptotic $\sim 1/\xi$ convergence of $\mathcal{R}^{A/B}$ towards its maximal value obtained in the limit of $\xi \to \infty$, which we show by plotting it as dashed curves in figures 3(b) and 4(b).

In passing, we note that while describing the final position of the particle after the recoil, the limits $v \to 0$ and $\xi \to \infty$ do not commute, see equation (43). In fact, it turns out that the final particle position depends on the driving velocity via the variable

$$ \tilde{v} \tilde{\xi}^{z-1} = v \xi^{z-1} / D_{\alpha}, \qquad (44) $$

which effectively quantifies the extent to which the medium (whose relaxation rate is controlled by $\xi$, see equation (16)) is driven out of equilibrium. Hence, the asymptotic behavior for $\tilde{\xi} \gg 1$ is characterized by a $\tilde{v}^{-2}$ scaling, which is the same as that of the scaling function $\mathcal{F}^{A/B}$ in the regime of large driving velocity—see the denominator of the integral in equation (33).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Scaling function $\mathcal{F}^B$ of the recoil for model B in spatial dimension $d=3$, see equation (33). Panel (a) shows $\mathcal{F}^B$ as a function of the dimensionless time $\tilde{t}$ with $\tilde{v}=1$, and various values of the dimensionless correlation length $\xi/R$. The dashed brown line corresponds to the long-time behavior at the critical point reported in equation (40), where $z=4$, and the dashed purple line corresponds to the long-time behavior in the non-critical case reported in equation (42). Panel (b) shows the dependence of the long-time limit $\mathcal{R}^B$ of the scaling function $\mathcal{F}^B$ of the recoil (see equation (39)) as a function of $\xi/R$, for various values of the scaled driving velocity $\tilde{v}$. The behaviors of these curves for large values of $\xi/R$, given by equation (43), are reported as dashed lines.</figcaption></figure>

The dimensionless recoil range $\tilde{v} \times \mathcal{R}^{A/B}$ is plotted as a function of the dimensionless velocity $\tilde{v}$ for $d=3$ in figure 5(a) for model A. In particular, one observes that it vanishes at $\tilde{v}=0$ as long as the correlation length $\xi$ remains finite, while it reaches a finite value as $\xi \to \infty$, according to table 1. Note that the value $\tilde{v}_{max}$ of the velocity $\tilde{v}$ at which $\tilde{v} \times \mathcal{R}^A$ is maximal decreases as $\tilde{\xi}$ increases. In figure 5(b) we plot $\tilde{v} \times \mathcal{R}^B(\tilde{v})$ with solid lines; additionally, we plot $\tilde{v} \times \mathcal{R}^A(\tilde{v})$ with dashed lines to demonstrate that the dimensionless recoil range of the particle differs considerably between the cases of model A and model B dynamics for $\tilde{v} \ll 1$. We recall that $\tilde{v}$—defined in equation (37)—describes the relative significance of the effects of the deformation compared with the relaxation of the field. Given that the expressions for the relaxation times depend on the dynamics of the field, as manifested by their $\alpha$-dependence in equation (16), we expect the final position of the recoil to differ considerably in the relaxation-dominated ($\tilde{v} \ll 1$) regime. In fact, the long-time limit of the scaling function introduced in equation (39) in the limit of $v \to 0$ takes the form

$$ \mathcal{R}^{A/B} (1/\tilde{\xi}, \tilde{v}) \xrightarrow{\tilde{v} \to 0} \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^{2\alpha} (q^2 + 1/\tilde{\xi}^2)^3} |\hat{V}_{\mathbf{q}}|^2, \qquad (45) $$

which manifestly depends on the dynamics considered (via the factor $q^{2\alpha}$ in the denominator<fnref n="9" />). This limit of weak driving out of equilibrium may be regarded as a linear

<footnote n="9">⁹ We recall that the integral in equation (45) may possibly exhibit infrared divergences, as already discussed after equation (40), and in table 1.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Dimensionless recoil range $\tilde{v} \times \mathcal{R}^{A/B}$ (see equations (32) and (39)) as a function of the dimensionless driving velocity $\tilde{v}$ in spatial dimension $d=3$, for (a) model A or (b) model B, and various values of the correlation length $\xi$. In order to ease the comparison between these two dynamics, the dashed lines in panel (b) correspond to the behavior of the recoil range for model A at small or large velocities. As indicated in the inset, for large velocities, the curves corresponding to models A and B coincide.</figcaption></figure>

response description of the recoil—indeed, using equations (32) and (39), one deduces that $\langle \Delta\mathbf{X}^{(2)}(t \to \infty) \rangle \propto v\mathcal{R}^{A/B}(1/\tilde{\xi}, \tilde{v} = 0)$ for $v \to 0$, i.e. a linear dependence on $v$.

Finally, in the opposite limit $\tilde{v} \gg 1$, the curves corresponding to model A and model B collapse onto each other (this is better visible in the inset of figure 5(b)). Indeed, for large $\tilde{v}$, the expression for $\mathcal{R}^{A/B}$ (see equations (33) and (39)) becomes $\alpha$-independent:

$$ \mathcal{R}^{A/B}(1/\tilde{\xi}, \tilde{v}) \sim \frac{1}{\tilde{v}^2} \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{1}{q^2 + 1/\tilde{\xi}^2} |\hat{V}_{\mathbf{q}}|^2, \quad \text{for} \quad \tilde{v} \gg 1. \quad (46) $$

## 4. Non-equilibrium steady-state configuration of the field

Having discussed the properties of the recoil in the previous sections, below we turn our attention to the slow-relaxing field, since the force acting on the tracer is related to its gradient via equation (8). In particular, we shall focus on the average field profile in the non-equilibrium steady state that the system reaches at long times in the co-moving frame when the tracer is driven by the moving trap. The discussion below follows the one presented in [31], and builds upon the results reported therein and in [32]. Beyond these references, we present here a detailed characterization of the steady-state configuration of the field for model A at criticality and for model B both at the critical point and in its vicinity. Such a characterization provides insight into the mechanism for the divergence of the recoil range and the algebraic asymptotic long-time behavior of the scaling function reported in equation (40).
We focus on the part of the protocol during which the particle is driven, i.e. on $t \le 0$, and measure the particle position $\mathbf{X}$ from the minimum of the moving trap by defining $\mathbf{Z} = \mathbf{X} - \mathbf{v}t$. It is also convenient to introduce the translated field $\varphi(\mathbf{x}, t) = \phi(\mathbf{x} + \mathbf{v}t, t)$. Note that both $\mathcal{H}_\phi$ and $\mathcal{H}_{\text{int}}$ given by equations (3) and (4) are invariant under spatial translation, i.e. $\mathcal{H}_\phi[\phi] = \mathcal{H}_\phi[\varphi]$ and $\mathcal{H}_{\text{int}}[\phi, \mathbf{X}] = \mathcal{H}_{\text{int}}[\varphi, \mathbf{Z}]$. Moreover, the noise acting on the particle and the field is homogeneous in space. Hence, the equations of motion (12) and (14) can be expressed in terms of $\mathbf{Z}$ and $\varphi_q$ as [31]

$$ \dot{\mathbf{Z}} = -\mathbf{v} - \kappa\nu\mathbf{Z} + \lambda\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} i\mathbf{q} V_{-\mathbf{q}}\varphi_q e^{i\mathbf{q}\cdot\mathbf{Z}} + \boldsymbol{\xi}, \quad (47) $$

$$ [\partial_t + \alpha_q - i\mathbf{q} \cdot \mathbf{v}] \varphi_q(t) = \lambda D_\alpha q^\alpha V_q e^{-i\mathbf{q} \cdot \mathbf{Z}} + \eta_q, \quad (48) $$

where $\alpha_q$ is defined in equation (16), and $\varphi_q(t)$ is the Fourier transform of $\varphi(\mathbf{x}, t)$ (according to the definitions in equation (11)). Note that changing the frame of reference essentially amounts to changing $\alpha_q \mapsto (\alpha_q - i\mathbf{q} \cdot \mathbf{v})$ on the left-hand side of the equation of motion of the field.

For a particle dragged through the medium for a sufficiently long time, the system is expected to reach a stationary state characterized by $\langle \dot{\mathbf{Z}} \rangle_{\text{ss}} = 0$ and $\langle \partial_t \varphi_q \rangle_{\text{ss}} = 0$, which is satisfied for (see equations (47) and (48))

$$ \langle \mathbf{Z} \rangle_{\text{ss}} = -\mathbf{v} / (\kappa\nu) + \frac{\lambda}{\kappa} \int \frac{d^d\mathbf{q}}{(2\pi)^d} i\mathbf{q} V_{-\mathbf{q}} \langle \varphi_q e^{i\mathbf{q} \cdot \mathbf{Z}} \rangle_{\text{ss}}, \quad (49) $$

$$ \langle \varphi_q \rangle_{\text{ss}} = \frac{\lambda D_\alpha q^\alpha V_q \langle e^{-i\mathbf{q} \cdot \mathbf{Z}} \rangle_{\text{ss}}}{\alpha_q - i\mathbf{q} \cdot \mathbf{v}}, \quad (50) $$

where $\langle \cdots \rangle_{\text{ss}}$ indicates the expectation value taken in the stationary state. The fact that the system is driven out of equilibrium and that the probability density function in the stationary state is not known a priori (being generically different from the Boltzmann distribution), in conjunction with the coupling between the field $\varphi_q$ and the particle position $\mathbf{Z}$, makes it very difficult to evaluate $\langle e^{-i\mathbf{q} \cdot \mathbf{Z}} \rangle_{\text{ss}}$ and $\langle \varphi_q e^{i\mathbf{q} \cdot \mathbf{Z}} \rangle_{\text{ss}}$ on the right-hand side of equations (49) and (50). Previous studies [18, 29, 32, 52] considered a perturbative expansion of these expectation values in powers of $\lambda$ in order to derive predictions in the weak-coupling limit at the lowest non-vanishing order. Qualitatively, assuming $\lambda V(\mathbf{x}) > 0$, the resulting steady-state field profile $\langle \varphi(\mathbf{x}) \rangle_{\text{ss}}$, which is the inverse Fourier transform of $\langle \varphi_q \rangle_{\text{ss}}$ in equation (50), is enhanced in the vicinity of the particle and stretched in the direction opposite to the dragging velocity (see figure 1). We will refer to this field configuration as the shadow. Its maximum turns out to be located behind the particle, and the presence of $\langle \varphi(\mathbf{x}) \rangle_{\text{ss}}$ induces an additional friction on the tracer, causing a shift of its steady-state position, as reported in equation (25).

## 4.1. The noiseless limit

In order to investigate the qualitative features of the steady-state field configuration $\langle \varphi(\mathbf{x}) \rangle_{\text{ss}}$, it is convenient to neglect thermal noise by setting $T=0$. In this case, the dynamics in equations (12) and (14) becomes deterministic. Correspondingly,
we denote $\langle \mathbf{Z} \rangle_{ss} = \mathbf{Z}^{(ss)}$ and $\langle \varphi_q \rangle_{ss} = \varphi_q^{(ss)}$, which satisfy $\langle \varphi_q e^{i\mathbf{q}\cdot\mathbf{Z}} \rangle_{ss} = \varphi_q^{(ss)} e^{i\mathbf{q}\cdot\mathbf{Z}^{(ss)}}$ and $\langle e^{-i\mathbf{q}\cdot\mathbf{Z}} \rangle_{ss} = e^{-i\mathbf{q}\cdot\mathbf{Z}^{(ss)}}$. Accordingly, equations (49) and (50) become

$$ \mathbf{Z}^{(ss)} = -\mathbf{v}/(\kappa\nu) - \frac{\lambda^2}{\kappa} \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} D_\alpha q^\alpha |V_q|^2 \frac{\mathbf{q} \cdot \mathbf{v}}{\alpha_q^2 + (\mathbf{q} \cdot \mathbf{v})^2}, \quad (51) $$

$$ \varphi_q^{(ss)} = \frac{\lambda D_\alpha q^\alpha V_q e^{-i\mathbf{q}\cdot\mathbf{Z}^{(ss)}}}{\alpha_q - i\mathbf{q}\cdot\mathbf{v}}. \quad (52) $$

In passing, we note that equation (26)—which describes the shift of the steady-state particle position obtained at $T=0$ in a perturbative expansion in $\lambda$—coincides with the $\lambda$-dependent term in equation (51), in which no assumption regarding the magnitude of $\lambda$ was made, yielding an exact expression. Accordingly, higher-order terms in the perturbative expansion are generated by thermal effects.

We note that, since the coordinate system is chosen such that the minimum of the optical trap is at the origin at time $t=0$, the steady-state configuration of the field provides the initial condition for the field dynamics for $t > 0$, according to

$$ \phi(\mathbf{x}, t=0) = \varphi^{(ss)}(\mathbf{x}). \quad (53) $$

Moreover, as we discuss in equation (A.25), the part $\lambda \tilde{\phi}_q^{(1)}(t)$ of the field contributing to the force (8) within the perturbative expansion in equations (21) and (22) is given by

$$ \lambda \tilde{\phi}_q^{(1)}(t) = e^{-\alpha_q t} \varphi_q^{(ss)}, \quad (54) $$

for $t>0$ and in the noiseless limit $T=0$. We conclude that it is the slow relaxation of $\lambda \tilde{\phi}_q^{(1)}(t)$ that gives rise to a non-zero recoil: indeed, if the field relaxed instantaneously (corresponding to the formal limit $\alpha_q \to \infty$ taken in equation (54)), then no systematic drift of the particle would be observed after its release from the trap. We note also that plugging equation (52) into the expression for the field-induced force in equation (8) renders a time-independent force, from which the short-time expression for the recoil in equation (30) with $T=0$ can be easily retrieved.

## 4.2. The steady-state field configuration

At distances from the particle much larger than the particle size $R$, the field profile in the driven steady state can be described analytically by approximating the interaction kernel in equation (4) by a delta distribution which corresponds to the limit of a point-like particle, i.e. $R \to 0$. In particular, we focus on the behavior of the field configuration in the steady state as a function of the coordinate $y = (\mathbf{x} - \mathbf{Z}^{(ss)}) \cdot \hat{\mathbf{e}}_1$ along the dragging axis $\hat{\mathbf{e}}_1$, with $x_i = 0$ for $i \ge 2$ (where, as previously done, we assumed $\mathbf{v} = \hat{\mathbf{e}}_1 v$). For
model A, one obtains the field profile in the case of $V(\mathbf{x}) = \delta^{(d)}(\mathbf{x})$ in the form (see appendix C.1 for details)

$$ \varphi^{(\text{ss})}(y) = \begin{cases} \frac{\lambda D_0}{v} \frac{1}{\sqrt{1+(2D_0/v\xi)^2}} \exp\left[-\frac{v}{2D_0}\left(|y|\sqrt{1+(2D_0/v\xi)^2} + y\right)\right], & \text{for } d=1, \\ \frac{\lambda}{2\pi} K_0\left(\frac{|y|v}{2D_0}\sqrt{1+(2D_0/v\xi)^2}\right) \exp\left[-\frac{v}{2D_0}y\right], & \text{for } d=2, \\ \frac{\lambda}{4\pi|y|} \exp\left[-\frac{v}{2D_0}\left(|y|\sqrt{1+(2D_0/v\xi)^2} + y\right)\right], & \text{for } d=3, \end{cases} \quad (55) $$

where $K_n(x)$ is the modified Bessel function of the second kind of order $n$. Equation (55) allows us to define two length scales $l_+$ and $l_-$, describing the spatial decay of the tails of the shadow away from the particle along the direction of the driving or opposite to it, respectively, such that

$$ \varphi^{(\text{ss})}(y \to \pm\infty) \approx \exp(-|y|/l_{\pm}). \quad (56) $$

The inverse of these lengths is given by

$$ (l_{\pm})^{-1} = \frac{v}{2D_0} \left( \sqrt{1 + (2D_0/v\xi)^2} \pm 1 \right), \quad (57) $$

which coincide with the expressions reported in the supplemental material of [32]. Note that, generically, $l_- > l_+$, i.e. the field profile in front of the particle decreases more rapidly than that on its back. Moreover, upon approaching criticality, one has

$$ l_+ \xrightarrow{\xi \to \infty} D_0/v \quad \text{and} \quad l_- \xrightarrow{\xi \to \infty} +\infty. \quad (58) $$

For a particle with finite size, the maximum of the field profile is located to the left of the steady-state position of the particle, as shown in figure 6.

We note that, according to equation (57), $l_+$ is always finite, and thus the shadow in front of the particle decays exponentially for any value of $\xi$. In contrast, at the critical point, the lengthscale $l_-$ describing the decay of the field profile behind the particle diverges, and the profile exhibits an algebraic decay according to<fnref n="10" />

$$ \varphi^{(\text{ss})}(y \to -\infty) \approx \frac{\lambda}{(4\pi|y|)^{(d-1)/2}} \left(\frac{v}{D_0}\right)^{(d-3)/2}, \quad (59) $$

as shown in appendix C. Interestingly, in $d=1$, the field profile approaches a finite value as $y \to -\infty$. This behavior can be traced back to the fact that a scalar Gaussian field theory is peculiar in $d \le 2$ at the critical point, in that the long-distance behavior of its correlation function

$$ \langle \phi(\mathbf{x}) \phi(0) \rangle \sim |\mathbf{x}|^{2-d} \quad (60) $$

<footnote n="10">¹⁰ Note that equation (59) agrees with equation (55), since $e^x K_0(x) \sim \sqrt{\pi/(2x)}$ for large $x$.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Steady-state configuration of the field in the co-moving reference frame for model A dynamics, as a function of the coordinate $y$ along the driving direction, and in the noiseless limit ($T=0$). Panel (a) corresponds to $d=1$ and various values of $r$, while panel (b) corresponds to $d=3$ at criticality ($r=0$). In both panels $R=0.5$, while $\lambda$, $D_0$, and $v$ are set to one. The dashed lines in both panels correspond to the expressions in the case of a point-like particle, i.e. for $R \to 0$, reported in equation (55). Note that for $d=1$ at the critical point, the field configuration attains a non-vanishing value for $y \to -\infty$. For $d=3$, instead, the field profile decays as $\sim 1/|y|$ for $y \to -\infty$, in agreement with equation (59).</figcaption></figure>

is characterized by a linear growth of correlations in $d=1$ [55]. Similarly, we attribute the divergence of the critical recoil range reported in equation (43) to this fact.

In passing, we note that in the limit in which the (non-critical) field equilibrates quickly as the particle moves (e.g. $D_0 \to \infty$), both length scales defined in equation (57) converge to the correlation length $\xi$:

$$l_{\pm} \to \xi \quad \text{for} \quad D_0 \to \infty. \qquad (61)$$

Accordingly, in this adiabatic limit, the shape of the field profile becomes independent of the dragging velocity and, in the co-moving frame, it is the same as the equilibrium one at rest. Since the mobility of the field enters equation (57) only via the combination $v/D_0$, the adiabatic limit is actually equivalent to the limit $v \to 0$ with fixed $\xi$ and $D_0$.

In the case of model B, determining the features of the field profile is technically more demanding, and thus we postpone the discussion of the non-critical case to appendix C. Here, we report the critical behavior and the qualitative features for finite correlation lengths. At the critical point, for $d=1$ and for the delta-like interaction potential, the stationary field profile $\varphi^{(ss)}(y)$ is given by

$$\varphi^{(ss)}(y) = \begin{cases} \frac{\lambda D_2^{1/3}}{3v^{1/3}} e^{-|y|(v/D_2)^{1/3}}, & \text{for } y < 0, \\ \frac{2\lambda D_2^{1/3}}{3v^{1/3}} e^{-\frac{y}{2}(v/D_2)^{1/3}} \sin\left(\frac{\pi}{6} - \frac{\sqrt{3}}{2}y(v/D_2)^{1/3}\right), & \text{for } y > 0. \end{cases} \qquad (62)$$

We notice that, since the dynamics locally conserves the field, the integral over the whole space of the steady-state field configuration is constant; in particular, here it vanishes, because the dynamics conserves the $\phi_{q=0}$ mode, which in the calculation we
had set to zero at $t_0 \to -\infty$. In spatial dimension $d=3$ and in the limit of a point-like particle, the field profile in the direction of the driving force is given by

$$ \varphi^{(ss)}(y) = \begin{cases} \frac{\lambda}{8\pi|y|} \left(1 + e^{-|y|(v/D_2)^{1/3}}\right), & \text{for } y < 0, \\ \frac{\lambda}{4\pi y} e^{-\frac{y}{2}(v/D_2)^{1/3}} \cos\left(\frac{\sqrt{3}}{2}y(v/D_2)^{1/3}\right), & \text{for } y > 0. \end{cases} \quad (63) $$

Accordingly, at criticality, a qualitative difference between model A and model B is observed in the behavior of $\varphi^{(ss)}(y)$ for $|y| \gg R$. In fact, for model A, the order parameter decays monotonically upon increasing the distance from the particle, as in equation (55). For model B, instead, the field profile in the direction opposite to the drive still decreases monotonically, while in the direction of the drive it displays damped oscillations, as in equation (63) above. Note that, according to equations (58) and (62), when the limits $R \to 0$ and $\xi \to \infty$ are taken simultaneously—i.e. when one considers a point-like particle in a critical medium—the field profile varies in space on the lengthscale defined in equation (38), which describes the relative significance of the effects due to the drive and the relaxation of the medium. In particular, this ratio is related to the Weissenberg number $Wi$ introduced in equation (37) via $2Wi = (R/l_{v,\alpha})^{\alpha+1}$.

In the case of a non-critical medium governed by model B dynamics, the steady-state field configuration tuns out to display damped oscillations away from the particle in the direction of the drive if the velocity $v$ is sufficiently large, i.e. for

$$ \frac{l_{v,2}}{\xi} < \frac{3^{1/2}}{2^{1/3}} \equiv a^* \approx 1.37, \quad (64) $$

irrespective of the spatial dimensionality of the system. In particular, the behavior of the field for $y > 0$ in the case of model B, in $d=1$, and in the limit of a point-like particle, as a function of

$$ a \equiv l_{v,2}/\xi, \quad (65) $$

can be written as

$$ \varphi^{(ss)}(y) \propto \begin{cases} \exp[-\gamma(a)y/l_{v,2}] \cos[k(a)y/l_{v,2} - \psi(a)], & \text{for } a < a^*, \\ \exp[-\gamma_1(a)y/l_{v,2}] - A(a) \exp[-\gamma_2(a)y/l_{v,2}], & \text{for } a > a^*, \end{cases} \quad (66) $$

where $1 > A(a) > 0$, $\gamma(a) > 0$, $\gamma_1(a) > \gamma_2(a) > 0$, and $k(a)$ is a dimensionless wavevector to be discussed below. Accordingly, for $l_{v,2}/\xi > a^*$, the field profile is described by a superposition of two exponential functions with amplitudes of opposite signs. Notably, the function with the negative amplitude decays more slowly than the other upon increasing the distance from the particle. Explicit expressions for the field profile are presented in appendix C.2. The dimensionless wavevector $k$ introduced in equation (66) is plotted as a function of $a$ (see equation (65)) in figure 7. In the case of $d > 1$, the transition from the damped-oscillations regime to that described by a superposition of two exponential functions also occurs at $a = a^*$ (see equation (64)).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The dimensionless wavevector $k$ of the damped oscillations of the field profile in the direction of the drag for model B and a delta-like interaction kernel, defined in equation (66) as a function of $a \equiv l_{v,2}/\xi$ (see equation (65)). Spatial oscillations are observed for $l_{v,2}/\xi < a^*$ (see equation (64)), and the wavevector $k$ attains its maximal value $k_{max} = \sqrt{3}/2$ in the limit $l_{v,2}/\xi \to 0$. When $l_{v,2}/\xi$ approaches the value $a^*$ from below, the wavevector $k$ vanishes according to $k \sim \sqrt{a^* - (l_{v,2}/\xi)}$.</figcaption></figure>

Representative plots of the field profiles for model A and B are reported in figures 6 and 8, respectively, which assume the Gaussian interaction kernel reported in equation (5). In particular, in panels (a) of both plots, the one-dimensional case is reported with various values of the parameter $r$, which controls the correlation length of the medium. The dashed lines in the plots correspond to the expressions in the $R \to 0$ limit, given in equations (55), (C.8) and (C.7), respectively. The plots of the steady-state field profile along the driving direction in spatial dimension $d=3$ at criticality are reported in panels (b) of figures 6 and 8, together with the corresponding asymptotic behaviors given in equations (55) and (63).

Finally, we note that in $d=3$, irrespective of the choice of the dynamics of the field and of the interaction kernel, taking the limit $\xi \to \infty$ and subsequently $v \to 0$ yields, for the steady-state field profile along the direction of the driving,

$$ \varphi^{(ss)}(y \to \pm\infty) \approx \frac{\lambda}{4\pi|y|}, \qquad (67) $$

which, as expected, corresponds to the equilibrium configuration of the field with a point-like source. This indicates that in spatial dimension $d=3$, the limits $\xi \to \infty$ and $v \to 0$ commute, which is not the case for $d=1$ or $2$. Again, the emergence of this issue can be traced back to the previously-discussed peculiarity of the scalar Gaussian field theory at the critical point, see equation (60).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Steady-state configuration $\varphi^{(ss)}(y)$ of the field in the co-moving reference frame for model B, as a function of the coordinate $y$ along the driving direction, and in the noiseless limit ($T=0$). Panel (a) corresponds to $d=1$ with $R=0.5$ and various values of $r$, while panel (b) corresponds to $d=3$ with $R=0.2$ and at criticality. In both panels $\lambda, D_0$, and $v$ are set to one. The dashed lines in both panels correspond to the expressions in the case of a point-like particle, i.e. for $R \to 0$, reported in equations (C.7), (C.8) and (63). Note that for $d=1$, the field configuration is characterized by overdamped oscillations for $r=0$ and $0.5$, and by a superposition of two exponential functions of opposite signs for $r=2$. These different qualitative behaviors are found in agreement with the condition for the presence of oscillations given in equation (64).</figcaption></figure>

## 5. Conclusions

In this work, we described the dynamics of an overdamped Brownian particle coupled to a fluctuating scalar Gaussian field, characterized by a correlation length $\xi$. The field is assumed to have relaxational dynamics, with the possible presence of a local conservation law. The particle and the field are linearly coupled to each other, with the strength of their interaction determined by the coupling constant $\lambda$, and both are in contact with a thermal bath at temperature $T$. At times $t<0$, the particle is driven across the medium by a harmonic trap moving with constant velocity $v$, and it is released from the trap at $t=0$. In describing the effective dynamics of the tracer particle, which is non-linear and non-Markovian, we resorted to a perturbative expansion in the coupling $\lambda$. We demonstrated that, after the trap is turned off at time $t=0$, the tracer exhibits recoil, i.e. it moves in the direction opposite to that of the dragging velocity.

After deriving an expression for the average recoil in equation (28), we focused on the noiseless limit corresponding to $T=0$, for which, in equation (32), we introduced a scaling function $F^{A/B}(R/\xi, t/\tau_R, v\tau_R/R)$—where $\tau_R$ is the typical relaxation time of a critical field over a distance of the order of the effective particle size $R$. For small velocities $v$, the amplitude of the final recoil is $\propto v$ as long as $t/\tau_R$ and $\xi/R$ remain finite. The non-trivial small-$v$ behavior arising for $\xi \to \infty$ and $t \to \infty$ is summarized in table 1. We showed that, generically, the recoil eventually vanishes upon increasing $v$, due to an effective decoupling of the particle from the field. In the critical limit $\xi \to \infty$, the eventual recoil has a finite amplitude for $d>2$, and it reaches its maximal
value following an algebraic law $\sim t^{(2-d)/z}$, where $z$ is the dynamical critical exponent of the field, as plotted for $d=3$ in figures 3(a) and 4(a) for dissipative (model A) and conservative (model B) dynamics, respectively. For $d \le 2$, as the system approaches its critical point, we described the divergence of the recoil range $\propto \xi^{2-d}$. Conversely, in the non-critical case, the recoil has a finite range for any spatial dimensionality $d$. For non-conserved dynamics, the final value of the recoil is approached exponentially in time, whereas for conserved dynamics, the approach is algebraic $\propto t^{-d/2}$. A comprehensive summary of the long-time behavior of the scaling function $F^{A/B}$ is reported in table 2.

As it drives the initial evolution of the recoil, we investigated the spatial configuration of the field which develops around the particle in the steady state generated by the driving at time $t<0$, considering it in the reference frame co-moving with the trap. In the case of non-conserved dynamics and at $T=0$, two lengthscales $l_\pm$ naturally emerge, which describe the rate of the exponential decay towards zero of the field upon moving away from the particle in the direction ($l_+$) or opposite ($l_-$) to the driving, i.e. in front or behind the dragged particle. In the critical case $\xi \to \infty$, $l_+$ remains finite whereas $l_-$ diverges, and the decay of the field profile is characterized by an algebraic law $\propto 1/|y|^{(d-1)/2}$, where $y$ is the distance along the axis of the drive from the steady-state position of the tracer in the co-moving frame (see equation (59)). We noted that, in the adiabatic limit in which the field equilibrates instantaneously around the current position of the particle, both $l_+$ and $l_-$ reduce to the correlation length $\xi$. In the case of conserved dynamics, depending on the values of the various parameters, it is possible to observe damped spatial oscillations of the average field profile in front of the dragged particle, whose wavevector is plotted in figure 7.

This work may motivate an experimental study of the recoil of colloidal particles dragged through a near-critical binary mixture, in a manner similar to the investigations reported for example in [21] for viscoelastic media. Colloidal particles in critical media are, by now, easily accessible experimentally [35–39], and the values of the correlation length that can be reached in these experiments are of the order of hundreds of nanometers, which is comparable with the typical size of the colloidal particles. Also, in experiments in viscoelastic media [8], the dragging velocities considered spanned the range between $0.4 \ \mu\text{m s}^{-1}$ and $40 \ \mu\text{m s}^{-1}$. In these experimental conditions, it should be possible to test, at least qualitatively, the behavior of the recoil predicted in this work. However, a more quantitative description would require accounting for the effect of hydrodynamics and of the self-interaction $\propto \phi^4$ of the fluctuating order-parameter field. These issues are left for future investigations.

Finally, the framework adopted in this work could be adapted to describe other phenomena exhibited by a tracer particle coupled to a slowly relaxing structured medium—such as the *Magnus effect* recently observed with rotating particles in viscoelastic media [21]. This possibility calls for an extension of the model discussed here in order to describe a particle and a medium that are subject not only to a translational, but also to a rotational shear. Moreover, the discussion may be extended to systems in which the fluctuation-dissipation relation is not satisfied, e.g. by considering an active field, modeling a non-equilibrium bath of active particles (as investigated, for example, in [56–60]), or conversely by considering an active particle evolving in contact with a passive field.
Acknowledgments

We thank Clemens Bechinger, Vincent Démery, Aljaz Godec, Matthias Krüger, Christian Maes, and Gennaro Tucci for stimulating discussions. AG acknowledges support from the MIUR PRIN project ‘Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST)’ n. 201798CZL.

Appendix A. Perturbative calculation of the correction to the particle’s position

In this appendix we present the derivation of the correction to the particle position X due to the interaction (8) with the field. In particular, in appendix A.1 this is done while the particle is driven by the trap, resulting in the time-independent shift reported in equation (25), and in appendix A.2, after the particle has been released from the confinement, yielding the recoil reported in equation (28).

A.1. Driven particle

For $t \le 0$, the particle is dragged by the moving trap. As we are interested in the steady state attained at long times in the co-moving frame of the particle, we formally take the limit $t_0 \to -\infty$. When the particle and the field are decoupled (i.e. with $\lambda = 0$), according to equations (12) and (14), their position and configuration are given by

$$ \mathbf{X}^{(0)}(t) = \mathbf{v}\left(t - \frac{1}{\nu\kappa}\right) + \mathbf{O}(t), \qquad (A.1) $$

$$ \phi_{\mathbf{q}}^{(0)}(t) = \int_{-\infty}^{t} dt' e^{-\alpha_{q}(t-t')} \eta_{\mathbf{q}}(t'), \qquad (A.2) $$

where $\mathbf{O}(t) = \int_{-\infty}^{t} dt' e^{-\nu\kappa(t-t')} \xi(t')$ is the Ornstein–Uhlenbeck process. This Gaussian process is characterized by a vanishing mean and two-point correlator [51]

$$ \langle O_i(t_1) O_j(t_2) \rangle = \frac{T}{\kappa} e^{-\nu\kappa|t_2-t_1|} \delta_{i,j}. \qquad (A.3) $$

Note that, since the initial conditions are imposed at $t_0 \to -\infty$, they are irrelevant in determining the steady-state configuration of the field in equation (A.2). Moreover, in equation (A.1), we chose the coordinate system such that the minimum of the trap is at the origin at $t=0$.

Using the expressions for the decoupled case (see equations (A.1) and (A.2)), we solve equations (12) and (14) term-by-term within the perturbative expansion in $\lambda$ introduced in equations (21) and (22). At the first order in $\lambda$, we obtain

$$ \mathbf{X}^{(1)}(t) = i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} \int_{-\infty}^{t} dt' e^{-\nu\kappa(t-t') + i\mathbf{q} \cdot [\mathbf{v}(t'-\frac{1}{\nu\kappa}) + \mathbf{O}(t')] } \phi_{\mathbf{q}}^{(0)}(t'), \qquad (A.4) $$

$$ \phi_{\mathbf{q}}^{(1)}(t) = D_{\alpha} q^{\alpha} V_{\mathbf{q}} \int_{-\infty}^{t} dt' e^{-\alpha_{q}(t-t')} - i \mathbf{q} \cdot [\mathbf{v}(t'-\frac{1}{\nu\kappa}) + \mathbf{O}(t')]. \qquad (A.5) $$
At the second order in $\lambda$, equation (14) for the particle position renders the expression reported in equation (24). We evaluate the average of each term therein separately, and write

$$i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} \left\langle e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} \phi_{\mathbf{q}}^{(1)}(t) \right\rangle \\ = i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} |V_{\mathbf{q}}|^2 D_{\alpha} q^{\alpha} \mathbf{q} \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t') + i\mathbf{q} \cdot \mathbf{v}(t-t')} \left\langle e^{i\mathbf{q} \cdot [\mathbf{O}(t)-\mathbf{O}(t')]} \right\rangle. \quad (A.6)$$

We now recall the property of the Ornstein–Uhlenbeck process (see e.g. appendix C.1 of [18]),

$$\left\langle e^{i\mathbf{q} \cdot [\mathbf{O}(t)-\mathbf{O}(t')]} \right\rangle = e^{-q^2(T/\kappa)[1-\exp(-\nu\kappa|t-t'|)]}, \quad (A.7)$$

and arrive at

$$i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} \left\langle \mathbf{q} e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} V_{-\mathbf{q}} \phi_{\mathbf{q}}^{(1)}(t) \right\rangle \\ = i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} |V_{\mathbf{q}}|^2 \frac{\alpha_q}{q^2 + r} \mathbf{q} \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t') + i\mathbf{q} \cdot \mathbf{v}(t-t') - q^2(T/\kappa) [1 - e^{-\nu\kappa(t-t')}]} . \quad (A.8)$$

For the second term in equation (24), we first average over the noise $\eta$ and find

$$i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} \mathbf{q} e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} V_{-\mathbf{q}} i\mathbf{q} \cdot \left\langle \mathbf{X}^{(1)}(t) \phi_{\mathbf{q}}^{(0)}(t) \right\rangle_{\eta} = -i\nu^2 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \int \frac{d^d \mathbf{q}'}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} V_{-\mathbf{q}'} (\mathbf{q} \cdot \mathbf{q}') \\ \times \int_{-\infty}^{t} dt' e^{-\nu\kappa(t-t') + i\mathbf{q} \cdot [\mathbf{v}(t-\frac{1}{\nu\kappa}) + \mathbf{O}(t)] + i\mathbf{q}' \cdot [\mathbf{v}(t'-\frac{1}{\nu\kappa}) + \mathbf{O}(t')]} \left\langle \phi_{\mathbf{q}}^{(0)}(t) \phi_{\mathbf{q}'}^{(0)}(t') \right\rangle_{\eta} \\ = i\nu^2 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \mathbf{q} |V_{\mathbf{q}}|^2 q^2 \frac{T}{q^2 + r} \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t') - \nu\kappa(t-t') + i\mathbf{q} \cdot [\mathbf{v}(t-t') + \mathbf{O}(t) - \mathbf{O}(t')]} , \quad (A.9)$$

where we used the correlator of the Fourier transform of the field's noise in equation (13) to derive, in the long-time limit [51],

$$\left\langle \phi_{\mathbf{q}}^{(0)}(t) \phi_{\mathbf{q}'}^{(0)}(t') \right\rangle_{\eta} = (2\pi)^d \delta^{(d)}(\mathbf{q} + \mathbf{q}') \frac{T}{q^2 + r} e^{-\alpha_q |t-t'|}. \quad (A.10)$$

Using again the property (A.7) of the Ornstein–Uhlenbeck process, we arrive at

$$i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} \left\langle e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} i\mathbf{q} \cdot \mathbf{X}^{(1)}(t) \phi_{\mathbf{q}}^{(0)}(t) \right\rangle \\ = i\nu^2 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \mathbf{q} |V_{\mathbf{q}}|^2 q^2 \frac{T}{q^2 + r} \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t') - \nu\kappa(t-t') + i\mathbf{q} \cdot \mathbf{v}(t-t') - q^2(T/\kappa) [1 - e^{-\nu\kappa(t-t')}]} . \quad (A.11)$$
Combining this with equation (A.8) allows us to write

$$i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} \mathbf{q} \left\langle e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} V_{-\mathbf{q}} \left[ \phi_{\mathbf{q}}^{(1)}(t) + i\mathbf{q} \cdot \mathbf{X}^{(1)}(t) \phi_{\mathbf{q}}^{(0)}(t) \right] \right\rangle \\ = i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} |V_{\mathbf{q}}|^2 \frac{1}{q^2 + r} \int_0^\infty ds e^{-\alpha_q s - q^2(T/\kappa)(1-e^{-\nu\kappa s})} [\alpha_q + q^2 \nu T e^{-\nu\kappa s}] e^{i\mathbf{q} \cdot \mathbf{v}s}. \quad (A.12)$$

Noting that

$$\frac{d}{ds} e^{-\alpha_q s - q^2(T/\kappa)(1 - e^{-\nu\kappa s})} = -[\alpha_q + q^2 \nu T e^{-\nu\kappa s}] e^{-\alpha_q s - q^2(T/\kappa)(1 - e^{-\nu\kappa s})}, \quad (A.13)$$

the integral over $s$ in equation (A.12) can be computed by parts. The boundary terms vanish, due to the fact that the $\mathbf{q}$-integrand evaluated at $s=0$ is odd under $\mathbf{q} \mapsto -\mathbf{q}$. This yields

$$i\nu \int \frac{d^d \mathbf{q}}{(2\pi)^d} \mathbf{q} \left\langle e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} V_{-\mathbf{q}} \left[ \phi_{\mathbf{q}}^{(1)}(t) + i\mathbf{q} \cdot \mathbf{X}^{(1)}(t) \phi_{\mathbf{q}}^{(0)}(t) \right] \right\rangle \\ = -\nu v \hat{\mathbf{e}}_1 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + r} \int_0^\infty ds \, e^{-(\alpha_q - i q_1 v)s - q^2(T/\kappa)(1 - e^{-\nu\kappa s})}, \quad (A.14)$$

where we assumed $\mathbf{v} = v\hat{\mathbf{e}}_1$. Using equation (24), this leads to

$$\left\langle \dot{\mathbf{X}}^{(2)}(t) \right\rangle + \nu\kappa \left\langle \mathbf{X}^{(2)}(t) \right\rangle = -\nu v \hat{\mathbf{e}}_1 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + r} \int_0^\infty ds \, e^{-(\alpha_q - i q_1 v)s - q^2(T/\kappa)(1 - e^{-\nu\kappa s})}. \quad (A.15)$$

This equation can be readily integrated by using the corresponding Green's function, leading to the expression for the time-independent shift reported in equation (25).

Sign of the shift—In order to prove that the shift reported in equation (25) is indeed in the direction opposite to the drag, it is sufficient to show that

$$\int_0^\infty ds \cos(q_1 vs) e^{-\alpha_q s - (T/\kappa)q_1^2(1 - e^{-\nu\kappa s})} > 0, \quad (A.16)$$

because the imaginary part of the $\mathbf{q}$-integral in equation (25) vanishes since the integrand is odd under $\mathbf{q} \mapsto -\mathbf{q}$. Let us denote

$$f(s) = e^{-as - b(1 - e^{-cs})}, \quad (A.17)$$

with positive $a$, $b$, and $c$. By expanding the exponential in terms of its exponential argument, it is easy to prove that the function $f(s)$ can be expressed in the form

$$f(s) = \int_{0}^{\infty} dt \, e^{-st} e^{-b} \sum_{n=0}^{\infty} \frac{b^n}{n!} \delta(t - a - nc), \quad (A.18)$$
from which it follows that

$$ \int_0^\infty ds \cos(q_1 vs) f(s) = e^{-b} \sum_{n=0}^\infty \frac{b^n}{n!} \frac{a+nc}{(a+nc)^2 + (q_1 v)^2} > 0. \quad (A.19) $$

Choosing $a = \alpha_q$, $b = T/\kappa$, and $c = \kappa\nu$ proves the validity of equation (A.16). Accordingly, the shift of the average particle position in the driven phase of the motion reported in equation (25) is in the direction opposite to the drag.

## A.2. Recoil after the particle is released

In deriving the expression of the particle position $X$ after the particle is released from the trap ($t > 0$), the solutions obtained for $t \le 0$ serve as initial conditions. In particular, the field configuration at $t = 0$ is the one given by equation (A.5). When the particle and the field are decoupled (i.e. for $\lambda = 0$), one obtains from equations (12) and (14),

$$ \mathbf{X}^{(0)}(t) = \mathbf{X}^{(0)}(0) + \mathbf{W}(t), \quad (A.20) $$

$$ \phi_{\mathbf{q}}^{(0)}(t) = \int_{-\infty}^{t} dt' e^{-\alpha_{\mathbf{q}}(t-t')} \eta_{\mathbf{q}}(t'), \quad (A.21) $$

where $\mathbf{W}(t) = \int_0^t dt' \xi(t')$ is the Wiener process. This Gaussian process has a vanishing mean and a variance given by

$$ \langle W_i(t_1)W_j(t_2) \rangle = d\nu T \min(t_1, t_2) \delta_{i,j}. \quad (A.22) $$

Again, using the results obtained in the decoupled case in equations (A.20) and (A.21), we solve equations (12) and (14) term-by-term according to the $\lambda$-expansion in equations (21) and (22). At the first order in $\lambda$, we obtain

$$ \mathbf{X}^{(1)}(t) = \mathbf{X}^{(1)}(0) + i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} e^{i\mathbf{q}\cdot\mathbf{X}^{(0)}(0)} \int_0^t dt' e^{i\mathbf{q}\cdot\mathbf{W}(t')} \phi_{\mathbf{q}}^{(0)}(t'), \quad (A.23) $$

$$ \phi_{\mathbf{q}}^{(1)}(t) = e^{-\alpha_q t} \phi_{\mathbf{q}}^{(1)}(0) + D_{\alpha} q^{\alpha} V_{\mathbf{q}} e^{-i\mathbf{q}\cdot\mathbf{X}^{(0)}(0)} \int_{0}^{t} dt' e^{-\alpha_q(t-t')-i\mathbf{q}\cdot\mathbf{W}(t')}, \quad (A.24) $$

We note that the first term in equation (A.24) is the contribution to the field from the relaxation over a time $t$ of its steady-state configuration, characterized by the average reported in equation (50). To proceed, we recall equation (24) and write

$$ i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} \langle e^{i\mathbf{q}\cdot\mathbf{X}^{(0)}(t)} \phi_{\mathbf{q}}^{(1)}(t) \rangle = i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} e^{-\alpha_q t} \langle e^{i\mathbf{q}\cdot\mathbf{W}(t)} \rangle \langle e^{i\mathbf{q}\cdot\mathbf{X}^{(0)}(0)} \phi_{\mathbf{q}}^{(1)}(0) \rangle \\ + i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} |V_{\mathbf{q}}|^2 D_{\alpha} q^{\alpha} \int_{0}^{\infty} dt' e^{-\alpha_q(t-t')} \langle e^{-i\mathbf{q}\cdot[\mathbf{W}(t)-\mathbf{W}(t')]} \rangle, \quad (A.25) $$
where we used the fact that the Wiener process performed by the particle after the trap is turned off is statistically independent from the dynamics at previous times. We now note that

$$ \langle e^{i\mathbf{q} \cdot [\mathbf{W}(t) - \mathbf{W}(t')]} \rangle = e^{-q^2 \nu T |t-t'|}, \qquad (A.26) $$

from which

$$ \langle e^{i\mathbf{q} \cdot \mathbf{W}(t)} \rangle = e^{-q^2 \nu T t}, \qquad (A.27) $$

which is readily obtained by taking the limit $\kappa \to 0$ in equation (A.7). Note that the integrand in the second term of equation (A.25) is odd under $\mathbf{q} \mapsto -\mathbf{q}$, and therefore it vanishes. Accordingly, only the first term (stemming from the relaxation of the steady-state configuration) contributes to the force acting on the particle, as we mentioned above equation (54). We may now use the result obtained for $t \le 0$ in equation (A.5) and write

$$ \begin{aligned} i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \langle \mathbf{q} e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} V_{-\mathbf{q}} \phi_{\mathbf{q}}^{(1)}(t) \rangle \\ = i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} |V_{-\mathbf{q}}|^2 \frac{\alpha_q}{q^2+r} e^{-(\alpha_q+q^2\nu T)t} \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t') + i\mathbf{q} \cdot \mathbf{v}(t-t') - q^2(T/\kappa) [1-e^{-\nu\kappa(t-t')}]}, \end{aligned} \qquad (A.28) $$

where we used $D_\alpha q^\alpha = \alpha_q/(q^2+r)$ and equation (A.7). Similarly, the second term in equation (24) can be evaluated as

$$ \begin{aligned} i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} \langle e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} i\mathbf{q} \cdot \mathbf{X}^{(1)}(t) \phi_{\mathbf{q}}^{(0)}(t) \rangle \\ = i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} \langle e^{i\mathbf{q} \cdot \mathbf{W}(t)} \rangle \langle e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(0)} i\mathbf{q} \cdot \mathbf{X}^{(1)}(0) \phi_{\mathbf{q}}^{(0)}(t) \rangle \\ + i\nu^2 \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} |V_{-\mathbf{q}}|^2 q^2 \frac{T}{q^2+r} \int_0^t dt' e^{-(\alpha_q+q^2\nu T)(t-t')}. \end{aligned} \qquad (A.29) $$

Again, the second term vanishes because its integrand is odd under $\mathbf{q} \to -\mathbf{q}$. We also note that

$$ \phi_{\mathbf{q}}^{(0)}(t) = e^{-\alpha_q t} \phi_{\mathbf{q}}^{(0)}(0) + \int_0^t ds \, e^{-\alpha_q (t-s)} \eta_{\mathbf{q}}(s), \qquad (A.30) $$

where the second term is statistically independent of all other quantities evaluated at $t=0$. Hence, using equation (A.27), we can write

$$ \begin{aligned} i\nu \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} V_{-\mathbf{q}} \langle e^{i\mathbf{q} \cdot \mathbf{X}^{(0)}(t)} i\mathbf{q} \cdot \mathbf{X}^{(1)}(t) \phi_{\mathbf{q}}^{(0)}(t) \rangle = i\nu^2 \int \frac{d^d\mathbf{q}}{(2\pi)^d} \mathbf{q} |V_{-\mathbf{q}}|^2 q^2 \frac{T}{q^2+r} e^{-(\alpha_q+q^2\nu T)t} \\ \times \int_{-\infty}^{t} dt' e^{-\alpha_q(t-t') - \nu\kappa(t-t') + i\mathbf{q} \cdot \mathbf{v}(t-t') - q^2(T/\kappa) (1-e^{-\nu\kappa(t-t')})}. \end{aligned} \qquad (A.31) $$
This allows us to cast the average of equation (24) in the form

$$ \langle \dot{\mathbf{X}}^{(2)}(t) \rangle = -\nu v \hat{e}_1 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + r} e^{-(\alpha_q + q^2 \nu T)t} \int_0^\infty ds \, e^{-(\alpha_q - i q_1 v)s - q^2(T/\kappa)(1 - e^{-\nu \kappa s})}. \quad (A.32) $$

Simple integration over $t$ yields the final result for the position of the particle, reported in equation (28).

## Appendix B. Analysis of the recoil

In this appendix, we further characterize the recoil discussed in sections 3.3 and 3.4, by deriving explicit expressions and evaluating asymptotic behaviors.

### B.1. Small-v behavior of $\mathcal{F}^{A/B}$

In this appendix, we derive the small-$v$ behavior of the scaling function $\mathcal{F}^{A/B}$ defined in equations (32) and (33), characterizing the recoil in the noiseless limit ($T \to 0$). The non-trivial dependence of its scaling properties on the spatial dimensionality $d$, reported in table 1, stems from the fact that the limit $v \to 0$ may not commute with the limits $\xi \to \infty$ and $t \to \infty$.

In this and the following appendices, we drop the symbol tilde on top of variables $\tilde{\xi}, \tilde{t}$ and $\tilde{v}$ defined in equations (36) and (37), understanding that henceforth $t, \xi$, and $v$ are dimensionless. In order to obtain the asymptotic small-$v$ behavior of the scaling functions defined in equation (33), we use the Mellin transform. For a function $f(u)$ defined for $u \ge 0$, the Mellin transform $\mathcal{M}\{f\}(s)$ is defined as [61]

$$ \mathcal{M}\{f\}(s) = \int_0^\infty du \, u^{s-1} f(u), \quad (B.1) $$

where $s \in \mathbb{C}$ within a certain domain $\mathcal{D}$ of analiticity, with the inverse

$$ f(u) = \int_{c-i\infty}^{c+i\infty} \frac{ds}{2\pi i} u^{-s} \mathcal{M}\{f\}(s), \quad (B.2) $$

where the integral runs within $\mathcal{D}$. We perform the transform (from the variable $v$ to the variable $s$) of the denominator of the function on the right-hand side of equation (33), obtaining

$$ \mathcal{M}\left\{\frac{1}{\alpha_q^2 + q_1^2 v^2}\right\}(s) = \frac{\pi}{2} \frac{1}{\sin(\pi s/2)} \frac{1}{\alpha_q^2} \left(\frac{\alpha_q}{q_1}\right)^s, \quad \text{with } 0 < \text{Re}(s) < 2, \quad (B.3) $$

where $\alpha_q = q^\alpha (q^2 + 1/\xi^2)$ (see equation (16)). Note that $1/\sin(\pi s/2)$ has simple poles when $s$ equals an even integer, in particular for $s=0$ and $s=2$ which are at the boundary of the domain of analiticity of the transform above. According to equation (B.2), these latter two poles determine the behavior $\sim v^0$ for $v \ll 1$, and $\sim v^{-2}$ for large $v$ of the scaling function defined in equation (33). The Fourier transform of the Gaussian interaction kernel introduced in equation (5), given by $\hat{V}_q = e^{-q^2/2}$, regularizes the large-$q$
<table><caption>**Table B1.** Behavior for small **q** of the integrands of the Mellin transforms (from the variable *v* to *s*) of the scaling function in equations (B.5) and (B.6). The poles that emerge after integration of these functions, according to the inverse Mellin transform in equation (B.2), correspond to the small-*v* behavior of the scaling function reported in table 1.</caption><thead><tr><th>α = 0</th><th>ξ &lt; +∞</th><th>ξ → ∞</th></tr></thead><tbody><tr><td>t &lt; +∞</td><td>~ q<sup>d-s+1</sup></td><td>~ q<sup>d+s-3</sup></td></tr><tr><td>t → ∞</td><td>~ q<sup>d-s+1</sup></td><td>~ q<sup>d+s-5</sup></td></tr></tbody><thead><tr><th>α = 2</th><th>ξ &lt; +∞</th><th>ξ → ∞</th></tr></thead><tbody><tr><td>t &lt; +∞</td><td>~ q<sup>d+s-1</sup></td><td>~ q<sup>d+3s-5</sup></td></tr><tr><td>t → ∞</td><td>~ q<sup>d+s-3</sup></td><td>~ q<sup>d+3s-9</sup></td></tr></tbody></table>

behavior of the integrand in equation (33). Accordingly, any non-trivial analytic behavior can arise solely from possible infrared divergences ($q \to 0$) in the $v \to 0$ limit. This means that, for any finite $t$, in order to study the small-$q$ behavior of equation (33), one can use $1 - e^{-\alpha_q t} \approx \alpha_q t$ for critical model A, and for model B (for any value of $\xi$). The integrand of equation (33) in the case of non-critical model A turns out to exhibit the same small-$q$ behavior both for finite and infinite $t$—this follows from the fact that its modes are characterized by a finite relaxation time at low wavevectors, see equation (16). Also, the angular integral in equation (33) is inconsequential in the discussion of possible infrared divergences present for $v \to 0$. In fact, note that in

$$ \mathcal{F}^{A/B} \left( \frac{1}{\xi}, t, v \right) \underset{v \to 0}{\longrightarrow} \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + 1/\xi^2} |\hat{V}_{\mathbf{q}}|^2 \frac{1 - e^{-tq^\alpha (q^2+1/\xi^2)}}{[q^\alpha (q^2+1/\xi^2)]^2}, \quad (B.4) $$

the angular integral only involves the term $q_1^2$ at the numerator and it always yields a finite result. Hence, we focus on the radial integral in equation (33), i.e.

$$ \mathcal{M} \left\{ \mathcal{F}^{A/B} (1/\xi, t) \right\} (s) \sim \frac{t}{\sin(\pi s/2)} \int_0^\infty dq \, q^{d-1} e^{-q^2} q^\alpha \left( \frac{q_1}{\alpha_q} \right)^{2-s}, \quad t < +\infty, \quad (B.5) $$

where we used the expression for $\alpha_q$ in equation (16), and

$$ \mathcal{M} \left\{ \mathcal{F}^{A/B} (1/\xi, t) \right\} (s) \sim \frac{1}{\sin(\pi s/2)} \int_0^\infty dq \, q^{d-1} \frac{e^{-q^2}}{q^2 + 1/\xi^2} \left( \frac{q_1}{\alpha_q} \right)^{2-s}, \quad \text{for } t \to \infty. \quad (B.6) $$

We recall that, for small-$q$, one has $\alpha_q \sim q^\alpha$ for finite $\xi$, and $\alpha_q \sim q^{\alpha+2}$ for $\xi \to \infty$. This yields the small-$q$ asymptotics of the integrands reported in table B1. Note that whenever the exponent governing the small-$q$ behavior of the integrand is equal to -1, the corresponding function has a simple pole. This observation, in conjunction with the expression for the inverse Mellin transform in equation (B.2), leads to the small-$v$ behavior of the scaling function reported in table 1.
## B.2. Long-t and large-ξ behavior of $\mathcal{F}^{A/B}$

In this appendix, we derive the asymptotic behavior of the scaling function of the recoil defined in equations (32) and (33) at long times, found in equations (40)–(42), and for large correlation lengths, after taking $t \to \infty$, as reported in equation (43). In the derivation presented below, a finite value of $v$ is assumed in order to avoid the issues discussed in appendix B.1.

Here we adopt the short-hand notation introduced in appendix B.1 and we focus on the case of finite $v$. Assuming that $\mathcal{R}^A(0, v) = \mathcal{F}^A(0, t \to \infty, v)$ exists, we can write

$$ \begin{aligned} \mathcal{R}^A(0, v) - \mathcal{F}^A(0, t, v) = & \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2} |\hat{V}_{\mathbf{q}}|^2 \frac{e^{-q^2 t}}{q^4 + q_1^2 v^2} = t^{1-d/2} \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2} |\hat{V}_{\frac{q}{\sqrt{t}}}|^2 \frac{e^{-q^2}}{q^4/t + q_1^2 v^2} \\ & \sim \frac{t^{(2-d)/2}}{v^2} |\hat{V}_0|^2 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{e^{-q^2}}{q^2} = \frac{t^{(d-2)/2}}{v^2} \frac{2|\hat{V}_0|^2}{(4\pi)^{d/2} (d-2)}. \end{aligned} \quad (B.7) $$

From the scaling $\sim t^{(2-d)/2}$ in the second line of this expression, it follows that $\mathcal{F}^A(0, t \to \infty, v)$ is finite only for $d \ge 2$. Moreover, the $q$-integral in the last line is divergent in the vicinity of $q \sim 0$ for $d \le 2$. The cases of $d=1$ or $2$ are studied separately in appendix B.3, for the specific choice of the Gaussian interaction kernel introduced in equation (5). Note that the right-hand side of equation (B.7) is proportional to $1/v^2$ for any dimensionality. Moreover, the factor $2|\hat{V}_0|^2/[(4\pi)^{d/2}(d-2)]$ in front of equation (B.7) changes sign as a function of the dimensionality $d$. In particular, it is positive for $d > 2$ and negative for $d < 2$. Accordingly, the analytic continuation to $d < 2$ describes the asymptotic long-time divergence of the critical recoil which was verified by a direct calculation of the critical scaling function in the case of model A in $d=1$ dimensions, see equation (B.13).

Similarly, in the case of model B at criticality, one finds

$$ \begin{aligned} \mathcal{R}^B(0, v) - \mathcal{F}^B(0, t, v) = & \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2} |\hat{V}_{\mathbf{q}}|^2 \frac{e^{-q^4 t}}{q^8 + q_1^2 v^2} \\ = & t^{1/2-d/4} \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2} |\hat{V}_{\frac{q}{t^{1/4}}}|^2 \frac{e^{-q^4}}{q^8/t^{3/2} + q_1^2 v^2} \\ \sim & \frac{t^{(2-d)/4}}{v^2} |\hat{V}_0|^2 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{e^{-q^4}}{q^2} = \frac{t^{(2-d)/4}}{v^2} \frac{4\pi^{1/2} |\hat{V}_0|^2}{(8\pi)^{d/2} \Gamma(d/4) (d-2)}, \end{aligned} \quad (B.8) $$

where we assumed $d > 2$ in order to evaluate the integral in the last line. Similarly to the case of model A discussed above, the analytic continuation of the expression for $d < 2$ yields a correct description of the long-time divergence of the recoil in the case of critical model B at $d=1$. In order to simplify the notation, we introduce

$$ \mathcal{C}_{\alpha}(d) = \begin{cases} 2/[(4\pi)^{d/2}(2-d)], & \text{for } \alpha = 0, \\ 4\pi^{1/2}/[(8\pi)^{d/2}\Gamma(d/4)(2-d)], & \text{for } \alpha = 2, \end{cases} \quad (B.9) $$
as we did in equation (40). Both $C_0(d)$ and $C_2(d)$ have simple poles at $d=2$ with the same residue $-1/(2\pi)$. At $d=2$, instead of power laws, one observes a logarithmic growth in time of the scaling function $\mathcal{F}^{A/B}$ in equation (33) for large times (see cf equation (B.15)).

For non-critical model A, the long-$t$ behavior is obtained as

$$ \begin{aligned} \mathcal{R}^A(1/\xi, v) - \mathcal{F}^A(1/\xi, t, v) &= \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + 1/\xi^2} |\hat{V}_{\mathbf{q}}|^2 \frac{e^{-(q^2+1/\xi^2)t}}{(q^2+1/\xi^2)^2 + q_1^2 v^2} \\ &= \frac{e^{-t/\xi^2}}{t^{d/2+1}} \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2/t + 1/\xi^2} |\hat{V}_{\mathbf{q}/\sqrt{t}}| \frac{e^{-q^2}}{(q^2/t + 1/\xi^2)^2 + q_1^2 v^2/t} \\ &\sim_{t \gg 1} \frac{\xi^6 e^{-t/\xi^2}}{t^{d/2+1}} |\hat{V}_0|^2 \int \frac{d^d\mathbf{q}}{(2\pi)^d} q_1^2 e^{-q^2} = \frac{|\hat{V}_0|^2}{2(4\pi)^{d/2}} \frac{\xi^6 e^{-t/\xi^2}}{t^{d/2+1}}, \end{aligned} \quad (B.10) $$

which corresponds to the expression reported in equation (41). Finally, for non-critical model B, we find

$$ \begin{aligned} \mathcal{R}^B(1/\xi, v) - \mathcal{F}^B(1/\xi, t, v) &= \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + 1/\xi^2} |\hat{V}_{\mathbf{q}}|^2 \frac{e^{-q^2(q^2+1/\xi^2)t}}{q^4(q^2+1/\xi^2)^2 + q_1^2 v^2} \\ &= \frac{\xi^{d+2}}{t^{d/2}} \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2\xi^4/t + 1} |\hat{V}_{\mathbf{q}/\sqrt{t}}| \frac{e^{-q^2-q^4\xi^4/t}}{q^4\xi^2/t(q^2\xi^2/t+1/\xi^2)^2 + q_1^2 v^2} \\ &\sim_{t \gg 1} \frac{\xi^{d+2}}{t^{d/2}v^2} |\hat{V}_0|^2 \int \frac{d^d\mathbf{q}}{(2\pi)^d} e^{-q^2} = \frac{|\hat{V}_0|^2}{(4\pi)^{d/2}} \frac{\xi^{d+2}}{v^2 t^{d/2}}, \end{aligned} \quad (B.11) $$

which is the expression reported in equation (42).

For a finite correlation length $\xi$, the limit $\mathcal{F}^{A/B}(1/\xi, t \to \infty, v)$ exists for any dimensionality. We can thus write

$$ \begin{aligned} \mathcal{R}^{A/B}(1/\xi, v) &= \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{q_1^2}{q^2 + 1/\xi^2} \frac{|\hat{V}_{\mathbf{q}}|^2}{[q^\alpha(q^2 + 1/\xi^2)]^2 + q_1^2 v^2} \\ &= \xi^{2-d} \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{q_1^2}{1 + q^2} \frac{|\hat{V}_{\mathbf{q}/\xi^2}|^2}{[q^\alpha(1 + q^2)]^2 / \xi^{2(\alpha+1)} + q_1^2 v^2} \\ &\sim_{\xi \to \infty} \frac{\xi^{2-d}}{v^2} |\hat{V}_0|^2 \int \frac{d^d\mathbf{q}}{(2\pi)^d} \frac{1}{1 + q^2} = |\hat{V}_0|^2 \frac{\Gamma(1-d/2)}{(4\pi)^{d/2}} \frac{\xi^{2-d}}{v^2}, \end{aligned} \quad (B.12) $$

leading to equation (43). Notice that the integral in the last line yields a finite result for $d < 2$, while it exhibits divergences at large $q$ for $d \ge 2$.
### B.3. Exact expressions for the critical recoil with model A dynamics

It turns out that the integral in the expressions (32) and (33) for the scaling function of the recoil can be calculated exactly in the case of model A with $\xi \to \infty$ when the Gaussian interaction kernel (see equation (5)) is considered. Here we report the results of such an integration, and discuss its main features.

For $d=1$, the analytic expression for the scaling function of the critical recoil turns out to be

$$ \mathcal{F}^A(0, t, v) = \frac{1}{v^2} \frac{\sqrt{1+t}-1}{\sqrt{\pi}} + \frac{1}{2v^3} \left[ e^{v^2(1+t)} \text{Erfc}(v\sqrt{1+t}) - e^{v^2} \text{Erfc}(v) \right], \quad (B.13) $$

where $\text{Erfc}(x) = 1 - (2/\sqrt{\pi}) \int_0^x du e^{-u^2}$ is the complementary error function. Accordingly, the long-time behavior of equation (B.13) is given by

$$ \mathcal{F}^A(0, t, v) \approx \frac{\sqrt{t}}{\sqrt{\pi} v^2} - \left[ \frac{1}{\sqrt{\pi} v^2} + \frac{1}{2v^3} e^{v^2} \text{Erfc}(v) \right], \quad \text{for } t \gg 1, \quad (B.14) $$

where we used the fact that $e^{v^2} \text{Erfc}(x) \approx 1/(\sqrt{\pi}x)$ for large $x$. Accordingly, upon increasing $t$, this expression diverges like $\sqrt{t}$. Note that the prefactor $1/\sqrt{\pi}$ corresponds to $C_0(1)$ given by equation (B.9). Moreover, expanding the expression in the square parenthesis in equation (B.13) around $v=0$ shows that the critical scaling function (see equation (33) with $R/\xi \to 0$) behaves like $v^{-1}$ for $v \ll 1$. This is in agreement with the expression reported in table 1.

Next, in the case of $d=2$, we find

$$ \mathcal{F}^A(0, t, v) = \frac{\ln(1+t)}{4\pi v^2} + \frac{1}{4\pi v^2} \left[ e^{v^2(1+t)/2} K_0 \left( \frac{v^2(1+t)}{2} \right) - e^{v^2/2} K_0 \left( \frac{v^2}{2} \right) \right], \quad (B.15) $$

where $K_0$ is the modified Bessel function of the second kind. Since $e^x K_0(x) \approx \sqrt{\pi/(2x)}$, the recoil exhibits a logarithmic divergence at long times. Also, the critical scaling function in equation (33) behaves asymptotically as $\propto \ln v$ for small $v$, in agreement with table 1.

Finally, in the case of $d=3$,

$$ \mathcal{F}^A(0, t, v) = \frac{1}{8\pi v^3} \left[ \frac{1 - e^{v^2(1+t)} \text{Erfc}(v\sqrt{1+t})}{1+t} - (1 - e^{v^2} \text{Erfc}(v)) \right] \\ - \frac{1}{4\pi^{3/2}v^2} \left( \frac{1}{\sqrt{1+t}} - 1 \right), \quad (B.16) $$

hence,

$$ \mathcal{F}^A(0, t, v) \sim \frac{1}{4\pi^{3/2}v^2} \left[ 1 - \frac{\sqrt{\pi}}{2v} \left( 1 - e^{v^2} \text{Erfc}(v) \right) \right] - \frac{1}{4\pi^{3/2}v^2} \frac{1}{\sqrt{t}}, \quad \text{for } t \gg 1, \quad (B.17) $$
where we used the behavior of Ercf(x) at large x reported after equation (B.14). Accordingly, the scaling function in equation (33)—and consequently the recoil— approaches asymptotically a finite value from below as $t^{-1/2}$. Again, the prefactor $-1/(4\pi^{3/2})$ agrees with the expression of $C_0(3)$ given in equation (B.9). Expanding equation (B.16) around $v=0$ for finite $t$ shows that the critical scaling function exhibits the asymptotic $\sim v^0$ behavior for small $v$, in agreement with table 1. However, taking first the limit $t \to \infty$ in equation (B.16) and then expanding around $v=0$ yields an asymptotic $\sim v^{-1}$ behavior, again in agreement with the scaling reported in table 1.

# Appendix C. Field configuration with a driven particle

In this appendix, we analyze the steady-state configuration of the average field in the reference frame co-moving with the harmonic potential. The Fourier transform of this field in the noiseless limit is given by equation (52). In particular, we consider separately the cases of model A and model B dynamics for the field, in appendices C.1 and C.2, respectively, focusing on the behavior far from the particle.

## C.1. Model A

In order to determine the behavior of the field profile $\varphi^{(ss)}(y)$ in the noiseless limit $T=0$, given by equation (52), away from the particle, i.e. for $|y| \gg R$, one can conveniently approximate the interaction potential $V(\mathbf{x})$ with that of a point-like particle, i.e. $V(\mathbf{x}) \propto \delta^{(d)}(\mathbf{x})$. We thus evaluate in Cartesian coordinates

$$ \begin{aligned} \varphi^{(ss)}(y) &= \lambda D_0 \int \frac{d^d \mathbf{q}}{(2\pi)^d} \frac{e^{i q_1 y}}{D_0 q^2 - i q_1 v} = \lambda D_0 \int_0^\infty du \int \frac{d^d \mathbf{q}}{(2\pi)^d} e^{-u(D_0 q^2 - i q_1 v) + i q_1 v} \\ &= \frac{2}{(4\pi)^{d/2}} \lambda \left( \frac{v}{D_0 |y|} \right)^{(d-2)/2} e^{-y v/(2D_0)} K_{(d-2)/2} \left( \frac{v|y|}{2D_0} \right). \end{aligned} \quad (C.1) $$

Using the asymptotic behavior of the modified Bessel function of the second type

$$ K_n(x) \approx e^{-x} \sqrt{\frac{\pi}{2x}}, \quad \text{for } x \gg 1, \quad (C.2) $$

mentioned previously in the case of $n=0$ (in the footnote<fnref n="10" /> on Page 22), one obtains

$$ \varphi^{(ss)}(y) \approx \frac{\lambda}{(4\pi|y|)^{\frac{d-1}{2}}} \left(\frac{v}{D_0}\right)^{\frac{d-3}{2}} e^{-v(y+|y|)/(2D_0)}, \quad \text{for } |y| \gg 1. \quad (C.3) $$

Note that the left tail (corresponding to $y<0$) of the steady-state configuration of the field exhibits an algebraic decay $\sim |y|^{-(d-1)/2}$. Moreover, the right tail (for $y>0$) is characterized by the lengthscale reported in equation (58), which characterizes the exponential decay of the profile upon moving away from the particle.

<footnote n="10">¹⁰ Note that equation (59) agrees with equation (55), since $e^x K_0(x) \sim \sqrt{\pi/(2x)}$ for large $x$.</footnote>
## C.2. Model B

For model B, it is convenient to introduce the dimensionless variable $\tilde{r} = 2^{2/3}rD_2^{2/3}/(3v^{2/3})$, which is connected to the lengthscale introduced in equation (38) via

$$ \tilde{r} = \frac{2^{2/3}}{3} \left( \frac{l_{v,2}}{\xi} \right)^2 . \qquad (C.4) $$

For $\tilde{r} < 1$, we define

$$ \beta(\tilde{r}) = \left(1 + \sqrt{1 - \tilde{r}^3}\right)^{1/3}, \qquad (C.5) $$

and evaluate the inverse Fourier transform of equation (52) which, for $y < 0$, is given by

$$ \varphi^{(\text{ss})}(y) = \frac{\lambda(2D_2)^{1/3} (\beta + \tilde{r}\beta^{-1}) e^{-|y|(v/2D_2)^{1/3}(\beta + \tilde{r}\beta^{-1})}}{3v^{1/3} \beta^2 + \tilde{r} + \tilde{r}^2\beta^{-2}}, \qquad (C.6) $$

while, for $y > 0$, it is given by

$$ \varphi^{(\text{ss})}(y) = \frac{\lambda(2D_2)^{1/3} (\beta + \tilde{r}\beta^{-1}) e^{-y(v/2D_2)^{1/3}(\beta + \tilde{r}\beta^{-1})/2}}{(\beta^2 + \tilde{r} + \tilde{r}^2\beta^{-2})} \left[ \cos\left(\frac{\sqrt{3}}{2}y(v/2D_2)^{1/3}(\beta - \tilde{r}\beta^{-1})\right) - \sqrt{3}\frac{\beta^2 + \tilde{r}^2\beta^{-2}}{\beta^2 - \tilde{r}^2\beta^{-2}} \sin\left(\frac{\sqrt{3}}{2}y(v/2D_2)^{1/3}(\beta - \tilde{r}\beta^{-1})\right) \right], \qquad (C.7) $$

where the argument $\tilde{r}$ has been omitted from $\beta(\tilde{r})$ for brevity. Note that, as a function of $y$, the term in the second line is characterized by spatial oscillations with a typical wavevector $k$ which vanishes as $\sim \sqrt{1-\tilde{r}}$ when $\tilde{r} \to 1$ from below. For $\tilde{r} > 1$, we additionally define

$$ \theta(\tilde{r}) = \frac{1}{3} \arctan \left( (\tilde{r}^3 - 1)^{-1/2} \right). \qquad (C.8) $$

Then, for $d=1$ and $\tilde{r} > 1$, a standard evaluation of the inverse Fourier transform of equation (52) by the means of a complex integration yields

$$ \varphi^{(\text{ss})}(y) = \frac{\lambda\xi}{2\sqrt{3}} \frac{\cos\left(\frac{\pi}{6} - \theta\right)}{\cos\theta \sin\left(\frac{\pi}{6} + \theta\right)} e^{-\frac{|y|}{\xi} \frac{2}{\sqrt{3}} \cos\left(\frac{\pi}{6} - \theta\right)}, \quad \text{for } y < 0, \qquad (C.9) $$

and

$$ \varphi^{(\text{ss})}(y) = \frac{\lambda\xi}{2\sqrt{3}} \frac{\cos\left(\frac{\pi}{6} + \theta\right)}{\cos\theta \sin\left(\frac{\pi}{6} - \theta\right)} e^{-\frac{y}{\xi} \frac{2}{\sqrt{3}} \cos\left(\frac{\pi}{6} + \theta\right)} \\ - \frac{\lambda\xi}{2\sqrt{3}} \frac{\sin\theta}{\sin\left(\frac{\pi}{6} + \theta\right) \sin\left(\frac{\pi}{6} - \theta\right)} e^{-\frac{y}{\xi} \frac{2}{\sqrt{3}} \sin\theta}, \quad \text{for } y > 0, \qquad (C.10) $$
where we dropped the argument $\tilde{r}$ in $\theta$ for brevity. Note that the $y$-integral over $\mathbb{R}$ vanishes as expected because of the local conservation of the field. In passing, we note that taking the limit $v \to 0$ ($\tilde{r} \to \infty$) in equations (C.9) and (C.10) recovers the field configuration at equilibrium. In addition, taking $\xi \to \infty$ yields the result reported in equation (62). Importantly, the qualitative behavior of the steady-state field configuration for $y > 0$ in the case of model B dynamics, which is either characterized by a superposition of two exponential functions as in equation (C.10), or by damped oscillations as in equation (C.7), holds for any spatial dimensionality $d$ of the system. In fact, these two qualitatively distinct behaviors are determined by the value of the parameter $\tilde{r}$ and follow from the analysis of the poles of equation (52), which holds for any $d$.

## ORCID iDs

Marcin Piotr Pruszczyk ID 0000-0003-2344-3333

Davide Venturelli ID 0000-0001-6952-9568

Andrea Gambassi ID 0000-0003-3450-6125

## References

[1] Einstein A 1906 *Ann. Phys.* **324** 371–81

[2] Langevin P 1908 Sur la théorie du mouvement brownien *Compt. Rendus* **146** 530–3

[3] Dhont J K 1996 *An Introduction to Dynamics of Colloids* (*Studies in Interface Science*) vol 2 (Elsevier)

[4] Larson R 1999 *The Structure and Rheology of Complex Fluids* (*Engineeringpro Collection*) (Oxford University Press) (available at: https://books.google.it/books?id=Vt9fwpf1LUC)

[5] Squires T M and Brady J F 2005 *Phys. Fluids* **17** 073101

[6] Gazuz I, Puertas A M, Voigtmann T and Fuchs M 2009 *Phys. Rev. Lett.* **102** 248302

[7] Gomez-Solano J R and Bechinger C 2014 *Europhys. Lett.* **108** 54008

[8] Gomez-Solano J R and Bechinger C 2015 *New J. Phys.* **17** 103032

[9] Jain R, Ginot F and Krüger M 2021 *Phys. Fluids* **33** 103101

[10] Jain R, Ginot F, Berner J, Bechinger C and Krüger M 2021 *J. Chem. Phys.* **154** 184904

[11] Mori H 1965 *Progr. Theor. Phys.* **33** 423–55

[12] Zwanzig R 2001 *Nonequilibrium Statistical Mechanics* (Oxford University Press)

[13] Caspers J, Ditz N, Krishna Kumar K, Ginot F, Bechinger C, Fuchs M and Krüger M 2023 *J. Chem. Phys.* **158** 1089–7690

[14] Loos S A M, Monter S, Ginot F and Bechinger C 2024 *Phys. Rev. X* **14** 021032

[15] Berner J, Müller B, Gomez-Solano J R, Krüger M and Bechinger C 2018 *Nat. Commun.* **9** 999

[16] Mayer J, Khairy K and Howard J 2010 *Am. J. Phys.* **78** 648–9

[17] Daldrop J O, Kowalik B G and Netz R R 2017 *Phys. Rev. X* **7** 041065

[18] Basu U, Démery V and Gambassi A 2022 *SciPost Phys.* **13** 078

[19] Robertson-Anderson R M 2018 *ACS Macro Lett.* **7** 968–75

[20] Ginot F, Caspers J, Reinalter L F, Krishna Kumar K, Krüger M and Bechinger C 2022 *New J. Phys.* **24** 123013

[21] Cao X, Das D, Windbacher N, Ginot F, Krüger M and Bechinger C 2023 *Nat. Phys.* **19** 1904–9

[22] Démery V and Dean D S 2010 *Phys. Rev. Lett.* **104** 080601

[23] Démery V and Dean D S 2010 *Eur. Phys. J. E* **32** 377–90

[24] Démery V and Dean D S 2011 *Phys. Rev. E* **84** 010103

[25] Démery V and Dean D S 2011 *Phys. Rev. E* **84** 011148

[26] Dean D S and Démery V 2011 *J. Phys.: Condens. Matter* **23** 234114

[27] Démery V 2013 *Phys. Rev. E* **87** 052105

[28] Gross M 2021 *J. Stat. Mech.* 063209

[29] Venturelli D and Gambassi A 2022 *Phys. Rev. E* **106** 044112

[30] Venturelli D and Gross M 2022 *J. Stat. Mech.* 123210

[31] Venturelli D and Gambassi A 2023 *New J. Phys.* **25** 093025
[32] Venturelli D, Loos S A M, Walter B, Roldán E and Gambassi A 2024 *Europhys. Lett.* **146** 27001

[33] Démery V and Gambassi A 2023 *Phys. Rev. E* **108** 044604

[34] Muzzeddu P L, Venturelli D and Gambassi A 2025 *New J. Phys* **27** 105003

[35] Hertlein C, Helden L, Gambassi A, Dietrich S and Bechinger C 2008 *Nature* **451** 172-5

[36] Gambassi A, Maciolek A, Hertlein C, Nellen U, Helden L, Bechinger C and Dietrich S 2009 *Phys. Rev. E* **80** 061143

[37] Paladugu S, Callegari A, Tuna Y, Barth L, Dietrich S, Gambassi A and Volpe G 2016 *Nat. Commun.* **7** 11403

[38] Martínez I A, Devailly C, Petrosyan A and Ciliberto S 2017 *Entropy* **19** 77

[39] Magazzù A, Callegari A, Staforelli J P, Gambassi A, Dietrich S and Volpe G 2019 *Soft Matter* **15** 2152-62

[40] Krech M 1994 *The Casimir Effect in Critical Systems* (World Scientific) (available at: www.worldscientific.com/doi/pdf/10.1142/2434)

[41] Brankov J G, Danchev D M and Tonchev N S 2000 *Theory of Critical Phenomena in Finite-Size Systems* (World Scientific) (available at: www.worldscientific.com/doi/abs/10.1142/4146)

[42] Gambassi A 2009 *J. Phys.: Conf. Ser.* **161** 012037

[43] Gambassi A and Dietrich S 2024 *Soft Matter* **20** 3212-42

[44] Hohenberg P C and Halperin B I 1977 *Rev. Mod. Phys.* **49** 435-79

[45] Reister E and Seifert U 2005 *Europhys. Lett.* **71** 859

[46] Reister-Gottfried E, Leitenberger S M and Seifert U 2010 *Phys. Rev. E* **81** 031903

[47] Camley B A and Brown F L H 2014 *J. Chem. Phys.* **141** 075103

[48] Gompper G and Hennes M 1994 *Europhys. Lett.* **25** 193

[49] Hennes M and Gompper G 1996 *Phys. Rev. E* **54** 3811-31

[50] Gonnella G, Orlandini E and Yeomans J M 1997 *Phys. Rev. Lett.* **78** 1695-8

[51] Täuber U C 2014 *Critical Dynamics: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior* (Cambridge University Press)

[52] Venturelli D, Ferraro F and Gambassi A 2022 *Phys. Rev. E* **105** 054125

[53] Démery V and Fodor E 2019 *J. Stat. Mech.* 033202

[54] Poole R 2012 *Brit. Soc. Rheol. - Rheol. Bull.* **53** 32-39

[55] Le Bellac M 1991 *Quantum and Statistical Field Theory* (Clarendon)

[56] Maes C 2020 *Phys. Rev. Lett.* **125** 208001

[57] Granek O, Kafri Y and Tailleur J 2022 *Phys. Rev. Lett.* **129** 038001

[58] Guevara-Valadez C A, Marathe R and Gomez-Solano J R 2023 *Physica A* **609** 128342

[59] Santra I 2023 *J. Phys. Complexity* **4** 015013

[60] Pei J H and Maes C 2025 *Phys. Rev. E* **111** L032101

[61] Davies B 1978 *Integral Transforms and Their Applications* (Springer) (available at: https://doi.org/10.1007/978-1-4757-5512-1)