Memory-induced oscillations of a driven particle in a dissipative correlated medium

The overdamped dynamics of a particle is in general affected by its interaction with the surrounding medium, especially out of equilibrium, and when the latter develops spatial and temporal correlations. Here we consider the case in which the medium is modeled by a scalar Gaussian field with relaxational dynamics, and the particle is dragged at constant velocity through the medium by a moving harmonic trap. This mimics the setting of an active microrheology experiment conducted in a near-critical medium. When the particle is displaced from its average position in the nonequilibrium steady state, its subsequent relaxation is shown to feature damped oscillations. This is similar to what has been recently predicted and observed in viscoelastic fluids, but differs from what happens in the absence of driving or for an overdamped Markovian dynamics, in which cases oscillations cannot occur. We characterize these oscillating modes in terms of the parameters of the underlying mesoscopic model for the particle and the medium, confirming our analytical predictions via numerical simulations.


I. INTRODUCTION
Complex media with macroscopic relaxation timescales and spatial correlations are not expected to be a source of white noise for the stochastic dynamics of a particle immersed into them. In fact, the assumption of timescale separation, underlying the use of a simple Langevin equation [1] for describing the dynamics of the mesoscopic particle, is no longer valid when the particle motion and the evolution of the medium occur over comparable timescales [2]. Moreover, in the presence of an external driving force acting on the particle, the surrounding medium to which the particle is (weakly) coupled can no longer be assumed to remain in equilibrium, as in the case of the undriven Brownian motion. In fact, the medium is expected to react to the passage of the particle, and thus to be generically out of equilibrium.
In this respect, a notable example of such complex media is provided by viscoelastic fluids [3]: their non-Newtonian behavior originates from the storage and dissipation of energy within their complex microstructure, which translate into a macroscopically long stressrelaxation time. Dragging a colloidal particle through such a fluid -as it is typically done in active microrheology experiments [4][5][6][7][8][9] -drives the medium out of equilibrium. In turn, this affects the statistics of the particle position [2]. At a coarse-grained scale, the resulting particle dynamics is often described by an overdamped generalized Langevin equation (GLE [10,11]). In this equation, the effect of the interaction between the particle and the medium is encoded in a friction kernel K(t) acting on the particle velocity as t −∞ du K(t − u)Ẋ(u) = F (X, t) + ζ(t), where F includes the forces exerted on the particle at position X, while ζ is a colored Gaussian noise.
Recently, it has been experimentally shown that viscoelasticity can give rise to oscillating modes in the overdamped motion of colloidal particles driven through the * dventure@sissa.it medium [12]. This is somewhat unexpected and noteworthy, because oscillations (which typically occur in systems with underdamped dynamics) are strictly forbidden at equilibrium, as shown, e.g., in Ref. [12]. Heuristically, one may note that integrating by parts the retarded friction in the GLE above formally renders a term t −∞ du M(t − u)Ẍ(u), where M(t) ∼ t du K(u) can be readily interpreted (if positive) as a memory-induced inertia [11], which is generally absent from the description of Markovian overdamped systems.
Memory terms in the effective evolution equation of a particle actually appear quite naturally in many physical systems, after integrating the slow degrees of freedom out of the original, microscopic dynamics in which they are coupled to those describing the tracer particle [13][14][15]. For example, a minimal model for diffusion in a thermally fluctuating correlated medium can be formulated in terms of the joint overdamped dynamics of a particle and of a scalar Gaussian field ϕ(x, t), the latter being characterized by a correlation length ξ and a finite relaxation time [16][17][18][19][20][21][22][23][24][25][26]. If the coupling between the field and the particle is chosen to be linear, then the field can be integrated out exactly, resulting into an effective evolution equation for the particle. This equation provides insight on the connection between the emerging memory kernel and the features of the original microscopic model. The system described here may be viewed as a toy model for a colloidal particle in contact with a fluid medium in the vicinity of a critical point, such as a binary liquid mixture, which displays long-range spatial correlations and long relaxation times. In this specific example, the field ϕ represents the order parameter associated with the second-order phase transition, while hydrodynamic effects and other slow variables that should be taken into account when describing real fluids [27] are neglected for simplicity. In recent years, this kind of physical systems have been the subject of various experimental investigations [28][29][30][31][32], especially concerning the emergence of effective, critical Casimir forces mediated by the fluctuations of the medium [33][34][35]. More generally, analogous spatio-temporal correlations also characterize the dynam-ics of such diverse physical systems as inclusions in lipid membranes [36][37][38], microemulsions [39][40][41], or defects in ferromagnetic systems [16][17][18][19][20][21].
In this context, it is natural to ask whether the memory kernel in the effective evolution equation of the tracer particle, which originates from the spatio-temporal correlations of the field, may give rise to oscillating modes similar to those observed in viscoelastic fluids [12] -which, instead, are primarily due to the mechanical response of the medium. In order to address this question, we consider here the simple setting of a particle driven through the medium at a constant velocity v by a moving harmonic potential, while being in contact with a scalar Gaussian field subject to an overdamped relaxational dynamics [42]. (The driving considered here can be practically realized via optical tweezers [43].) We first integrate out the field degrees of freedom, thus obtaining an effective (non-linear) equation which describes the motion of the particle in the steady state reached by the system at long times. By linearizing this equation and inspecting the analytic structure of the field-induced memory kernel, we demonstrate that damped oscillations are indeed displayed by the particle during the relaxation that occurs after it has been displaced from its steady-state position. These oscillations are confirmed via numerical simulations of the model. The simplicity of our model allows us to study in detail how the interplay between the various timescales of the system dictates the emergence of the particle oscillations, and to determine their frequency and typical decay time.
The rest of the presentation is organized as follows. In Sec. II we introduce the model and we characterize the steady state attained at long times by the particle in the moving trap. In Sec. III we analyze the relaxation of a particle initially displaced from its steady-state position, under the assumption that thermal fluctuations are negligible. The effect of these fluctuations is then assessed in Sec. IV, where we also compare our analytical predictions with numerical simulations. We finally summarize our findings and present our conclusions in Sec. V.

II. THE MODEL
We consider a particle at position Y(t) ∈ R d in d spatial dimensions, which is confined by means of a moving harmonic potential of stiffness κ, the center of which is dragged at a constant velocity v. The particle (solely described by the position Y of its center) is in contact with a correlated medium, which we model for simplicity as a fluctuating scalar order parameter field ϕ(x, t) ∈ R. The fluctuations of the latter are assumed to be characterized by the quadratic Schematic representation of the system under investigation. A particle is spatially confined by a harmonic potential (realized, e.g., by optical tweezers), the center of which is dragged at constant velocity v. The particle interacts with a thermally fluctuating order parameter ϕ (red background) according to the Hamiltonian in Eq. (4). The particle position X(t) is measured in a comoving frame of reference, chosen such that ⟨X⟩ = 0 in the steady state (see Section III for details). At time t = 0 the particle is suddenly displaced by a small amount X0 from its current position, and the ensuing relaxation is observed.

