Stochastic thermodynamics of a probe in a fluctuating correlated field

We develop a framework for the stochastic thermodynamics of a probe coupled to a fluctuating medium with spatio-temporal correlations, described by a scalar field. For a Brownian particle dragged by a harmonic trap through a fluctuating Gaussian field, we show that near criticality (where the field displays long-range spatial correlations) the spatially-resolved average heat flux develops a dipolar structure, where heat is absorbed in front and dissipated behind the dragged particle. Moreover, a perturbative calculation reveals that the dissipated power displays three distinct dynamical regimes depending on the drag velocity.

In this work we investigate the stochastic thermodynamics of a system consisting of a mesoscopic, externallydriven probe coupled to a fluctuating medium.The latter is here represented by a scalar field obeying nonconserved or conserved dynamics [51], and the whole system (probe+field) is immersed in a homogeneous heat bath at a fixed temperature T that induces thermal fluctuations.We provide suitable definitions for the heat, work and entropy exchanges between the probe, the field and the thermal bath, which are consistent with the first and second laws of stochastic thermodynamics.This approach is simpler than addressing a fully microscopic model, while it goes beyond a description based on a generalized Langevin equation (GLE) [52][53][54][55][56][57], which may capture temporal, but not spatial correlations of the environment.Notably, our framework is particularly powerp-1 arXiv:2305.16235v3[cond-mat.stat-mech]30 Apr 2024 (a) Sketch of a probe particle dragged by a trap with velocity v through a correlated field (red surface).(b) Cartoon of the particle at position Y, the modes of the field ϕ (red waves), and the particle-field interaction (blue box), which can store elastic energy H int .Particle and field exchange heats dQy and dQ ϕ , respectively, with a bath at temperature T .An external agent exerts work dW on the particle.Particle and field may exchange energies dW int y and dW int ϕ , respectively, via their interaction.For clarity we do not indicate work exchanges between field modes.Heat fluxes are considered positive when they are supplied to the system.ful when dealing with a fluctuating medium that is close to a critical point, where one can replace a complex microscopic dynamics with the simplest model belonging to the same universality class [51,58].As a fruit of our approach, one can thereby investigate the interplay between the probe, the spatio-temporal correlated field, and the heat bath under the light of stochastic thermodynamics.We illustrate the theory for the minimal yet insightful example of a probe particle dragged by a harmonic trap through a fluctuating Gaussian field.This mimics the setup commonly used to study friction and viscosity in (active) microrheology experiments [59][60][61], and is prototypical in stochastic thermodynamics.
The model.-We study a system formed by a mesoscopic probe at position Y(t) ∈ R d in d spatial dimensions and a scalar field, whose value at x ∈ R d and time t is denoted by ϕ(x, t) ∈ R. We assume ϕ to represent the relevant slow degrees of freedom of the medium, after integrating out all other possible faster modes.For example, ϕ may represent the local magnetization in a spin system [22,51], the height of a fluctuating membrane over a reference plane [11][12][13][14], the relative concentration of two chemical species in a near-critical binary liquid mixture [62], etc., although in some of these cases the actual dynamics is more complex than the one considered here.
The energy of the system is described by the Hamiltonian where H ϕ denotes the energy of the field, U is an external potential acting only on the probe, and H int encodes the interaction between the probe and the field.The probe and the scalar field are assumed to obey the following coupled overdamped Langevin equations: where F ext (Y, t) in Eq. ( 2) accounts for non-conservative external forces acting on the probe, while in Eq. ( 3) a = 0 or 1 for locally non-conserved or conserved field dynamics [51,63].Moreover, γ y,ϕ are friction coefficients, while ν and η (a) are independent Gaussian white noises satisfying the fluctuationdissipation relations ⟨ν i (t)ν j (t ′ )⟩ = 2γ y T δ ij δ(t − t ′ ) and ⟨η (a) (x, t)η (a) (x ′ , t ′ )⟩ = 2γ ϕ T (−∇ 2 ) a δ d (x − x ′ )δ(t − t ′ ), where ⟨. ..⟩ denotes the average over the realizations of the noises (note that we set Boltzmann's constant k B = 1 throughout the paper).Accordingly, in the absence of external driving, i.e., F ext = 0 and time-independent U, the system reaches a state of thermal equilibrium.Figure 1(a) is a sketch of the model with a harmonic external potential U.
First law of stochastic thermodynamics.-To analyze the thermodynamic properties of this system, we generalize ideas from stochastic thermodynamics [1,2].For the fluxes of energy and entropy between particle, field and thermal bath sketched in Fig. 1(b), this yields the expressions below (derived in Ref. [64]).First, the only systematic energy input into the system (probe+field) is the work exerted by an external agent on the probe.During a time interval [t, t + dt], this work is given by dW where d denotes non-exact differentials, and • indicates the usage of Stratonovich calculus throughout [1].Requiring energy conservation for the probe leads to the first law dU = dW + dW int y + dQ y , according to which any change dU in its potential energy can either be due to dW , or caused by the elastic coupling to the field, or else dissipated into the bath in the form of heat − dQ y .From Eqs. (2) to (5) it directly follows that the latter is given by dQ which is identical to the Sekimoto expression [2] for a single particle in a heat bath (i.e., H int = 0).Likewise, requiring local energy conservation for the field leads to the first law for the field [64] (An analogous, spatially-resolved first law is given in Ref. [64].) Accordingly, all the work done on ϕ is either stored in the field configuration dH ϕ (its "bending"), or dissipated in the form of heat dQ ϕ .Here, is the work locally exerted on the field by the elastic coupling.In general, dW int , because H int itself can store energy.From Eqs. ( 3) and ( 7) we find that the stochastic heat dQ ϕ takes the form generalizing Sekimoto's expression [1].Importantly, dQ ϕ is a density (or field) of local heat dissipation.For a = 0, Eq. ( 9) formally resembles the heat exchange of a single particle.The case a = 1 of conserved dynamics involves the inverse Laplace operator (∇ 2 ) −1 , which is nonlocal in space, and has thus to be interpreted in terms of its Green's function [65].In confined geometries, the latter depends on the boundary conditions [66,67], but here we will only consider the system in the bulk.Second law of stochastic thermodynamics.-Next, we define the irreversibility [1,34,44] associated with an individual, joint (probe+field) trajectory {Y, ϕ} t f ti as Here P denotes the path probability of the trajectory, starting with the joint probability density , while P R denotes the path probability of the corresponding time-reversed process 1 .By construction, the second law ⟨S tot ⟩ ≥ 0 holds [68].A central result of this manuscript is that the irreversibility S tot defined in Eq. ( 10) equals the stochastic (total) entropy production [64] with the heat dissipation Q y and Q ϕ defined in Eqs. ( 6) and ( 9), respectively -indicating the thermodynamic consistency of our framework 2 .Here ∆S sh y,ϕ = 1 We assume that ϕ and Y are both even under time reversal.
In the time-reversed process, the time-dependent external driving protocols U (x, t) + Fext(t) are reversed in time.The backward trajectories are initialized with . 2 The definition of irreversibility given in Eq. ( 10) generalizes those of, e.g., Ref. [44] [see Eq. ( 7) therein] and Ref. [34] [see Eq. ( 16) therein] in the presence of a probe.Similarly, Eq. ( 11) is compatible with the expression given in Eq. ( 26) of Ref. [32] for a generic system strongly coupled to a heat bath, at the level of (total) spatially integrated quantities.
} is the change in the fluctuating Shannon entropy.The first and second terms on the r.h.s. of Eq. ( 11) are the fluxes of entropy to the thermal bath from the probe and from the field, respectively.
If the system reaches a steady state, then ⟨dH int ⟩, ⟨dH ϕ ⟩, and ⟨dU⟩ vanish.This simplifies the first laws for probe and field, and implies Together with the second law and Eq. ( 11), this implies (12) with Q := dQ/dt, Ẇ := dW/dt, and Ṡtot ≃ S tot /(t f − t i ).Thus, in such steady states, ⟨ Ṡtot ⟩ is proportional to the average power ⟨ Ẇ ⟩ injected into the system -which is partly dissipated into heat directly by the probe, ⟨Q y ⟩, and partly through the field, ⟨Q ϕ ⟩.
Example: Particle dragged through a Gaussian field.-Within this framework, we consider the typical setup [59] in which a probe particle is confined in a harmonic trap of stiffness κ moving at constant velocity v.
2 /2, and the dissipation rate follows from Eq. ( 4) as In the long-time limit, the system reaches a steady state in the comoving reference frame with velocity v.As a minimal model for a near-critical medium -and as the simplest approximation of various complex systems displaying long-range spatial correlations -we consider a Gaussian field with Hamiltonian [51] The correlation length ξ = r −1/2 ≥ 0 controls the spatial range of the field correlations at equilibrium, and diverges upon approaching the critical point r = 0. We model the particle-field interaction, as in Refs.[19][20][21][22][23][24], by where V (x) ≥ 0 is a rotationally invariant function, with unit spatial integral and rapidly vanishing as |x| increases beyond the "radius" R of the particle, while λ sets the coupling strength.For λ > 0 and at equilibrium, configurations in which the field is enhanced around the particle are energetically favored.We express the particle position in the comoving frame as Z ≡ Y − vt + vγ y /κ, so that Z = 0 is the mean steady-state position for λ = 0. Similarly, we introduce the field φ(z, t) ≡ ϕ(z + vt − vγ y /κ, t) in the comoving frame.

Heat dissipation field. -
The first thermodynamic quantity we analyze is the spatially-resolved heat dissipation Q φ in the comoving reference frame [64].In Fig. 2 we show numerical results for Qφ of a field in d = 2 with a = 0.For small values of ξ [see panel (a)], we observe that Qφ is essentially negligible (within numerical uncertainties) and displays no discernible spatial structure.In contrast, if ξ reaches the order of the particle size R [panel (b)], regions of average heat dissipation ( Qφ < 0) or absorption ( Qφ > 0) start developing.Close to criticality, with ξ ≫ R [panel (c)], a dissipation dipole forms, with a region of heat absorption in front of the particle, whose spatial extent is approximately given by ξ.Hence, surprisingly, in front of the particle the heat bath supplies net energy as if it was coupled to a cooler object.Note that the second law implies Qy + d d x Qφ (x) < 0, but it does not preclude local heat absorption, i.e., Qφ (x) > 0 for some x.To further elucidate the origin of this effect, below we analytically investigate the statistics of particle and field for various values of v = |v| and ξ, and the dissipated power.
Particle statistics and bending of the field.-We now assume λ to be small, and use it as a perturbative parameter [66,[69][70][71][72]. Expressing Eqs.(2) and (3) in terms of Z and φ, we first obtain where α p ≡ p 2a (p 2 + r)/γ ϕ , and φ p is the Fourier transform3 of φ(x, t).From Eq. ( 17) we can calculate the mean stationary profile ⟨φ p ⟩ st , which is comoving with the particle, and which we call "shadow" [64].The latter is shown in Fig. 3(a) for a field in d = 1 with a = 0: the field is strongly bent around the particle, while ⟨φ(z)⟩ st ∝ exp(−|z|/ℓ ± ) for z → ±∞, with . Far from criticality (ξ → 0, see inset), the shadow vanishes, rationalizing the corresponding vanishing of ⟨ Qφ (z)⟩ in Fig. 2(a).Via a perturbative approach, we can investigate the particle fluctuations analytically.The moment generating function g(q) ≡ ⟨exp(−iq • Z)⟩ of the particle position at the lowest nontrivial order in λ reads [64] where σ 2 (u) ≡ T 1 − e −κu/γy /κ, and Interestingly, Eq. ( 18) predicts a non-Gaussian statistics of the particle position.In addition, the variance ⟨Z l Z m ⟩ c ∝ δ lm changes anisotropically compared to the case v = 0, so that the position distribution is elongated in the direction parallel to v [64,73].Furthermore, i∇ q g| q=0 gives the average displacement which turns out to be directed along −v; hence, the field induces a shift of the probability density of the particle position, which lags behind the average stationary value in the absence of the field [see the inset of Fig. 3(b)].Such a lag is the footprint of an underlying additional source of dissipation, which we analyze next.
Power fluctuations.-From the moment generating function in Eq. ( 18), we can access the distribution of the dissipated power.Indeed, rewriting Eq. ( 13) in terms of Z gives Ẇ = γ y v 2 − κv • Z, and thus encoding all moments of the dissipated power.To study the impact of the field, we focus on the average power where we identify the dissipation rate in the absence of the field ⟨ Ẇ ⟩ 0 = γ y v 2 ≥ 0, while ⟨ Ẇ ⟩ λ = −κv • ⟨Z⟩ ≥ 0 encodes the additional dissipation due to the field.According to Eq. ( 12), these are further equal to the entropy production at λ = 0 and λ > 0, i.We present a thorough analysis of these regimes in Ref. [64], for the case a = 0, and summarize it here.First, by inserting into Eq.( 16) the formal solution of Eq. (17) for φ p (t), an effective equation for Z(t) can be obtained, which is non-Markovian and nonlinear.However, in the limit of small v, the field is sufficiently fast to equilibrate around the particle at any instant in time [66,71].Conversely, for large v the evolution of the field is so slow that the particle encounters an effectively static field configuration.Accordingly, the first two regimes can be quantitatively captured by adiabatically replacing the field φ p (t) in Eq. ( 16) with its mean (comoving) profile ⟨φ p ⟩ st , i.e., with the shadow shown in Fig. 3(a), resulting in an approximately Markovian evolution of Z(t).In contrast, at intermediate values of v, the timescales of the particle dy- 4 The first and the third among these regimes are consistent with the scaling of drag forces reported in Refs.[19][20][21], where particles moving with constant velocity (i.e., without positional fluctuations) were studied.
namics are comparable with the relaxation time τ ξ ≃ γ ϕ ξ 2 of the field, and the adiabatic approximation is no longer accurate: the particle dynamics within the crossover between the first two regimes is dominated by the memory effects caused by the mutual influence of the particle and the field [64].Finally, in the third regime, the shadow becomes negligible compared to the (critical) fluctuations of the field, and the particle effectively encounters a rough landscape resulting from them.Notably, as the field approaches criticality (ξ → ∞), the amplitude of its fluctuations diverges [63], and thus the last (non-Stokesian) regime extends to low v.
Conclusions.-We developed a thermodynamically consistent framework to study the energetic and entropic flows for a probe in a fluctuating medium with spatiotemporal correlations, modeled here by a scalar field immersed in a heat bath.We showed that the mutual influence of the probe and the correlated environment leads to unusual thermodynamic properties, even for the simple example of a particle dragged by a harmonic trap through a Gaussian field.We showed that, close to criticality or, more generally, when the correlation length ξ of the field exceeds the particle size, a dipolar structure develops in the local heat dissipation field -with systematic heat absorption in front of the particle, and whose extent is determined by ξ.We thus expect such heat dipoles to be possibly revealed also in long-range interacting systems of diverse nature [74].Furthermore, the additional power required to drag the particle in the presence of the field features three regimes with distinct scaling in the drag speed v, among which are two non-Stokesian regimes -a feature that cannot be captured by a linear GLE.Far from criticality (ξ → 0), the medium is only weakly correlated, and indeed both the heat dipole and the additional dissipation vanish.The framework developed here can be readily applied to study more complex scenarios, e.g., non-quadratic Hamiltonians, or extended to systems of multiple particles in a common correlated (active) environment.Moreover, the recent and intense experimental investigations of colloidal particles in correlated fluids -such as viscoelastic media [60,75,76], binary liquid mixtures [25][26][27][28][29][30][31], and living cells [77,78] -could provide access to the phenomena described here.

Supplemental Material for "Stochastic thermodynamics of a probe in a fluctuating correlated field"
Davide Venturelli, 1, 2, * Sarah A. M. Loos, 3, 4, * Benjamin Walter, 5, 1, * Édgar Roldán, 4 and Andrea Gambassi In this Supplemental Material we provide additional analytical derivations and numerical results which support the discussion in the main text, and justify some of the results stated therein.

I. STOCHASTIC ENERGETICS FOR THE PARTICLE AND FIELD
In order to address the thermodynamic properties of the system described in our manuscript -composed of a probe and a field in interaction -we need to introduce suitable definitions of the energy flows and entropy changes associated with their stochastic dynamics.To this end, we utilize ideas from the framework of stochastic thermodynamics [1], and generalize them to the present case.In addition to the equations of motion, given in Eqs. ( 2) and (3) in the main text, the only assumption we make is energy conservation for both the probe and field dynamics.We use Stratonovich calculus throughout.
First, we assume that the only systematic energy input into the system is due to the work done on the probe by the action of an external agent.This can either result from changing in time the potential U(Y, t) -thereby increasing the potential energy of the probe -or from the application of an external force F ext .In agreement with classical mechanics, the work dW within the time frame [t, t + dt] is simply given by The total work W along a certain stochastic trajectory {Y(t), φ(x, t)} t f ti in the configuration space of the field and probe, from time t i to t f , is obtained by integrating over the infinitesimal increments, i.e., W [{Y, φ} Next, we consider the total energy U + H φ + H int of the system and how it may change, in order to identify all possible types of energy (ex-)changes within the system.First, the probe can store potential energy U(Y, t).The latter can change due to a probe displacement, i.e., a change of Y, or due to the motion of the trap itself, giving rise to the differential Note that the last term on the r.h.s. of this equation can be expressed, via Eq.( 1), in terms of the work dW as dW − F ext • dY.This fact will be used further below in Sec.I A. Second, the field can store energy H φ in its configuration, characterized by its local energy density h φ (x), such that Since we assume that H φ is not explicitly time-dependent, the total change dH φ of the internal energy is accordingly given by with Here dh φ (x) denotes the local change of internal energy of the field.Notably, depending on the choice of H φ , dh φ (x) may contain contributions from the neighboring points φ(x = x).For example, for a Gaussian Hamiltonian as given in Eq. ( 14) in the main text, the Laplacian term induces energy exchanges between neighboring locations, as one realizes by discretizing space on a lattice.
Third, the interaction H int between field and probe can store elastic energy.Again, defining the density of interaction energy h int such that the differential of the latter reads hence Here, we have defined the scalar dW int y and dW int φ (x) (the latter being a field) as the infinitesimal changes of the interaction energy due to the fluctuations of Y or φ(x), respectively.We interpret these energy flows as the work done by Y and φ, respectively.

A. First law of thermodynamics
As a final step, we connect the various energy flows, thereby obtaining the balance equations and the appropriate expressions for the heat.We start with dU(Y, t) given in Eq. ( 2): using Eq. ( 1), one finds Substituting the equation of motion of the probe dynamics [given in Eq. ( 2) of the main text] into the last equation, we find In the last step we have identified dW int y , given by Eq. ( 8).Since we assume energy conservation for the probe dynamics, all the work dW int y + dW done on the system that is not used to increase the potential energy of the probe (by an amount dU) must be dissipated into the heat bath in the form of the heat − dQ y .Thus, it follows from Eq. ( 11) that the first law for the probe dynamics reads dU = dQ y + dW int y + dW, (12) with the heat flow dQ y defined by as given in Eq. ( 6) of the main text.We now turn to the energetics of the field.We start with dh φ (x) given in Eq. ( 5) and we insert the equation of motion of the field [given in Eq. ( 3) of the main text], finding Identifying the last contribution on the r.h.s. of this equation with dW int φ [see Eq. ( 9)], we find Again, energy conservation requires that all the work dW int φ that is locally done on the field is either stored in the field configuration in the form of a dh φ (x), or dissipated by the field in the form of heat − dQ φ (x, t).Accordingly, the (local) first law for the field, valid at any point in space, reads We note that the l.h.s.describes the total change of internal energy of the field at a given point x, which generally also contains nonlocal contributions, as mentioned after Eq. ( 5).These terms can be interpreted as the energy exchanges within the field, i.e., between φ at different points in space.The integral form of Eq. ( 16) is given in Eq. ( 7) in the main text.This first law allows us to identify from Eq. ( 15) the heat dissipation field of φ as which is reported in Eq. ( 9) in the main text.
The accumulated heat and work along a joint trajectory {Y, φ} t f ti from t i to t f are obtained by integrating over the infinitesimal increments defined above, i.e.,

e[{Y, φ} t
with e ∈ {Q y , W, dQ φ (x)}.Explicit expressions of Qφ (x) for a scalar Gaussian field coupled to a dragged particle [see Eqs. ( 14) and ( 15) in the main text] will be reported in a future work.

B. Long-time limit
For t → ∞, the system approaches a state in which the probability densities of the probe's degree of freedom Y and the field value φ in the comoving frame of reference [e.g., x → x − vt w.r.t. the lab frame for a uniformly dragged particle] are time-independent.In this nonequilibrium steady state, the thermodynamic laws simplify due to the presence of additional constraints.In particular, the potential energy of the probe and of the field, as well as the interaction Hamiltonian, are on average conserved, such that dU = 0, dh φ = 0, and dh int = 0. Thus, the first laws for the probe and field given in the main text and reported here in Eqs. ( 12) and ( 16), respectively, simplify to and Further, from Eq. ( 7), we find Then, by using Eqs.( 19) and (20), this equation implies that the total dissipation of the entire process is given by the work applied to the probe, i.e., indicating the physical consistency of our approach.

II. STOCHASTIC ENTROPY PRODUCTION
In accordance with the framework of stochastic thermodynamics for a finite number of degrees of freedom coupled to a heat bath, and with the framework of irreversibility for fields [2,3], we define the total entropy production along a trajectory {Y(t), φ(x, t)} t f t=ti of the system from time t i to t f as where P denotes the path probability of the forward path of the joint trajectory of field and probe, starting from the probability density field ρ Y,φ [Y(t i ), φ(x, t i )].P R denotes the probability of the corresponding backward path in the time-reversed process, initialized with ρ Y,φ [Y(t f ), φ(x, t f )].In the time-reversed process, the time is running in the reverse order from t f to t i , hence the time-dependent external driving protocols U(x, t) + F ext (t) are reversed in time.(In the example of a particle dragged by a harmonic trap discussed in the Main Text, this corresponds to reversing the dragging velocity from v to −v.) Due to the Gaussian statistics of the noise, the path weights are proportional to the exponential of the Onsager-Machlup action A [2, 4], i.e., Since the noise terms of the probe and field are assumed to be independent from each other, the total action of the joint process is the sum of the actions of probe and field, reading [2] All the stochastic integrals in this section [such as the one appearing in Eq. ( 25)] are intended in the Stratonovich sense, but we will omit the • sign to ease the notation.Accordingly, in Eq. ( 25) we are actually adopting the Freidlin-Wentzell discretization convention; note that the final result for the entropy production [see, c.f., Eq. ( 30)] is however insensitive to the choice of such convention -see, for instance, Section III in Ref. [5], or Section 5.2.2.2 in Ref. [6].
In Eq. ( 25) we defined with a = 0 for non-conserved dynamics and a = 1 for conserved dynamics.Note that for a = 1 the operator (−∇ 2 ) a is nonlocal, hence, in the last term on the r.h.s. of Eq. ( 25), we use d 26) is the deterministic dynamical operator of the field dynamics in Eq. ( 3) of the main text, i.e., −∇ • J d = γ φ φ − η (a) .The action of the time-reversed process is accordingly given by Inserting the path probabilities into Eq.( 23) yields where the first term on the r.h.s.can be identified as the change, along the considered time interval, of the fluctuating Shannon entropy associated with the probability density field ρ Y,φ of the joint process, i.e., considering both the probe and the field.Furthermore, the second term on the r.h.s. of Eq. ( 28) can be simplified to In the last step, we have inserted the equations of motion ( 2) and (3) of the main text, and we converted the integrals over time into integrals along the trajectory {Y, φ} t f ti .As a final step, we identify in Eq. ( 29) the heat flows as given in Eqs. ( 13) and (17), and therefore obtain consistently with the total thermodynamic entropy production, given in Eq. (11) in the main text.From Eq. ( 10) of the main text [or equivalently, Eq. ( 23)] it automatically follows that, on ensemble average, the total entropy production is constrained by the second law S tot ≥ 0 [7,8].The fluctuations of S tot obey, by construction, the integrated and detailed fluctuation theorems, e −Stot = 1 and p[S tot ]/p[−S tot ] = e −Stot [7,8].Along long trajectories with t f t i , this further implies that the total dissipation Q = dx Q φ (x) + Q y asymptotically also fulfills these relations, as ∆S sh y,φ -contrary to Q φ and Q y -does not grow with time.

III. COMPARISON WITH THE DISSIPATED POWER PREDICTED BY A GENERALIZED LANGEVIN EQUATION
An established approach [9][10][11][12][13][14] to study the thermodynamics of driven particles in complex environments is based on using generalized Langevin equations (GLEs) [15,16].Within this description, the effect of the (slowly relaxing) medium on the particle is described via a friction kernel and a colored noise.However, GLEs generally capture only temporal correlations of the medium, but not dynamically-varying spatial correlations.In this section, we consider the conditions under which a GLE may capture the different scaling regimes in the mean dissipation rate, which we describe in the main text (see, e.g., Fig. 3 of the main text).
In particular, consider an arbitrary linear (overdamped) GLE of the form with a certain friction kernel Γ (which we do not need to specify further here) and a zero-mean Gaussian colored noise µ.We consider a nonequilibrium steady-state of the system with v = 0 and are particularly interested in the mean position (relative to the center of the trap), from which we can access the mean power [see Eq. ( 13) in the main text].
By averaging the GLE (31) over the possible realizations of the noise, one finds with u(t) ≡ Y (t) .In the steady state, the average velocity u(t) of the particle does not depend on time t.Accordingly, we can further simplify the previous equation as where we have introduced the constant The solution of Eq. ( 33) is Together with Eq. ( 21) in the main text, this implies that, for any linear GLE, with arbitrary friction kernel, the correction to the average dissipated power is This expression, being quadratic in v, corresponds to the first regime among those shown in Fig. 3(b) in the main text, and to the usual Stokes friction.We conclude that the other regimes discussed therein can only possibly emerge in the presence of nonlinearities in the effective evolution equation of the probe.Now, let us consider a generalization of Eq. ( 31) to the case of a nonlinear friction kernel [17,18], of the generic form: with some arbitrary nonlinear function f [19].Repeating the same steps as above, we find that the steady-state solution is again of the form given in Eq. ( 35), with Clearly, if we allow f to be nonlinear in Ẏ (i.e., "non-Stokesian friction"), this equation can also reproduce the different regimes we found in the friction (equivalently, mean power) described in the main text.For example, the scaling in the second regime in Fig. 3(b) in the main text formally corresponds to a situation where Eq. ( 37) has a contribution to the friction of the form . This is a rather unusual friction term that seems difficult to justify on heuristic or phenomenological grounds alone.In contrast, by systematically integrating out the field degrees of freedom φ(x, t) from the joint field-particle dynamics [given in Eqs. ( 2) and (3) in the main text], we can obtain an effective description for the motion of the particle in the form of a nonlinear GLE akin to Eq. (37) [see, c.f., Eq. ( 91)] -as discussed in the next sections.

IV. BROWNIAN PARTICLE AND GAUSSIAN FIELD: INDEPENDENT PROCESSES
We summarize here the well-known solutions of the non-interacting processes for Z(t) and ϕ p (t) [equivalently, Y(t) and φ p (t)] which are obtained by setting the coupling constant λ = 0 in Eqs. ( 16) and ( 17) of the main text.This case forms the basis for the perturbative calculation of the moment generating function of the probe particle position, which we derive in Section V.

A. Brownian particle in a harmonic potential
The equation of motion for a Brownian particle in a harmonic potential with is that of the Ornstein-Uhlenbeck (OU) process, which we briefly revise here for clarity and to set the notation for the following sections.Each component Z j (with j ∈ {1, . . ., d}) of the particle position Z is a Gaussian and Markovian process.Accordingly, the propagator of the process is where is the expectation value of the particle position Z(t) given that the process started at Z 0 at time t 0 , and is the corresponding variance.Because of isotropy, σ 2 0 (t, t 0 ) is independent of the component j ∈ {1, . . ., d}.For the process we are currently interested in, they read Note that the time-translational invariance of the system implies P (Z, t|Z 0 , t 0 ) = P (Z, t−t 0 |Z 0 , 0).Using the Langevin equation (39) or, alternatively, Eq. ( 41), one can readily calculate the connected two-time correlation function which, again, is independent of the value of j, while correlations Z j (t 1 )Z i =j (t 2 ) c vanish.At long times, the system attains the thermal equilibrium described by the probability distribution

B. Dynamics of the free field
We now consider the Langevin equation for the free field, i.e., the one obtained from Eq. ( 3) in the main text by setting H int = 0, and using the Gaussian H φ given in Eq. ( 14) In Fourier space, this equation reads where, as in the main text, we conveniently introduced α p = p 2a (p 2 + r)/γ φ , with a = 0 or 1 corresponding to non-conserved or conserved dynamics, respectively.The zero-average Gaussian noise η p above is characterized by the correlations where Ω φ (p) = 2T p 2a /γ φ .We adopted the Fourier convention f 48) is the same as the one describing the position of a Ornstein-Uhlenbeck particle, whence In the previous expressions, we introduced the superscript (0) in order to distinguish these expressions from those derived in, c.f., Section V, which pertain to reference frame moving at velocity v.Note that, by construction, p (s 2 , s 1 ).In the limit in which the time t 0 of initialization of the field is sent to t 0 → −∞, one recovers the equilibrium correlator which is a function of the time difference τ = s 2 − s 1 only.It is also customary to define the response function and linear susceptibility of the free field [20] respectively, where Θ(t) is the Heaviside theta function.The functions χ (0) p (t) and G (0) p (t) are time-translational invariant (also in a generic, possibly nonequilibrium, stationary state, since the equation of motion is).Moreover, in equilibrium, they are linked to the correlator in Eq. ( 52) by the fluctuation-dissipation theorem [20] T χ (0)  p (τ ) = −Θ(τ )

V. DERIVATION OF THE MOMENT GENERATING FUNCTION OF THE PARTICLE POSITION
Here we provide the derivation of the cumulant generating function of the particle position reported in Eq. ( 18) of the main text.We start from the coupled equations of motion for the variables Z = Y − vt + vγ y /κ and ϕ(z, t) ≡ φ(z + vt − vγ y /κ, t), i.e., Note that, for λ = 0, the second equation is formally the same as Eq. ( 48) with α p → α p − ip • v, which describes the free Gaussian field.Its solution in the stationary state thus remains the same upon replacing the equilibrium correlator C (0) p (t) and the free susceptibility χ (0) p (t) in Eqs. ( 52) and ( 53) by Starting from the Langevin equations ( 56) and ( 57), and assuming the coupling constant λ to be small, it is possible to derive a (forward) evolution equation which describes the marginal probability distribution P 1 (z, t) of the particle coordinate alone.The derivation is presented in Ref. [21], where a more general case is treated, in which a second particle also interacts with the field via an interaction potential V (z) p ; the present case corresponds to V (z) p → 0. We report here only the final result: Here is the Fokker-Planck operator for an Ornstein-Uhlenbeck particle [22], while the second term on the r.h.s. of Eq. ( 61) involves a convolution of the two-time probability distribution P 2 (z, t; x, s) with a memory kernel L(r, u).This typically emerges in non-Markovian processes, where one usually obtains a hierarchy of master equations relating the n-point distribution P n (x n , t n ; x n−1 , t n−1 ; . . .; x 1 , t 1 ) with P n+1 (see for instance Refs.[23][24][25][26]).This kernel reads where the notation ∇ j z is shorthand for ∂/∂z j , and the summation over the repeated indices j and l is implied.This last expression had already been reported in Ref. [21] but, unfortunately, with a confusing notation.
We want to find an expression for the characteristic function We start by taking the Fourier transform of the various terms in the master equation (61), which leads to where, from Eq. ( 62), and we have introduced The solution of Eq. ( 65) is formally given, at the lowest nontrivial order in λ, by where we introduced i.e., the moment generating function of the uncoupled (λ = 0) particle, and where we indicated by g 1|1 the propagator of the Ornstein-Uhlenbeck process in momentum space.The latter is easily found by solving the homogeneous part of Eq. ( 65) via the method of characteristics [27], or even by simply Fourier transforming the real space propagator given in Eq. ( 41): both ways lead to where σ 0 (t − t ) was introduced in Eq. ( 45).
Let us now focus on the stationary state, for which Equation ( 68) then becomes where in the second line we used Eq. ( 70).Similarly, manipulating Eqs. ( 63) and ( 67) renders in terms of the double Fourier transform P 2 of the two-point probability distribution P 2 , with respect to both its spatial arguments.Since we are neglecting terms of order higher than λ 2 , then P 2 above can be replaced by its expression for λ = 0, i.e., where P 1|1 is the Ornstein-Uhlenbeck propagator given in Eq. ( 41), with Fourier transform g 1|1 provided in Eq. ( 70).In the long-time limit, we can replace P 1 (x, s) by the stationary distribution P st 1 (x) eventually attained by the uncoupled particle in its harmonic trap, see Eq. ( 47).We recall that the Fourier transform of P st 1 (x) is the function g (0) (p) given in Eq. ( 69).This yields Equation ( 73) can be further simplified by noting that, for any u-independent value A [and, in particular, for A ≡ p • (p + k)], we can rewrite [28] where we recall that χ  58) to ( 60)], and where we introduced with σ 0 given in Eq. (45).Equation ( 76) readily follows as a direct consequence of the fluctuation-dissipation theorem for the free field, recalled in Eq. ( 55).From Eq. ( 73), after integrating by parts in u, we thus get where G p (u) is the field response propagator in the comoving frame, see Eq. (60).We finally substitute the expression of f (k) given in Eq. ( 78) into Eq.( 72) for g(k), and we simplify the integral over u by recognizing that This gives the moment generating function Upon taking the logarithm of Eq. ( 80) and expanding for small λ, we obtain the cumulant generating function reported in Eq. ( 18) of the main text.In particular, computing −∇ 2 q g(q)| q=0 gives the correlation function We note that the contribution ∝ λ 2 actually vanishes for l = m, as it follows from the fact that, otherwise, the integrand would be an odd function of some component of the vector p.It thus appears that no cross-correlations arise between the various orthogonal components of the displacement, i.e., Z l Z m c ∝ δ lm , but the variance is anisotropically modified along the directions parallel and perpendicular to the trap velocity v (see also Fig. 1).In the critical case, as well as for small velocities or small thermal fluctuations, it is possible to prove analytically that the correction in Eq. ( 81) is positive.Note that in the absence of the field (i.e., for λ = 0) the variance does not change, at finite v, compared to the case v = 0.