Hamiltonian [27]
where ξ = r −1/2 ≥ 0 is the correlation length, which controls the spatial range of the field correlations at equilibrium, and diverges upon approaching the critical point r = 0. The interaction between the particle and the field is chosen as [23][24][25]44] so that the system composed by the particle and the field is described by the total Hamiltonian The system under investigation is schematically described in Fig. 1. The coupling in Eq. (3) is linear and translationally invariant, while the interaction potential V (x) models the shape of the particle. We choose V (x) to be normalized so that its integral over all space is equal to one. With this normalization, the strength of the fieldparticle interaction is set only by the coupling constant λ. If λ V (x) in Eq. (4) is chosen to be positive, then field configurations are favored in which ϕ is locally enhanced, and therefore it assumes preferentially positive values in the vicinity of the particle. At the same time, the particle experiences an attractive force directed along the gradient of the field. We make the assumption that the interaction potential V (x) = V (x/R) is isotropic and characterized by a single length scale, namely the "radius" R of the particle. For example, we may choose an exponentially decaying potential where Ω d is the d-dimensional solid angle, and Γ E (z) is the Euler gamma function. The physical dimensions [ϕ] and [λ] of the field and the coupling, respectively, follow from the dimensional analysis of the Hamiltonian in Eq. (4). In units of energy E and length L, they are given by [ϕ] = E 1/2 L 1−d/2 and [λ] = E 1/2 L d/2−1 . These expressions facilitate the dimensional analysis of the quantities introduced further below.

A. Dynamics
The dynamics of the particle is here described by the overdamped Langevin equatioṅ where ν is the mobility of the particle, ϕ q = dx ϕ(x, t) exp(−iq · x), and analogously V q is the Fourier transform of V (x) [45]; finally, ξ(t) is a white Gaussian noise with zero mean and variance (k B ≡ 1) where T is the temperature of the thermal bath (see below). Similarly, we assume a purely relaxational dynamics for the field [42], i.e., Here α = 0 for a non-conserved dynamics of the order parameter ϕ, while α = 2 if ϕ is subject to a local conservation during the evolution. In this case, Eq. (8) can be cast in the form ∂ t ϕ(x, t) = −∇ · J(x, t) with a suitable current J(x, t). The two choices of α correspond to model A (α = 0) and model B (α = 2) in the classification of Ref. [27], in which we neglect the self-interaction term ∝ ϕ 4 (i.e., we consider the Gaussian approximation of these models). The particle and the field are assumed to be in contact with the same thermal bath at temperature T , so that η(x, t) is also a Gaussian white noise with zero mean and variance where D is the mobility of the field. In the absence of external dragging, i.e., for v ≡ ∥v∥ = 0, the coupled dynamics of the field and of the particle satisfies detailed balance, and therefore the stationary state is described by the equilibrium canonical distribution P eq [ϕ, Y] ∝ exp(−H[ϕ, Y]/T ) [23], with H given in Eq. (4). In view of deriving the effective dynamics of the particle, it is convenient to write Eq. (8) in Fourier space aṡ Upon setting λ = 0, Eqs. (6) and (10) reduce to a collection of non-interacting Ornstein-Uhlenbeck processesone for each of the d components of the position of the particle, and one for each of the field modes ϕ q (which form a continuum in the bulk, i.e., q ∈ R d ). These processes are characterized by the inverse relaxation timescales (see Appendix A) Accordingly, ϕ(x, t) is a medium which is correlated over both space and time, and in which the corresponding ranges are determined by ξ and τ ϕ , respectively. In particular, the relaxation time τ ϕ (q ∼ 0) of the long-wavelength modes of the field becomes arbitrarily long for model A dynamics at r = 0. The same happens for model B with generic values of r, i.e., also off-criticality, due to the presence of the conservation law [42]. These long-wavelength modes are always present in the bulk, while they are cut-off in a confined geometry such as that considered in Refs. [22,26].

B. Steady state in the comoving frame
We start by measuring the position Z ≡ Y − vt of the particle in the frame of reference that is comoving with the harmonic trap. In terms of the coordinate Z, the equations of motion (6) and (8) becomė where we introduced the translated field φ( , by translational invariance (which applies also to the white noises ξ and η). In Fourier space, these equations can be written aṡ with α q given in Eq. (13). Note that, for λ = 0, the evolution equation (16) for φ q (with v ̸ = 0) is formally the same as that for ϕ q in a fixed reference frame (i.e., with v = 0 -see, c.f., Eq. (A1)), up to a shift α q → (α q −iq·v). Accordingly, its solution in the steady state is the same as the equilibrium one reported in Appendix A, upon replacing the equilibrium correlator C   The field rearranges around the particle forming a shadow (see Eq. (21)), represented here in spatial dimensionality d = 1 (green line). The particle is subject to the attractive force F λ (Z) directed towards the shadow (and due to the field, see Eq. (15)), to the friction force which is, on average, ⟨Fν ⟩ = −v/ν, and to the restoring force Fκ = −κZ due to the harmonic trap (red parabola). The steady-state position of the particle (see Eq. (20)) results from the balance ⟨Fκ⟩ = ⟨Fν ⟩ + ⟨F λ ⟩. For λ = 0 the field and the particle are decoupled, so that F λ = 0 and the steady-state position reduces to ⟨Z⟩ = −v/γ.
where Θ(t) indicates the Heaviside theta function. We will make use of these expressions in what follows. At long times, we expect the system to reach a stationary state with ⟨Ż⟩ ss = 0 and ⟨∂ t φ q ⟩ ss = 0 in which, according to Eqs. (15) and (16), Due to the coupling between the field φ q and the particle coordinate Z, it is difficult in general to evaluate the terms φ q e iq·Z ss and e −iq·Z ss which appear in Eqs. (20) and (21). For example, in Refs. [23][24][25] this has been achieved by a perturbative expansion in increasing powers of the weak coupling λ. In particular, the properties of the nonequilibrium stationary state of the system investigated here and predicted with that approach are analyzed quantitatively in Ref. [44], also based on numerical simulations. Here we discuss only some of the qualitative features that emerge from these studies and which are relevant in the present context.
In the steady state, the average field profile ⟨φ(x)⟩ ss (obtained from the inverse Fourier transform of ⟨φ q ⟩ ss ) is enhanced in correspondence of the particle position, and is stretched in the direction opposite to the particle motion: we will refer to this field configuration as the shadow, and we represent it schematically in Fig. 2. Note that, by using Eq. (3), the term ∝ λ on the r.h.s. of Eq. (14) can be written as λν d d z∇φ(z)V (z − Z). Accordingly, for λ ̸ = 0, the particle is subject to a force that pushes it towards the maximum of the shadow. In the stationary state, this force adds up to the friction force in counterbalancing the restoring force exerted by the harmonic trap. Using perturbative arguments [44,46], one deduces that in general the field is responsible for the emergence of an additional (non-linear) friction acting on the dragged particle. Accordingly, the equilibrium position of the particle is further displaced to the left with respect to the value ⟨Z⟩ = −v/γ it would have in the absence of the field (i.e., for λ = 0 -see Eq. (20)). Note that the formation of the shadow is due to the response of the field to the passage of the particle, an aspect which is usually neglected in models used to describe the passive advection of a particle by a fluid flow [47,48].
In the present work we are primarily interested in exploring the effect of the field when the coupling λ is relatively strong, and thus we will adopt a different approach compared to that used in Refs. [23][24][25]44]. In particular, we will focus first on the noiseless limit of the dynamics, i.e., the limit in which the amplitude T of the stochastic noises ξ(t) and η(x, t) is set to zero. This allows one to determine an analytic expression of the particle trajectory X(t) for generic values of λ. The effect of thermal noise when T ̸ = 0 will then be added perturbatively in Section IV.