VI. AVERAGE PARTICLE POSITION AND DISSIPATION RATE WITH NON-CONSERVED DYNAMICS
In this Appendix, we analyze the behavior of the average particle position Z reported in Eq. ( 19) in the main text, together with its consequences on the dissipation rate discussed therein, focusing on the case of non-conserved dynamics.

A. Typical time and length scales
First, it is useful to introduce the following timescales: where z = 2 + 2a is the dynamical critical exponent of the field (z = 2 for non-conserved dynamics) and ξ the correlation length of its fluctuations.We recall that R is the length scale which characterizes the interaction potential and thus sets the particle "radius" [see Eq. ( 15) in the main text], κ is the stiffness of the harmonic trap, while the parameters γ y,φ are the friction coefficients of the particle and the field, respectively [see Eqs. ( 2) and (3) in the main text].
The first timescale τ κ is simply the relaxation time of the particle in its harmonic trap.The second is the time τ v taken by the moving trap to cover a distance corresponding to the particle radius; equivalently, τ −1 v estimates the shear rate near the driven particle [29].The next two scales, instead, characterize the field φ.In particular, τ R is the time that a critical field φ takes to relax across a distance ∼ R; similarly, τ ξ is the relaxation time of the field φ across its correlation length ξ, and it represents its slowest timescale in the vicinity of the critical point (for non-conserved dynamics).
Upon appropriately rescaling the integration variables, both the cumulant generating function in Eq. ( 18) of the main text and the moments of the particle position Z can be recast in the form of scaling functions which depend only on dimensionless ratios of the typical time and length scales of the system.Consider, for instance, the average position in Eq. ( 19) of the main text in the case of non-conserved dynamics, where the integral over momenta p can be computed upon choosing a Gaussian interaction potential V p = exp −p 2 R 2 /2 .One possible choice is to rescale momenta by ξ −1 , so that after some manipulation one gets where the second factor highlighted on the r.h.s. is dimensionless.Similarly, we introduced above the auxiliary function which depends on dimensionless parameters, while the parameter can be interpreted as the squared Péclet number of the system, i.e., χ = (Pe) 2 .Alternatively, we note that in the context of microrheology experiments conducted in viscoelastic media [29] one can identify the Weissenberg number Wi ≡ τ s /(2τ v ), where τ s is the (single) relaxation timescale of the medium.It follows that the parameter χ in Eq. ( 86) can be seen as the product of two Weissenberg numbers, one for each of the two dominant relaxation timescales of the field, i.e., τ R and τ ξ .The remaining dimensionless parameters {π i } encode the ratios of the typical length or time scales of the system.In the limit in which τ ξ τ κ (strong confinement [30,31], or close to criticality), the function f (a, b, {π i }) in Eq. (85) simplifies, and the dependence on the three parameters {π i } can be replaced by that on a single parameter π 4 defined as Note that T may be regarded as the effective particle radius, once its average thermal fluctuations have been taken into account: indeed, T /κ is the typical thermal fluctuation of the particle position in a trap at thermal equilibrium.We thus obtain in this limit with the function g defined as g in Eq. ( 84), but with f replaced by f above.

B. Behavior as a function of the drag speed v
The average position Z in the steady state is plotted in the inset of Fig. 3(b) in the main text as a function of the trap velocity v, for the case of non-conserved dynamics in d = 1.For any finite value of r, i.e., of the correlation length ξ = r −1/2 , one can generically identify three regimes and a crossover: • Low-velocity (adiabatic) regime, for 0 < v v I .When v = 0 the system is in equilibrium, and the average particle position is not modified by the presence of the field [28], thus yielding Z(v = 0) = 0. Inspecting Eq. ( 19) of the main text shows that Z ∝ −v, i.e., a linear dependence of Z on v for small v.
Here | Z | reaches a maximum and behaves in a possibly nonmonotonic way (depending on the choice of the various parameters).
• High-velocity (adiabatic) regime, for v II v v III .The amplitude of the average position starts decaying as • "Depinning" regime, for v v III .At very high speeds, one can check that Z tends to a small but finite value.These regimes correspond to those observed for the correction ∝ λ 2 to the dissipation rate Ẇ λ ∝ v Z in Fig. 3(b) in the main text, which grows like ∝ v 2 for small v (as it is the case for the usual Stokes friction), while it is v-independent for intermediate v, and finally grows like ∝ v in the last regime.
We can understand these behaviors starting from the effective equation of motion for the particle coordinate Z(t), which can be obtained from Eqs. ( 56) and ( 57) by solving for ϕ p (t) in the second (linear) equation, and by inserting this result into the first equation [21,30,32].This yields where the field susceptibility χ p (u) is given in Eq. ( 59), while we introduced the colored Gaussian noise which has zero mean and the correlator C p (t) [see Eq. ( 58)].Equation ( 91) is markedly non-Markovian because of the presence of a (nonlinear) memory term [i.e., the term that is proportional to χ p (u)], and a colored multiplicative noise term ζ p (t).In the limits of very small and very large drag speed v, however, we physically expect to recover an approximately Markovian description.First, for very small v, the field is fast enough to equilibrate with the particle being fixed at its actual position, so that an adiabatic approximation is applicable [21,28,33,34].Second, for very large v, the evolution of the field is so slow that the particle effectively encounters a field landscape that is effectively static, i.e., quenched.
To make these statements more substantial, it is interesting to investigate the properties of the nonequilibrium steady state of the field, characterized by the formation of a comoving stationary profile around the particle, which we dubbed shadow in the main text.Its analytical expression is found by averaging Eq. ( 17) of the main text over thermal fluctuations and by setting ∂ t ϕ p = 0, which yields, at the leading order in λ [31], where we used Eq. ( 18) of the main text, for λ = 0, in order to express the average exp(−ip • Z) = g(p) at the numerator of the resulting expression.The shadow is plotted in Fig. 3(a) of the main text for the case of d = 1 and non-conserved dynamics.Qualitatively, as the value of the speed v increases, the shadow becomes increasingly asymmetric (with its wavefront becoming steeper and its wake longer), and its overall amplitude decreases towards zero.In the limit of a point-like particle -which is obtained by choosing an interaction potential V (x) = δ(x)one can explicitly compute the inverse Fourier transform of Eq. ( 93) in d = 1 non-conserved dynamics by complex integration, finding Accordingly, the shadow decays exponentially, in front or behind the particle (i.e., for z → ±∞), over two generically distinct length scales respectively, as reported in the main text.From this expression of ± it appears that 1/(γ φ ξ) is a natural velocity scale of this problem and, in fact, we shall see below that it coincides with the velocity v I under which the dynamics is adiabatic.Close to equilibrium, i.e., for v v I = 1/(γ φ ξ), Eq. ( 94) reduces to which is symmetric as expected, and extends over distances of O(ξ).Conversely, for large speeds v v I the amplitude of the shadow decreases as 1/v, and decays in its front over a typical length 1/(vγ φ ) ξ, and in its wake over the length ξ 2 vγ φ ξ.Notably, by using a more realistic interaction potential V (x) that takes into account the finite particle radius R ξ, Eq. ( 94) would still describe the tails of the shadow, i.e., its behavior for |z| R. As we emphasized in the main text, choosing positive values of the coupling constant λ and V (x) in Eq. ( 15) results in an effective attraction between the particle and the shadow; we are interested here in how this attraction modifies the steady-state average particle position Z .Let us assume initially that this deterministic attraction dominates over the (non-Gaussian) thermal fluctuations induced by the field, which are encoded in the noise term ζ p (t) in Eq. ( 92).Focusing on the two Markovian limits described above, we can approximately replace the field ϕ p in Eq. ( 56) by its mean value, i.e., by the shadow given in Eq. ( 93).This approximation implies that The second term on the right hand side represents a field-induced effective force which, in the steady state defined by ∂ t Z(t) = 0, counterbalances the restoring force −κ Z of the harmonic trap.For small speeds v, the shadow is essentially symmetric [see Eq. (96)] around its center which we will denote as z p .Linearizing Eq. (97) around the unperturbed position Z = 0 (for simplicity, we consider below the case d = 1) renders in general where the expression of the effective linear strength κ φ and of z p introduced above can be found on the basis of the expression of the shadow in Eq. ( 93) -as they are not necessary for our argument, we do not provide them here.In the steady state and at leading order in the coupling λ, we thus find At small but nonzero v, one expects the center of the shadow to lag behind the average particle position (which effectively "drags" the shadow), i.e., z p = av + O v 2 with some a < 0. Accordingly, Eq. ( 99) with this expression for z p predicts that | Z | has a linear dependence on v for small velocities, as observed in Fig. 3(b) in the main text (see the inset), and as discussed above.
The simplified description presented above is expected to become less accurate as v increases, because the rearrangement of the field can no longer be assumed to be instantaneous, and the particle actually moves considerably (due to the dragging) while the shadow forms.The spatial extension of the shadow is of O(ξ) in this regime [see Eq. ( 96)], so it builds over a timescale τ ξ γ φ ξ 2 in the case of non-conserved dynamics [see Eq. ( 82)].This timescale should be compared with the time τ taken by the trap to span a distance of the order of the field correlation length ξ, i.e., τ ξ/v.The two timescales balance at which we identify as the end of the first regime, i.e., the beginning of the crossover in Fig. 3(b) of the main text.This estimate assumes also that τ ξ τ κ , i.e., that the field relaxes faster than the typical timescale set by the harmonic trap (which is typically the case in experimental realizations with colloidal particles -see Section VII B).
For intermediate velocities v I ≤ v ≤ v II , the behavior of Z is fully determined by non-Markovian effects which, contrary to the previous regimes, cannot be explained in terms of a simplified model.Conversely, for v > v II , we can again replace the field by its expectation value as in Eq. (97).In this high-velocity regime, the shadow is so asymmetric [see Fig. 3(a) of the main text] that it is no longer meaningful to treat it as a point-particle concentrated in its center z p [as we essentially did in Eq. ( 98)].In spite of this difficulty, it actually turns out that its form can still be determined analytically.In particular, from Eq. ( 93) one can first compute where the approximation introduced in the last equality is expected to be valid when v exceeds all the relevant ratios of length and time scales of the system.This requires, inter alia, that [see Eq. ( 82)], which justifies the replacement of the field ϕ p by its mean value in Eq. ( 97).Moreover, we assume (as done above and in the main text) that the interaction potential V (x) is characterized by a single length scale R (i.e., the particle radius), so that V p provides an effective cutoff over momenta p 1/R.In the denominator of Eq. ( 101), we can thus impose the large velocity condition α p pv in correspondence of p ∼ 1/R, so that this condition is satisfied for all values of p ≤ 1/R: this gives [with reference to Eq. ( 82)] Accordingly, the second threshold velocity v II is expected to be provided by the smallest velocity v for which both conditions in Eq. ( 102) and ( 103) are satisfied, which depends on the particular choice of parameters and specific details of the system.In Eq. ( 101) we also recognize the term V p ≡ V p exp −T p 2 /(2κ) , featuring the Fourier transform of the interaction potential V (x), corrected by a factor which takes into account the mean squared displacement of the particle in the harmonic trap induced by the thermal fluctuations [see also T in Eq. ( 88)].Accordingly, in this limit, the shadow becomes [Here we assumed that ϕ(x) st → 0 for |x| → +∞, as expected from Fig. 3(a) in the main text.]As a consequence, the effective force term in Eq. (97) scales in this regime as v −1 , which produces an analogous dependence on v for the average particle position Z in the steady state.Due to this very v −1 dependence, we expect this deterministic effect to die out for very large v.In this last ("depinning") regime, the field is so slow with respect to the particle motion that no shadow can build up, and the particle effectively sees a rough landscape whose features are solely determined by the thermal fluctuations of the field.In order to describe this situation, we have to go back to the particle effective equation (91), the average over thermal fluctuations of which yields By using Novikov's theorem [35,36] where F [ζ] is any functional of the (Gaussian) noise ζ, we can compute the expectation value where in the last step we used the equation of motion (91) of Z(t), we indicated by . . .0 the expectation value over the unperturbed process with λ = 0, and we discarded higher-order terms in λ.The leading-order expression for the dynamical structure factor was computed in Ref. [28], with σ(u) given in Eq. (77).Setting ∂ t Z = 0 in Eq. (105) and taking the limit for t → +∞ on its r.h.s. then yields the steady-state average position This last expression indeed coincides with Eq. ( 19) in the main text, for any finite temperature T > 0, which can be seen upon integrating by parts in du and by using the fluctuation-dissipation theorem in the form of Eq. ( 76).
It is simple to verify that the term proportional to χ p (u) in Eq. ( 109) decays as 1/v for large v, as it encodes the deterministic effect of the shadow.At large v, the only term that contributes is the one originating from the noise ζ p (t) and proportional to the field correlator C p (u), which explicitly becomes This term is proportional to the strength T of the thermal fluctuations and the coupling λ 2 to the field, in agreement with the interpretation we provided earlier in this section -i.e., that the features of the field landscape encountered by the particle in this regime are determined by the critical fluctuations of the former.The identification of the value v III of the velocity v at which this regime begins is, however, not straightforward: it marks the value at which the field-induced thermal fluctuations start to dominate over its deterministic attraction.Since these fluctuations exhibit a rather strong dependence on the distance r from the critical point [as it appears by inspecting the field correlator C p (u) in Eq. ( 58)], then v III also shows a similar dependence on the correlation length ξ = r −1/2 (as well as on the temperature T ).In particular, Fig. 3(b) in the main text shows that v III decreases upon approaching the critical point r = 0, while the extents of the first and second velocity regimes decrease.