III. NOISELESS LIMIT
In the absence of thermal noise (i.e., for T = 0), the equations of motion (15) and (16) become deterministic and no fluctuations occur. Accordingly, ⟨Z⟩ ss = Z (ss) , ⟨φ q ⟩ ss = φ (ss) q , while φ q e iq·Z ss = φ (ss) q e iq·Z (ss) and e −iq·Z ss = e −iq·Z (ss) . Then, by using Eqs. (20) and (21) in the steady state, one readily finds that Equation (23) provides the expression of the shadow in the absence of thermal noise. As anticipated in Section I, we aim to describe the motion of the particle after it is suddenly displaced, at time t = t 0 , from the position it assumes in the stationary state. In order to do this, one can solve Eq. (16) (where η q = 0 in the limit we are interested in) by assuming that the field configuration at time t = t 0 is the one corresponding to its stationary state -i.e., φ q (t = t 0 ) = φ (ss) q is used as the initial condition of the dynamics. The resulting evolution of the field is thus given by where χ q (t) and G q (t) are the field susceptibility and response propagator introduced in Eqs. (18) and (19), respectively. Equation (22) suggests the natural change of reference frame, in which the origin of the coordinate system corresponds to Z (ss) . Accordingly, we introduce X ≡ Z − Z (ss) , so that the resting position of the particle is X = 0 in the stationary state (as depicted in Fig. 1). By substituting φ q (t) found in Eq. (24) into Eq. (15) with ξ = 0, we obtain the effective equatioṅ This non-linear equation with memory cannot be generically solved. However, further analytical progress can be made by assuming that the particle is actually perturbed by a small, sudden displacement X 0 away from its resting position, as sketched in Fig. 1. Under this assumption, it is possible to linearize Eq. (25) around X = 0, which leads (upon using Eq. (22)) tȯ for j = 1, . . . , d. Let us now introduce the memory kernel and its Laplace transformΓ j (s) = ∞ 0 dt e −st Γ j (t); in terms of these quantities, the linearized equation of motion (26) can be written in the compact forṁ We recognize Eq. (28) as the noiseless limit of an overdamped generalized Langevin equation [10,11]. By setting t 0 = 0, the solution of the latter equation with initial condition X j (t = 0) = X 0 can be conveniently expressed in Laplace space aŝ where, as in the case ofΓ j after Eq. (27),X j (s) stands for the Laplace transform of X j (t).
A. The memory kernel Γ(t) The dynamics of X j (t) is determined by the analytic structure of the functionΓ j (s) in the complex plane, which we discuss here. For later convenience, it is useful to introduce the following timescales: The first timescale τ R is the time taken by a critical field to relax over a distance of order R: this can be seen by using Eq. (13) with r = 0, q ≃ 1/R, and z ≡ 2 + α.
(We recall that R enters as a length scale in V (x), and plays the role of the radius of the particle described by V (x).) The second timescale τ v represents, instead, the time taken by the moving trap to cover a distance of order R; equivalently, τ −1 v estimates the shear rate near the driven particle [12].
By rescaling momenta as p = qR in Eq. (27) and evaluating the Laplace transform, we easily obtain where the prefactor λ 2 ν/R d has the physical dimensions of an inverse time, whilef is a dimensionless scaling function defined as Note that V y/R is in fact R-independent by construction (see, e.g., Eq. (5)). Moreover, the timescale τ κ (see Eq. (12)) which determines the relaxation time of the particle (decoupled from the field) in the harmonic trap does not enter the memory kernel, which thus describes solely the interaction between the medium and the particle. By substituting Eq. (32) into Eq. (29), one eventually findŝ where we introduced the dimensionless coupling constant and where we simplified the notation by explicitly indicating only the dependence on s of f introduced in Eq. (32). The coupling constant g can be used to quantify the effect of the interaction with the medium on the particle dynamics. We note thatX(s) in Eq. (34) satisfies the initial value theorem for Laplace transforms [49], i.e., as expected -indeed, one can check that f (s) ∼ 1/s for large s.
FIG. 3. Evolution of the particle position X(t) as a function of time t, in the noiseless limit. At t = 0 the particle is released from the position X(t = 0 + ) = X0 away from the steadystate position X = 0. The displayed symbols are obtained from the numerical Laplace inversion ofX (sc) (s) in Eq. (38), corresponding to the strong-confinement limit. The plot refers to model A at criticality in spatial dimensionality d = 1 (see Section III B), and was obtained by fixing the Weissenberg number w = 0.75 (see Eq. (42)), while varying the coupling strength g (see Eq. (35)), with X0 = 1. The so-obtained X (sc) (t) differs from the actual γX(t) only at short times t < τκ, where the former diverges (while the latter tends to γX0 -see Eq. (36)).
In order to get physical insight into the dynamics of the particle, it is convenient to consider the case in which the timescales τ R and τ κ , whichX j depends on via f , are well separated. This is actually achieved in the strongconfinement limit [46], defined as the limit in which τ κ , determined by the harmonic trap, is shorter than the typical relaxation time τ R of the field, i.e., τ R ≫ τ κ or, equivalently, (see Eq. (12)). In this limit we will focus on the dynamics occurring at times t ≫ τ κ , so that τ κ is indeed the smallest timescale in the problem. A convenient way of singling out the behavior in this temporal regime is to consider, in Eq. (29), the formal limit γ → ∞ and thuŝ the analysis of which is simplified by the fact thatX (sc) (s) depends on s only via the function f (s). As a drawback of this approach,X (sc) (s) defined above no longer satisfies the initial value theorem and, as a consequence, its inverse Laplace transform X (sc) (t) diverges in the initial temporal region t ≤ τ κ . Beyond this initial regime, however, the functions X (sc) (t) and γX(t) are expected to agree quantitatively (as discussed in, c.f., Section III D).

B. The case of model A
To make further progress with our analysis, we focus here on the one-dimensional case d = 1, with the field poised at its critical point r = 0 (further below we consider also the case r > 0). In addition, we choose an exponential interaction potential as in Eq. (39), which takes a particularly simple form in Fourier space, namely This choice renders the expressions below more amenable to analytical manipulation. In fact, the resulting memory kernel in Eq. (32) becomeŝ where we dropped the subscript j fromΓ j (s) since we are considering d = 1. In the Gaussian model A, the dynamical exponent z equals 2, so that the integrand inΓ(s) presents two simple poles in q = ±i and two additional poles in For later convenience, we parameterized these latter two poles as indicated above, with (see Eqs. (30) and (31)), and β(s) ≡ √ s τ R + w 2 . In the context of microrheology experiments conducted in viscoelastic media, one usually identifies the Weissenberg number Wi ≡ τ s /(2τ v ), where τ s is the typical relaxation timescale of the medium. For a critical field this timescale is actually provided by τ R (see Eq. (30)), and therefore the parameter w introduced in Eq. (42) above is readily identified with the Weissenberg number Wi of the system under investigation here. By using complex integration, one then finds thatΓ(s) in Eq. (40) can be expressed aŝ which can be inserted into Eq. (29) together withΓ(s) given above in order to obtain an analytical expression forX(s). The latter can then be inverted numerically to determine X(t). An example of the resulting X(t) is shown in Fig. 3, which refers to the strong-confinement limit, while a comparison with numerical simulations is presented in, c.f., Section IV B. The oscillatory character of this X(t) is clearly visible from the figure and it can be amplified by increasing the coupling strength g (a systematic analysis of this dependence is presented in the next subsection).
By inspecting Eqs. (32) and (33), we finally note that the expression ofΓ(s) for model A away from criticality (i.e., with r > 0) can be obtained from Eq. (43) by means of the substitution where quantifies the relaxation timescale of the field φ over its correlation length ξ = r −1/2 (see Eq. (13) with q ≃ ξ −1 and z = 2).