C. The case of critical non-conserved dynamics
It is instructive, at this point, to inspect the case of critical non-conserved dynamics, i.e., to consider r = 0 and a = 0.In d = 1 and with a Gaussian interaction potential V p (with zero mean and variance R), the average particle position given in Eq. ( 19) of the main text can be simplified as The corresponding curve for the dissipation rate is plotted in Fig. 3(b) in the main text (solid green line with ξ = ∞).
As the critical point is approached, we note that the crossover velocity v I described in Section VI B decreases; exactly at criticality, the entire low-velocity (adiabatic) regime disappears, in agreement with the expression of v I given in Eq. (100).Furthermore, for r → 0, critical fluctuations are found to play a major role, as signaled by the fact that the crossover velocity v III also decreases considerably.Interestingly, we note that the expression of Z for r = 0 in Eq. ( 112) has a finite, non-vanishing limit for v → 0, while Eq. ( 19) of the main text, in the same limit but with a generic value of r, correctly vanishes independently of r = 0; in fact, at equilibrium, the presence of the field does not modify the equilibrium distribution of the particle [28].We are thus led to the conclusion that the limits r → 0 and v → 0 do not commute: the physical interpretation is that, when r → 0, the relaxation timescale τ ξ of the field diverges [see Eq. (82)] and therefore the system is out of equilibrium for any finite value of v, no matter how small.