C. Relaxation in the strong-confinement limit
In the strong-confinement limit introduced in Eq. (38), the analytic properties ofX (sc) (s) are completely determined by those of the memory kernelΓ(s) and, in particular, of the associated function f (s) introduced in Eq. (32). This kernel was specialized in Eq. (43) to the case of model A at criticality (r = 0), while for r > 0 one can use the change of variable indicated in Eq. (44). The latter impliesX r=0 (s + Dr), and therefore Above we highlighted the dependence of X (sc) (t) on r via a subscript. The inverse Laplace transform ofX (sc) (s) in the previous expression is obtained, as usual, by performing the Bromwich integral along a vertical line that is on the left of the leftmost pole of the integrand in the complex plane. Accordingly, in model A, the dynamical properties of X (sc) (t) in the off-critical case r ̸ = 0 are the same as in the critical case r = 0, up to an additional exponential decay factor exp(−t/τ ξ ) (see Eq. (45)). A second remarkable feature ofX (sc) (s) is that it depends on s only via the combination s τ R -see Eqs. (32) and (33), and the definition of τ R in Eq. (30). Taking the inverse Laplace transform ofX (sc) (s) as in Eq. (46) and changing the integration variable as s ′ ≡ sτ R , it follows that We deduce that rescaling s τ R → s inX (sc) (s) simply corresponds to measuring time t in units of τ R . Furthermore, the explicit dependence on R ofΓ(s), which occurs in Eq. (32) only via R/ξ, is lost at criticality (ξ → ∞). Accordingly, the resultingX (sc) (s) eventually depends only on the pair of parameters (w, g) -see Eqs. (35) and (41) to (43). Let us then focus on the analytic structure ofX (sc) (s) in the complex plane s ∈ C, for r = 0 and with τ R ≡ 1. First, from Eq. (41) we infer the presence of a branch cut along the real axis for Re{s} < −w, as shown in Fig. 4a. The exact position of the poles ofX (sc) (s) in Eq. (38) cannot be determined analytically; however, they are easily found numerically. Indeed, the plot of Re{X(s)} in the complex plane, shown in Fig. 4a as a colormap, reveals the presence of a pair of complex conjugate poles in s ± = −δ ± iΩ, with Ω ≥ 0. The red dashed line in the plot indicates the trajectory of these poles upon varying g at fixed w < 1. The poles appear for a small g = g * > 0, in the vicinity of the origin s = −w of the branch cut, and have a vanishing imaginary part Ω = 0; upon increasing g, they depart from the branch cut and acquire a nonzero imaginary part Ω > 0. As g is further increased, the two poles eventually move to the right of the branching point (i.e., |δ| < |w|), and thus they become the dominant singularities. The presence of a dominant complex pole in the analytic structure ofX (sc) (s) implies the emergence of an oscillatory behavior of X (sc) (t) at long times, with frequency Ω (see Appendix B for additional details). These are the oscillations featured in Fig. 3, where we plotted X (sc) (t) (obtained via numerical inversion of the analytical solution forX (sc) (s)) for increasing values of the coupling strength g, while keeping w fixed. Figure 4b shows the oscillation frequency Ω as a function of the values of the parameters (w, g). Within the white region of the plot there are no poles, and thus no oscillations occur. Even for small values of g, instead, complex poles appear and oscillations are seen to develop as soon as w ≳ 1; moreover, Ω is in general an increasing function of w, for any fixed value of g. We recall that w = Wi measures the ratio between the relaxation time of the medium and the timescale τ v set by the moving trap (see Eq. (42)). This suggests a way to rationalize the "dynamical phase diagram" in Fig. 4b. Indeed, for small g and sufficiently large values of the dragging speed v ∝ w (see Eqs. (31) and (42)), the field is no longer able to quickly rearrange around the instantaneous position assumed by the particle at a given time: the non-Markovian interplay between the dynamics of the particle and the field shadow is then at the origin of the complex oscillatory behavior of X(t). The effects of this interplay become increasingly prominent upon increasing the coupling strength g, so that at large g one observes an oscillatory behavior even for w ≲ 1.
Conversely, these oscillations are increasingly damped upon increasing the dragging speed, i.e., for w ≫ 1. This is shown in Fig. 4c, where we plot the real part δ of the dominant pole as a function of (w, g) -indeed, δ controls the long-time exponential decay of X(t) (see Appendix B). To understand the damping at large w ∝ v, we first note that the shape of the shadow φ (ss) (x) is given by the Fourier transform of φ (ss) q in Eq. (21): upon inspection, the latter shows that the amplitude of the shadow itself decreases upon increasing v [44]. This is expected, since the finite relaxation time of the field φ does not allow φ to react instantaneously to the passage of the particle, and thus at a very large speed v the shadow cannot build up at all. The damping of the oscillations at large values of w thus simply reflects these facts.
We emphasize that no poles emerge inX(s) within the white region in the (w, g)-plane in Figs. 4b and 4c. Correspondingly, the long-time behavior of X(t) in that region is determined solely by the branch cut (see Fig. 4a): as we recall in Appendix B, this generically implies that X(t) decays monotonically as X(t) ∼ t −a exp(−bt), for some positive constants a and b (see Eq. (B4)). Conversely, upon increasing g far beyond the values that Fig. 4a refers to, the real part −δ of the poles s ± eventually becomes positive. This would imply an unbounded (oscillatory) growth of X(t) at long times (see Appendix B for details), and thus it signals the breakdown of the linear-response approximation within which such solution has been derived.

D. How generic are these oscillations?
Beyond the strong-confinement limit discussed in Sec. III C above, i.e., upon decreasing the value of the parameter γ, new poles eventually appear in the complex-s plane shown in Fig. 4a. Although the precise value of Ω at a certain point (w, g) of the plane is in general modified compared to the value it has in the strong-confinement limit ρ ≫ 1 (see Eq. (37)), we find that the oscillatory nature of the solution X(t) persists, within the same range of values as in Fig. 4b, down to ρ ≳ 1. Note that, after rescaling s ′ ≡ s τ R in Eq. (34), the latter readŝ showing (as expected) that the strong-confinement limit becomes increasingly accurate as ρ ≫ 1 -compare with Eq. (38). Moreover, Figs. 3 and 4 (together with the numerical simulations presented in, c.f., Section IV B) show that X(t) typically decays to zero on a scale of a few tens of τ R . As a result, even for ρ ≲ 1, the strongconfinement limit well approximates the behavior of X(t) at times t > τ κ = γ −1 . Indeed, by taking the inverse Laplace transform ofX j (s) in Eq. (34) and calling z ≡ st, one obtains where the integration is intended along the Bromwich contour as in Eq. (46). The term z/(γt) at the denominator can be safely neglected as soon as γt ≫ 1, yielding in fact X j (t) ≃ X (sc) j (t)/γ (see Eq. (38)). Away from the critical point (i.e., for r > 0), the damped oscillations of X(t) persist, but they are additionally suppressed by the exponential factor exp(−Drt) = exp(−t/τ ξ ) (see Eqs. (45) and (46)). Taking into account all the trends highlighted above, we expect that the oscillatory behavior of X(t) is maximally amplified within the timescale window τ κ < τ R < τ ξ , where the second inequality corresponds to requiring ξ > R (see Eqs. (30) and (45)).
Note that increasing the trap strength κ has the effect of both increasing γ (thus pushing the system further into the strong-confinement regime), and decreasing the effective coupling g and therefore reducing the amplitude of the oscillations (see Figs. 4b and 4c). Accordingly, oscillations generically develop at intermediate values of κ, while they vanish both at very large and very small values of κ. This was also the case in experiments performed on colloidal particles dragged in viscoelastic media (see Ref. [12] and Fig. 5 therein). Figure 4 additionally confirms that no oscillations occur if the trap is not dragged, i.e., for v = 0 (hence w = 0). This was also the case in the experiments of Ref. [12] involving a viscoelastic medium (see Fig. 3 therein). This fact also agrees with the analytical and numerical results of Ref. [23], where the relaxation towards equilibrium of a trapped particle in contact with a near-critical Gaussian field was investigated perturbatively in the coupling λ. In particular, it was found that ⟨X(t)⟩ decreases algebraically upon increasing time t for model A at criticality, and generically for model B. For completeness, in Appendix C we reconsider this problem within the noiseless but nonperturbative approach presented in this Section, and we re-derive the exponents of the long-time algebraic decay of ⟨X(t)⟩ originally reported in Ref. [23].
Finally, the emergence of oscillations in viscoelastic media reported in Ref. [12] was rationalized therein in terms of the stochastic dynamics of an underdamped harmonic oscillator. In fact, it was shown that such a simplified model (with a positive, memory-induced mass term) is able to reproduce quantitatively the oscillations displayed at long times by the dragged colloidal particle. While such an effective model turns out to be inappropriate in our case due to the non-analytic behavior of the memory kernel, in Appendix D we discuss in detail the comparison between our model and the theoretical description of viscoelastic fluids used in Ref. [12]. In particular, it turns out that the memory kernels Γ(t) emerging in the present case and in viscoelastic media appear to be both negative at long times t, confirming that the negative response of the surrounding medium (whose origin in our model has been clarified in the previous Sections) is actually essential for the emergence of the oscillating modes exhibited by the overdamped particle.