VII. NUMERICAL SIMULATION OF THE STOCHASTIC DYNAMICS
In this Appendix we provide further details regarding the numerical simulation of the coupled stochastic dynamics of field and particle.

A. Numerical integration scheme
The joint stochastic evolution of the field and the particle, given by Eqs. ( 2) and (3) in the main text, is integrated in time by using the stochastic Runge-Kutta method [37] with a fixed discretization time interval ∆t.While the particle evolves in continuous space (i.e., off-lattice), the field is defined on Λ = ∆x d i=1 i e i , 1 ≤ i ≤ L , a cubic grid of L d lattice points spaced ∆x apart.We set out with φ(x, t) and Y(t) and define the first auxiliary increment Here, the field gradient is approximated by its lattice discretization, Further, ν(t) denotes an uncorrelated Gaussian random vector with unit variance.The first increment K 1 is added to Y(t) to obtain the auxiliary value Ŷ(t) = Y(t) + K 1 (t).The stochastic increment of the field is then evaluated by using the auxiliary value as where the Laplacian is approximated via Hereafter, the position of the particle is updated using the new field configuration by defining the second auxiliary increment The noise ν(t) used in this second increment is identical to the one previously used in Eq. ( 113).Finally, Y(t) is updated as This scheme can be suitably adapted in order to to include, e.g., conserved field dynamics.As a second-order stochastic Runge-Kutta integration scheme, the local truncation error in time is of order O(∆t 3 ), and can be improved by using higher-order schemes [37].The initial condition φ(x, t = 0) ≡ φ(x) of the field φ(x, t) is drawn from the Gaussian equilibrium distribution of the decoupled field (i.e., with λ = 0, see Section IV B).By applying a d-dimensional Fourier transform to φ, one obtains where p sums over the reciprocal lattice Λ We randomly draw the complex modes φ p ∈ C for momenta belonging to the positive half-cube p ∈ p ∈ Λ −1 |p 1 ≥ 0 , containing all vectors in the reciprocal lattice with positive first component; while the magnitude is drawn from a Gaussian random distribution N , i.e., φ p ∼ N 0, 2π (L∆x) 2 p 2 + r , the complex phase θ ∈ U([0, 2π)) is drawn from a uniform random distribution U, leading to a random initial mode φ p = φ p e iθ .If r = 0, we set φ p=0 = 0.Then, to ensure that φ(x) is real, we set φ −p = φ * p for the remaining vectors in the negative half-cube.The initial condition then follows from inserting the random modes into Eq.(119).The particle is initialised at the centre of the trap at t = 0. Since neither the particle nor the field start in their joint equilibrium distribution (λ = 0), we allow for thermalization by evolving the system for typically 10 4 timesteps at a time discretization of ∆t = 10 −2 , after which we assume the system to have settled into a (nonequilibirium) steady state.