IV. EFFECTS OF THERMAL FLUCTUATIONS
Thermal fluctuations act on the field and the particle, via the noise terms ξ(t) and η q (t) in Eqs. (15) and (16), whenever T ̸ = 0. The presence of thermal noise represents an obstacle to the analytical derivation of the time-dependent relaxation of the particle, because it modifies the steady-state average of both the particle position ⟨Z⟩ ss in Eq. (20), and the field profile ⟨φ q ⟩ ss in Eq. (21). Once incorporated into the effective equation of motion of the particle, the field-induced fluctuations turn out to be non-Gaussian, as we will verify shortly; in order to account for them, we shall resort below to a perturbative expansion in the coupling constant λ (as previously done in related investigations of this model [22][23][24][25]44]). We emphasize that the (noiseless) effective equation (25) is actually non-perturbative in λ, and so is its solution in Eq. (34). Expanding the dynamics for small λ is just a computational tool to take fluctuations into account analytically, but the qualitative conclusions we reach are valid beyond the perturbative regime, as we confirm in Section IV B by using numerical simulations.

A. Weak-coupling approximation
The effective equation (25) in Section III was determined first by choosing the shadow state in Eq. (23) as the initial condition for the field φ at time t = t 0 , and then by moving to a reference frame in which the resting position of the particle corresponds to X = 0. In this Section we adopt a different strategy: instead of explicitly determining the stationary shadow configuration (which is difficult in the presence of thermal fluctuations), we first solve for φ q (t) as we did in Eq. (24), but we impose the flat initial condition φ q (t = t 0 ) = 0 at the initial time t 0 , and we take into account the contributions due to the noise. Plugging the result into Eq. (15) then yieldṡ where the field susceptibility χ q (u) was defined in Eq. (18), and where we introduced the Gaussian colored noise which has zero mean and correlator C q (t) (see Eq. (17)). Although we did not specify the stationary shadow configuration of the field (see Fig. 2) as the initial condition of its evolution, one can convince oneself that such a configuration is inevitably recovered by taking the limit t 0 → −∞, since it coincides with the nonequilibrium steady state of the system. The leading correction to the average particle position ⟨Z⟩, due to thermal fluctuations, can then be calculated from Eq. (50) by following the steps detailed in Appendix E. This eventually yields In analogy with the derivation in Section III, we now change reference frame to X ≡ Z − ⟨Z⟩ ss in Eq. (50), we take the average over thermal fluctuations, and we linearize the resulting equation. This way we find the evolution of the average position ⟨X(t)⟩ to be given by which is formally the same as Eq. (28), but where X j (t) is replaced by ⟨X j (t)⟩, the initial time t 0 is set to −∞, and the memory kernel is replaced by As expected, compared to the memory kernel for the noiseless case in Eq. (27) -which includes only the first term on the r.h.s. of Eq. (55) -the present one involves a second term ∝ C q (t) due to thermal fluctuations. Note that the integration in the variable u in Eq. (54) runs from −∞, and this fact prevents a direct solution of the equation of motion by using the Laplace transform [50]. However, in order to determine the response of the average particle position to a sudden displacement X 0 imposed at time t = 0 from its stationary value ⟨X⟩ = 0, one can look for a solution ⟨X(t)⟩ of Eq. (55) with ⟨X(t)⟩ ≡ 0 for t < 0, and ⟨X(t)⟩ = X 0 at t = 0. In this way, ⟨X(t)⟩ for t > 0 follows immediately from a Laplace transform as in Eq. (29), with ⟨X j (s)⟩ in place ofX j (s), and with the memory kernel Γ j (t) given by the new expression in Eq. (55). In this case, an expression of the functionΓ j (s) in closed form (such as the one found in Section III B in the noiseless limit) cannot be obtained. In spite of this complication, studying the strong-confinement limit provides already valuable information concerning the main effects of thermal fluctuations. In fact, in Section III D it was shown that this limit actually captures the particle evolution for times t > τ κ . Proceeding as in Eq. (38), we then inspect the formal limit γ → ∞, which has the effect of suppressing the term proportional to the field correlator C q (t) in the memory kernel given in Eq. (55). Accordingly, this results into where the function f (s) is the same as in Eq. (33) upon replacing Since the role of V q is essentially that of providing a largemomentum cutoff for q ≳ 1/R [23,25], we conclude that V q represents an effective renormalization of the particle radius R, which is replaced by a combination of R and the thermal length appearing in Eq. (57). Note that l coincides with the mean squared displacement of the particle in its harmonic trap due solely to thermal fluctuations. For instance, a choice of V q as in Eqs. (5) and (39) yields |V q | 2 ≃ 1 − 2q 2 R 2 for small q, so that | V q | 2 ≃ 1 − 2q 2 (R 2 + l 2 ), and therefore R is effectively renormalized as R → √ R 2 + l 2 .