B. Choice of parameters
Here we discuss the choice of the values of the parameters of the model to be used in the numerical simulations, so that they eventually correspond to the typical time and length scales observed in experiments with colloidal particles.As a prototypical example, we consider the case of a µm-sized colloidal particle immersed in a binary liquid mixture close to its demixing transition [38][39][40][41][42], so that the field φ represents the relative concentration of the two species that compose the mixture.It must be emphasized, however, that our model is not meant to reproduce quantitatively the results of such experiments (for which hydrodynamic effects would need to be taken into account), but rather to highlight the role played by spatial correlations.Moreover, other very diverse physical systems (such as inclusions in lipid membranes [43][44][45][46], microemulsions [47][48][49][50], or defects in ferromagnetic systems [51][52][53][54][55][56]) also fall within the scopes of our model, and they may entail very different time and length scales.
Binary liquid mixtures can undergo a demixing phase transition close to room temperatures [57], in correspondence of which the correlation length ξ = ξ(T ) diverges.In real experiments, achieving large values of ξ is generally challenging, as it requires a fine-tuning and highly accurate control of the temperature throughout the sample.However, in the following we will assume that ξ can indeed be made much larger than the particle radius R, in order to emphasize the qualitative effects of spatially extended correlations, which are the central topic of this paper.
Note that our model features several parameters, but only a limited number of physical units (i.e., mass, length, time and temperature).As detailed in Ref. [30], close to the critical point, the system can be conveniently described in terms of a reduced set of dimensionless parameters, which correspond to ratios of the typical time and length scales of the system described in Section VI A. These parameters are The first is the Weissenberg number w = Wi, which compares the shear rate due to dragging with the typical relaxation time of the field τ R , over length scales of the order of the particle radius R. The second parameter, ρ, compares the relaxation time of the field with the time scale τ κ set by the harmonic trap.The effective coupling g quantifies the strength of the interaction between the particle and the field.Finally, is the ratio between the thermal length l = k B T /(2κ) (i.e., the mean-squared displacement of the particle in the trap according to equipartition in equilibrium) and the particle radius R. Specializing these expressions to the case of critical non-conserved dynamics in d = 2 (see Section VI A and Ref. [30]) yields g = λ 2 /(κR 2 ), w = Rvγ φ /2, and ρ = κR 2 γ φ /γ y .
Having identified the relevant parameters in Eq. (120), we now discuss which values of them are within experimental reach.The typical magnitude of the thermal length l can be estimated at room temperature T 300 K by assuming a typical value of the optical trap strength κ 0.5 pN/µm [40], yielding l 100 nm, and hence 0.1 for a µm-sized particle.Next, the typical relaxation time scale τ φ of a near-critical binary mixture can be estimated by using model H [58]. Within mode-coupling theory and for momenta p such that pξ 1 (i.e., for small ξ, which has been the case in past experiments [38][39][40][41][42]), the relaxation timescale of the field is given by [59] τ −1 φ (p) D ξ p 2 , ( where and η is the fluid viscosity.By comparison with the diffusion coefficient D R of the particle [42] D R = k B T 6πηR 0.22(µm) 2 s −1 , (123) one can estimate the relaxation time scale of the field over distances of the order of the particle radius as where we chose ξ R/50 (which is common for water-lutidine mixtures [41,42]).Since the typical relaxation time τ κ of optical traps is of the order of a few tens of ms [41], this gives ρ 10 −1 .This value is expected to increase if the correlation length ξ can be made larger, since this generally entails longer relaxation times τ φ [20].Next, with typical drag speeds up to a few tens of µm/s [60], one can explore small Weissenberg numbers up to w 10 −1 .On the other hand, the effective coupling g quantifies the strength of the interaction between the particle and the field, FIG. 1. Upper panel: scatter plot of the particle position, in the comoving frame, as measured in numerical simulations in d = 2 and for a field with non-conserved dynamics.The average positions for v = 0 and v = 5 are indicated with a red dot and a blue square, respectively.Lower panels: histograms of the particle position along the directions parallel (z ) and perpendicular (z ⊥ ) to the trap velocity v.In the simulations we used the parameter set given in Eq. ( 125), with T given in Eq. (88).In the plot, we indicated with a solid black line the Gaussian distribution, which is recovered for v = 0 (boundary of the gray regions) or in the absence of particle-field interactions, i.e., for λ = 0 (leftmost curve in the left panel).
and its amplitude thus depends on the specific coupling mechanism realized in a certain experiment.At present, it is still unclear how to estimate its order of magnitude on the basis of past experimental data, but we clearly expect any overall qualitative effect to be enhanced if g can be made larger.We finally keep an eye for consistency on the ratio Θ ≡ τ v /τ κ = ρ/(2w) [see Eq.
where space is expressed in units of the lattice spacing ∆x, and we furthermore chose a discretized time step ∆t = 10 −2 .This choice corresponds, in fact, to w = 0.25, ρ = 0.2, = 0.66, Θ = 0.4, and g = 125.The lattice size L = 128 is chosen sufficiently large so as to accommodate the tails of the field shadow [see Fig. 3(a) in the main text], and to avoid spurious field currents due to the stirring that the particle would generate if it were dragged around the periodic boundaries [30].

C. Numerical measurement of the distribution of the particle position
Using the numerical simulation described above, we determined the statistics of the particle position, shown in Fig.1:(a) Sketch of a probe particle dragged by a trap with velocity v through a correlated field (red surface).(b) Cartoon of the particle at position Y, the modes of the field ϕ (red waves), and the particle-field interaction (blue box), which can store elastic energy H int .Particle and field exchange heats dQy and dQ ϕ , respectively, with a bath at temperature T .An external agent exerts work dW on the particle.Particle and field may exchange energies dW int

Fig. 2 :
Fig. 2: Average local heat dissipation rate ⟨ Qφ(z)⟩ by the field into the bath, in a frame comoving at velocity v with the harmonic trap (see Fig. 1).These results of numerical simulations [64] refer to a Gaussian field in d = 2 with nonconserved dynamics and, from (a) to (c), increasing values of ξ/R, where ξ is the correlation length and R the probe size.The dash-dotted circles indicate the size of the probe.Panel (d) shows ⟨ Qφ(z)⟩ along the drag direction z ∥ , for z ⊥ = 0 and the values of ξ considered in (a-c).We used λ = 5, v = 5, lattice size L = 128, time step ∆t = 10 −2 , a Gaussian potential V with variance R = 4, and we set all other parameters to unity.

Fig. 3 :
Fig. 3: (a) Steady-state expectation value ⟨φ(z)⟩st of the field in the comoving frame, for various drag velocities v > 0, in d = 1 and with non-conserved dynamics.The shadow ⟨φ(z)⟩st flattens upon increasing v (at fixed correlation length ξ = 10, main plot) or upon decreasing ξ (at fixed v = 5, inset).⟨φ(z)⟩st decays exponentially upon increasing |z| with different decay lengths I. Stochastic energetics for the particle and field 1 A. First law of thermodynamics 3 B. Long-time limit 4 II.Stochastic entropy production 4 III.Comparison with the dissipated power predicted by a generalized Langevin equation 5 IV.Brownian particle and Gaussian field: independent processes 6 A. Brownian particle in a harmonic potential 7 B. Dynamics of the free field 7 V. Derivation of the moment generating function of the particle position 8 VI.Average particle position and dissipation rate with non-conserved dynamics 11 A. Typical time and length scales 11 B. Behavior as a function of the drag speed v 12 C.The case of critical non-conserved dynamics 15 VII.Numerical simulation of the stochastic dynamics 15 A. Numerical integration scheme 16 B. Choice of parameters 17 C. Numerical measurement of the distribution of the particle position 18 D. Numerical measurement of the heat dissipation field 19 References 20

p
(u) differ from the corresponding χ p (u) and C p (u) only by a multiplicative factor exp(ip • v u) for u > 0 [see Eqs. ( (120)] of the shear rate to the trap time scale which, based on the estimates above for ρ and w, is typically of O(1).The experimentally accessible values of the parameters in Eq. (120) discussed above are reproduced in numerical simulations by choosing the following set of values of the model parameters:T = 0.7, κ = 0.2, v = 5, R = 2, λ = 10, r = 10 −4 , γ −1 φ = 20, and γ −1 y = 5, (

Fig. 1 .
In particular, for a field in d = 2 with non-conserved dynamics, the upper panel of Fig 1 presents a scatter ). ⟨φ(z)⟩st decays exponentially upon increasing |z| with different decay lengths ℓ± in front of or behind the particle -see the main text and Ref.[64].We used Vq = exp −q 2 R 2 /2 and R 2 + T /κ = 1, while the other parameters were set to unity.(b)Additionaldissipationrate[see Eq. (21)] as a function of v, for various values of ξ (with T = 0.1, while the other parameters were set to one).Symbols correspond to simulations.Inset: scaled correction − ⟨Z⟩ /λ 2 to the average particle position [see Eq. (19)].