B. Numerical simulations
In this section we present and discuss the results of numerical simulations of the system in one spatial dimension, which confirm our analytical predictions, also beyond the noiseless limit presented in Section III and the perturbation theory discussed in Section IV. In particular, the numerical data are obtained via a direct integration of the Langevin equations (6) and (8) for the particle and the field, respectively. The latter is evaluated by discretizing the field over a regular lattice with spacing a ≪ R, similarly to Refs. [23,25,51]. The coupled stochastic differential equations are then integrated by using a stochastic Runge-Kutta algorithm, as described in Ref. [52] (which is suited for investigating also cases with an explicitly time-dependent external drag).
The field is initially prepared, at time t = −T , in the flat configuration ϕ(x, t = −T ) = 0; the harmonic potential which traps the particle is dragged for a certain time T until the system reaches its steady state, in which the average particle position stops evolving in the comoving frame of reference. At time t = 0, the particle coordinate is suddenly displaced by an amount X 0 and its relaxation is recorded. Since the actual position of the particle at time t = 0 − depends on the realization of the noise, it fluctuates. Accordingly, the result of this displacement is equivalent to extracting the initial particle position at time t = 0 + from a distribution that is the same as the one in the steady state, but shifted in space by an amount X 0 . We repeat the whole process (including thermalization) several times, and we finally take the average over the various realizations. Simulations are performed with periodic boundary conditions in order to approximate the behavior of the particle in the bulk. The lattice extension L is chosen sufficiently large so as to avoid stirring effects: in fact, a particle dragged along a ring of finite length L soon generates spurious field currents, which in general modify the particle statistics. The value of the particle displacement X 0 is chosen within the linear-response regime, which is verified a posteriori by comparing simulations performed for various (small) values of X 0 , checking that the corresponding average particle trajectories ⟨X(t)⟩ collapse onto each other upon rescaling their amplitude by X 0 . Figure 5 presents the results of the numerical simulations described above. In particular, Fig. 5a corresponds to the case of critical model A, which we studied analytically in Section III B. For various values of the the drag velocity v (which determines the value of the Weissenberg number w indicated in the plot, see Eq. (42)), we plot ⟨X(t)⟩ (solid line) of a particle that is initially displaced from its steady-state position by an amount X 0 , as a function of the time t elapsed from the displacement. These numerical curves are compared with our analytical prediction (symbols), which is obtained by numerical Laplace inversion of Eqs. (34) and (43), and in which the particle radius R is replaced by the effective radius (R 2 + l 2 ) 1/2 to account for thermal fluctuations (see discussion at the end of Section IV A). The plots show an overall agreement within the entire time range, including the fast initial decay displayed at short times. In general, this decay develops over a timescale t ∼ τ κ , followed by an oscillating behavior which persists over a few tens of τ R . Following our discussion in Section III D, the latter region t > τ κ in Fig. 5a is essentially described by the strong-confinement limit. This limit turns out to describe accurately the numerical data even when the choice of parameters is not strictly into the strong-confinement regime ρ ≫ 1, as shown in Fig. 5a (see caption), which corresponds to τ κ = γ −1 = 1, τ R = 1, and therefore ρ = 1. This fact confirms the expectation that the phenomenology described by the dynamical phase diagram presented in Fig. 4 actually carries over moderately beyond the strong-confinement limit.
Our previous discussion in Section III revealed that the behavior of the noiseless model is completely determined by fixing the dimensionless numbers g = λ 2 /(κR), ρ = γR 2 /D, and w = Rv/(2D) -see Eqs. (35), (37) and (42), here specialized for model A in spatial dimension d = 1. Thermal fluctuations are, instead, perturbatively quantified by the ratio l/R of the thermal length l (see Eq. (58)) to the particle radius R -see Section IV A. Note that the effective particle dynamics at criticality r = 0 has been written in the previous Sections in terms of n = 10 physical variables (i.e., X, t, κ, R, ν, λ, D, v, X 0 , and T ), but only k = 4 distinct physical units (i.e., mass, length, time and temperature). The physics of the model is thus actually captured by the mutual dependence of the n − k = 6 dimensionless parameters as suggested by dimensional analysis [53]. In the simulations, we chose ρ = 1 and a large effective coupling g ≃ 25, while we varied the Weissenberg number w within the range 0 < w ≲ 1. We finally added thermal fluctuations of moderate strength by tuning the noise temperature T so that l/R ≃ 10 −1 − 10 −2 . This choice facilitates the numerical computation, and we detected no significant qualitative change at higher temperatures.
In Fig. 5b we show the results of simulations analogous to those presented in Fig. 5a, but for a field that evolves FIG. 6. Evolution of the average position ⟨X(t)⟩ of the particle which is suddenly displaced at time t = 0 by an amount ±X0 from its actual position in the steady state. The solid lines in both the main panel and the inset are obtained from the numerical simulation of the system with model A dynamics in spatial dimension d = 1. The inset refers to the case X0 = 1, so that the corresponding behavior is captured by the linear-response prediction in Eqs. (34) and (43) -the symbols correspond to its numerical Laplace inversion, and the scales on the axes are the same as in the main plot. This implies, inter alia, that ⟨X(t)⟩ starting from +X0 is the opposite of ⟨X(t)⟩ starting from −X0. The main plot, instead, refers to X0 = 10, which turns out to be beyond the linear regime. In fact, the relaxation occurring from the initial values +X0 and −X0 are no longer related by the symmetry highlighted above. The remaining simulation parameters are the same as in Fig. 5a, with w = 0.6.
according to the conserved dynamics prescribed by model B. In addition, we chose here a finite correlation length ξ ≃ 2R, so that the system is off criticality. Although analytical predictions cannot be derived in closed form for model B, the simulations show a behavior similar to that of model A. This is interesting in view of possible experimental investigations of the effects qualitatively predicted in this work, because off-critical model B more realistically represents, e.g., the case of colloidal particles immersed in a binary liquid mixture [30][31][32] (still assuming that hydrodynamic effects are negligible). Note that the parameters w, ρ and the ratio l/R used above are not far from those realistically achievable in experiments [54]; the magnitude of g depends, instead, on the specific mechanism that couples the medium with the particle and, from our discussion, one expects the overall effect to be enhanced if g can be made large in an experimental realization. Note also that the interaction potential V (x) in Fig. 5b is chosen to be Gaussian with variance R, rather than exponential as in Eq. (5). The overall qualitative behavior is thus shown to be robust against changing the details of V (x), as expected [23].
Finally, we can use the numerical simulations to explore qualitative features that are not captured by the linear-response analysis. In particular, our analytical so-lution in Eq. (29) depends linearly on the initial particle displacement X 0 , meaning that the evolution of ⟨X(t)⟩ after displacing the particle in the steady state by +X 0 is expected to be the opposite of that of ⟨X(t)⟩ after a displacement −X 0 . This is indeed the case in our numerical simulations performed at small X 0 , as we show in the inset of Fig. 6 (it is also mostly the case in the experiments of Ref. [12] -see Fig. 3 therein). However, the asymmetry between the two evolutions is expected to emerge upon increasing X 0 , as it is clearly shown in Fig. 6. This asymmetry is a consequence of the non-linearity of the field-particle coupling, and therefore of the effective evolution equation for the particle position.

V. SUMMARY AND CONCLUSIONS
In this work we considered an overdamped Brownian particle dragged at constant velocity by a harmonic trap through a correlated medium, modeled here by a scalar Gaussian field with an overdamped Langevin dynamics. We have demonstrated that, when displaced from its position in the steady state, the resulting average position of the particle can exhibit oscillations during relaxation, in spite of the system dynamics being overdamped. This is reminiscent of the oscillatory modes recently observed with colloidal particles dragged through a viscoelastic fluid [12], except that the medium considered in this work is not viscoelastic. Accordingly, we have shown that oscillating modes can be found in overdamped media characterized by spatial and temporal correlations, which is typically the case for systems close to a second-order phase transition -such as those involved when studying defects moving within spin systems [17], or colloidal particles immersed in binary liquid mixtures close to the critical point of their demixing transition [35].
In particular, in Section III we first neglected thermal fluctuations and we derived an analytic solution of the effective equation of motion for the coordinate X(t) of the particle, within the linear-response approximation, and when the particle is suddenly displaced at time t = 0 from its steady-state position -see Eq. (29). This approximation involves the field-induced memory kernel Γ(t) which appears into the effective equation (28) of the particle, once the field coordinate has been integrated out. We then focused on the case d = 1 of model A dynamics, and we characterized the analytic structure of the memory kernelΓ(s) in Laplace space (see Eq. (43) and Fig. 4), along with its implications for the dynamics of X(t). In particular, it turns out that ⟨X(t)⟩ generally exhibits oscillations if the relaxation timescale of the field τ R (over distances of the order of the particle size R, see Eq. (30)) exceeds the typical timescale τ κ set by the harmonic trap (see Section III C). These oscillations are damped (and eventually vanish) at large values of the trap strength κ, and whenever the correlation length of the field ξ ≪ R (i.e., far from the critical point).
Thermal fluctuations were then reinstated into the prob-lem by using a perturbative expansion in the field-particle coupling λ -see Section IV. Their main effect on the late-time particle dynamics is a renormalization of the particle radius R by its thermal mean squared displacement l in the harmonic trap (see Eq. (58)), while the qualitative features of ⟨X(t)⟩ remain the same as in the absence of the noise. The accuracy of our analytic predictions was tested via numerical simulations in Section IV B, finding good agreement (see Fig. 5). However, simulations can also be used to explore a range of parameters that are in principle out of reach of our analytical predictions. For example, in Fig. 5b we showed that the qualitative features of ⟨X(t)⟩ obtained by using a conserved field dynamics, i.e., model B, are similar to those of model A. Moreover, such features are robust against changing the particular shape of the field-particle interaction potential V (x) (see Eq. (3)), as expected. In addition, by choosing a sufficiently large value of the initial particle displacement X 0 , we can go beyond the linear-response approximation under which our analytical predictions were derived. In Fig. 6, the actual non-linearity of the field-particle coupling causes an asymmetry between the response of the system to a +X 0 or −X 0 initial particle displacement.
Our work opens the possibility of observing oscillatory modes with colloidal particles in near-critical binary liquid mixtures, in a fashion similar to that described in Ref. [12]. This type of systems is already accessible experimentally [28][29][30][31][32], and correlation lengths of the order of microns (which is the typical size of a colloidal particle) can nowadays be obtained by using, e.g., micellar solutions. It would be therefore desirable to test if the phenomena described in this work can also be observed in experiments, at least qualitatively. In fact, a more quantitative description requires to incorporate hydrodynamic effects into the model, which is left for future investigations.
Among the various additional aspects of the dynamics of the system considered here, its stochastic thermodynamics turns out to be particularly rich [44]. Moreover, we note that other interesting features displayed by viscoelastic fluids -such as those observed in the recoil experiments performed in Ref. [7] -are found to emerge also within the minimal model for correlated (but not viscoelastic) media studied here. These interesting similarities will be explored in future works. Finally, it would be interesting to consider the case in which the field interacting with the tracer particle is active: this might provide a model of a nonequilibrium bath made of active particles, such as those recently investigated in Refs. [55][56][57][58].

ACKNOWLEDGMENTS
We thank Clemens Bechinger and Félix Ginot for useful insights. DV would like to thank Guido Giachetti and Ignacio A. Martínez for fruitful discussions. We also thank Sarah A. M. Loos, Édgar Roldán, and Benjamin Walter for collaboration on related topics. AG acknowledges support from MIUR PRIN project "Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST)" n. 201798CZL.
Appendix A: Dynamics of the free field The Langevin equation (10) for the field reads, for where α q was defined in Eq. (13), and the noise correlations are given in Eq. (11). This equation is the same as that for the coordinate ϕ q of an Ornstein-Uhlenbeck particle [59]. Accordingly, by setting for simplicity ϕ q (t 0 ) ≡ 0 (a choice which is inconsequential in the steady state which we focus on in this work), one can easily derive [42] where (A3) is the free-field correlator. By formally taking the limit t 0 → −∞ in Eq. (A3), one obtains the equilibrium correlator which is a function of the time difference t = s 2 − s 1 only. The response function G  q (t) of the free field are usually defined [42] as These quantities are related to the equilibrium correlator C where Θ(τ ) is the Heaviside theta function. The features of the long-time behavior of a function f (t) can be inferred from the analytic structure of its Laplace transformf (s). Relations of this type are referred to as Haar's Tauberian theorems in the mathematical literature [60]. In this Appendix we recap and summarize some useful related results which are applied in Section III.
Simple poles. -Consider, first of all, the textbook case [49] in whichf (s) is a meromorphic function in the complex s plane:f where g(s) is an analytic function, and the n poles {s 1 , s 2 , . . . , s n } are located at s j = −δ j + iΩ j . We order the poles so that 0 ≤ δ 1 < · · · < δ n . By using the Cauchy residue theorem we then easily obtain where in the last step we retained the dominant term at long t > 0. This shows that the rightmost pole s 1 of f (s) (i.e., the closest to the imaginary axis) determines the long-time behavior of f (t), which exhibits damped oscillations with frequency Ω 1 whenever s 1 has a nonzero imaginary part. Branch cuts. -Next, assume thatf (s) is no longer meromorphic, but rather displays a branch cut with branch point s 0 (for instance, in Fig. 4a the branch cut develops along the real axis for Re{s} < −w, with s 0 = −w). In this case, we can generally expandf (s) around the branch point s 0 aŝ for some (possibly non integer) λ j , and take the inverse Laplace transform term by term as in Eq. (46) to obtain In the presence of poles alongside the branch cut (as in Fig. 4a), the contribution in Eq. (B4) simply adds up to that in Eq. (B2). Again, the long-time behavior of f (t) is determined by the rightmost among the poles s j and the branching point s 0 . Algebraic decays. -We describe for completeness the case in which f (t) does not exhibit an oscillatory behavior, but rather an asymptotic algebraic decay of the form where µ ≥ 0, and t c is a crossover time [61]. The strategy to obtain the corresponding Laplace transform is to divide the integration domain aŝ The first term is regular, i.e., it can be expanded in a power series containing only integer powers of s. The integration in the second term of Eq. (B6) can be further split into [62] ∞ tc dt e −st t −µ = The first term on the r.h.s. of Eq. (B7) can be shown to involve both a regular and a non-regular part. To see this, we expand the exponential in power series and integrate term by term to find (B8) For µ not integer, the second term in the series above is regular. The first term ∼ s µ−1 , together with the last term in Eq. (B7), reconstructs the Euler gamma function via its integral representation (for z ̸ = 0, −1, −2, . . . [63]) If instead µ ∈ N + , let us first set µ = p + ε, with p ∈ N + and ε ≪ 1, and isolate the diverging term in Eq. (B8) as where we used x ε = e ε ln x ≃ 1 + ε ln x. Note that the remaining infinite series may still produce terms proportional to s m ln s, with m > p − 1, but these are subleading for small s. Including all the terms in Eq. (B6) we thus finally get wheref r (s) is the regular part off (s). Thus, to unveil a long-time asymptotic power-law decay of f (t) as in Eq. (B5), one can expand its Laplace transformf (s) in series for small s, and check for the presence of a term ∼ ln(st c )s µ−1 (with µ a positive integer), or ∼ s µ−1 (with µ ≥ 0 not integer). An application of these last relations is presented in Appendix C.
We consider here the problem of the relaxation towards equilibrium of a particle in contact with a scalar Gaussian field, in a fixed harmonic trap, and which is subject to a small initial displacement X 0 ̸ = 0 at time t = 0. Indeed, in the absence of external dragging the steady state reached by the system at long times (see Section II B) is actually an equilibrium state [23].
The problem is analogous to the one we analyzed in Section III upon setting v = 0, and thus the solutionX j (s) for T = 0 (noiseless case), and within the linear-response approximation, is given by Eq. (29). In particular, the memory kernel in Eqs. (32) and (33) reduces tô At the critical point r = 0 of the medium one has α q = Dq z (see Eq. (13)), so that, using polar coordinates and changing variables to y = Dq z /s, one findŝ (C2) Here c d is a numerical constant accounting for the integration over the angular variables, while in the last step we expanded the expression for small s by using the normalization condition V q = 1 + O(q) of the interaction potential. Expanding the denominator of Eq. (29) in a geometric series now giveŝ Comparing with Eq. (C2), we deduce that the power series ofX j (s) contains a term ∼ s d/z , which is nonregular whenever the ratio d/z is not integer. From our discussion in Appendix B, this corresponds to an algebraic asymptotic decay of X j (t) ∼ t −µ , with µ = 1 + d/z. Note that the terms ∼ s nd/z contained in the series of Eq. (C3) may also be non-regular, but they correspond to subleading algebraic contributions to X j (t) at long times.
The powers of the algebraic decays of X j (t) agree with those previously found in Ref. [23] -see Eqs. (33) and (34) therein. The approach used in Ref. [23] in order to derive these predictions differs, however, from the one used here and in Section III: in particular, the former relies on a weak-coupling expansion for small λ, while it assumes T ̸ = 0 and it does not require the linear-response approximation (thus it also allows the investigation of an intermediate nonlinear dynamical crossover, see Ref. [23] for details). Although limited to the noiseless case T = 0, the analysis presented here and in Section III of the present work makes, instead, no assumption concerning the magnitude of λ. This suggests that the exponents of the algebraic decays determined above may in fact be nonperturbative in λ, as was conjectured in Ref. [23] based on the evidence provided by numerical simulations. The underdamped oscillations of a colloidal particle dragged through a viscoelastic fluid reported in Ref. [12] have been described therein in terms of the linear, generalized Langevin equation is a stochastic correlated noise term with vanishing average. As in Section III, X(t) denotes the component of the particle position along the direction of the trap displacement, in the comoving frame of reference. Using that both ⟨X(t)⟩ and ⟨Ẋ(t)⟩ vanish for t < 0 and integrating by parts, one can check that Eq. (D1) is formally equivalent to the effective equation in (54) upon identifying γ ∞ ≡ 1/ν, and Γ(t) ≡ −νK ′ (t). (D2) Finally, the memory kernel is assumed in Ref. [12] to be of the form where γ 0 > γ ∞ is the zero-frequency friction coefficient, while γ i and τ i are phenomenological parameters introduced to fit the experimental data. The kernel in Eq. (D3) reduces to that of a Jeffrey's fluid [64] for γ i ≡ 0, which is shown in Ref. [12] to appropriately describe the particle dynamics in a static trap, i.e., for v = 0. For v > 0, instead, one can rationalize the experimental data by considering one or more additional relaxation timescales τ i > τ , weighted as in Eq. (D3) by some suitable coefficients γ i < 0. It actually turns out that an increasing number of pairs (γ i , τ i ) is needed for fitting the experimental data upon increasing the dragging velocity v, and hence the Weissenberg number [65].
Integrating by parts the second term on the l.h.s. of Eq. (D1) and introducing M(t) ≡ − . The parameters of the former kernel were chosen as in Ref. [12] (see the main text). For illustrative purposes, in the latter kernel in Eq. (D11) we set τ ξ = ∞ (i.e., we considered the field at criticality), while we fixed the parameters τR and τv by matching the leading behavior of the former Γ(t) at short times, and the crossing point of the horizontal axis. Both curves eventually approach zero from below as t → ∞.
cast Eq. (D1) in the form (D4) At long times, the l.h.s. may be further approximated as One normally finds m < 0 when γ i ≡ 0 (corresponding to exponentially decaying solutions for ⟨X(t)⟩), while an appropriate choice of the coefficients γ i < 0 can render m > 0, i.e., a bona-fide inertia which may explain the emergence of oscillations within the system. An oscillating behavior of ⟨X(t)⟩ with frequency is then expected for γ 0 < 2 √ mκ. It is interesting to check if the analogy with a harmonic oscillator holds for our model as well, via the mapping in Eq. (D2). Considering for instance critical model A (see Section III B) and using Eqs.
as a function of the Weissenberg number w (see Eq. (42)). This effective mass is plotted in Fig. 7a and it appears to be always positive, while it increases and diverges upon reducing the value of w towards zero. This behavior can be rationalized by comparison with the dynamical phase diagram in the strong-confinement limit shown in Fig. 4b. In fact, we note that Ω ≃ 0 in the latter as soon as the complex poles appear for small values of w and g, after which Ω is a growing function of w (for any value of g). Inverting Eq. (D8) yields in fact, consistently, m ∼ 1/Ω 2 -i.e., heavier objects oscillate more slowly.
It is now tempting to use the condition stated in Eq. (D8) in order to predict the boundaries within the phase diagram in Fig. 4b. In the strong-confinement limit of critical model A, however, the argument outlined above renders a friction coefficient γ 0 in Eq. (D4) equal to which is negative for all values of w. Accordingly, this indicates that the approximation used in Eq. (D7)corresponding to keeping only the first order term in the expansion for small s in Eq. (D5) -is no longer accurate in our case. In fact, it turns out that the nontrivial analytic structure of the memory kernel in Eq. (43) (see Fig. 4) prevents us from simply expandingΓ(s) (and hencê M(s)) in a Taylor series around s = 0. In order to check this, one can numerically invertX(s) in Eq. (29) after replacingΓ(s) by its n-th order Taylor expansion: the amplitude of the corresponding approximation to X(t) turns out to diverge upon increasing t (unlike the actual solution, which is expected to be bounded).
In conclusion, contrary to the phenomenological model presented in Ref. [12] -which explains the origin of the observed oscillations in terms of the emergence of an effective harmonic oscillator -the dynamics investigated here does not admit such a simplified explanation. Still, it is interesting to compare the qualitative features of the memory kernels Γ(t) that emerge in these two cases. To be concrete, we consider a field with model A dynamics in d = 1 and we choose a Gaussian interaction potential V q = exp q 2 R 2 /2 , so that the integral in Eq. (27) can be computed in closed form, yielding We note that this function, which is positive for t = 0, becomes negative upon increasing t and, for t → ∞, it approaches zero from below, provided that the dragging speed v does not vanish (hence τ v < ∞, see Eq. (31)).
The kernel Γ(t) in Eq. (D11) is plotted in Fig. 7b, where we compare it to the one of Ref. [12]. The latter, which encodes the interaction with the viscoelastic fluid, can be readily obtained by combining Eqs. (D2) and (D3). As we show in Fig. 7b, choices of the values of the parameters γ i and τ i exist such that this second Γ(t) also becomes negative for sufficiently large t. In particular, in Fig. 7b we used two timescales τ i , with i ∈ {1, 2}, as reported in Tab. 1 of Ref. [12] for Wi = 0.17. We thus conclude that, in both models, the memory kernel Γ(t) features anti -correlations at long times. This is reminiscent of the negative memory often found in the context of rheology of complex fluids [66][67][68][69][70], and suggests that the underdamped modes displayed by the particle are indeed due to the negative response of the surrounding non-equilibrium environment [12], independently of its actual physical origin (i.e., due to either correlations or viscoelasticity).