Stochastic energetics of a colloidal particle trapped in a viscoelastic bath

We investigate the statistics of the fluctuations of the energy transfer between an overdamped Brownian particle, whose motion is confined by a stationary harmonic potential, and a surrounding viscoelastic fluid at constant temperature. We derive an analytical expression for the probability density function of the energy exchanged with the fluid over a finite time interval, which implicitly involves the friction memory kernel that encodes the coupling with such a non-Markovian environment, and reduces to the well known expression for the heat distribution in a viscous fluid. We show that, while the odd moments of this distribution are zero, the even moments can be explicitly expressed in terms of the autocorrelation function of the particle position, which generally exhibits a non-mono-exponential decay when the fluid bath is viscoelastic. Our results are verified by experimental measurements for an optically-trapped colloidal bead in semidilute micellar and polymer solutions, finding and excellent agreement for all time intervals over which the energy exchange takes place.


Introduction
Understanding energy exchanges in micron-and sub-micron-sized systems subject to fluctuations, e.g.colloidal particles, living cells, biopolymers, molecular motors, and small electronic circuits, is of prime importance in many disciplines of natural and applied sciences [1].In the last couple of decades, a significant progress has been made towards achieving this objective thanks to the advent of stochastic thermodynamics, which extends macroscopic concepts such as heat, work and entropy production to the level of single stochastic realizations of such processes [1,2,3].On this basis, fundamental constraints for the fluctuations of these thermodynamic quantities [4,5,6,7,8], as well as fluctuation-dissipation relations around out-of-equilibrium states [9,10,11,12,13,14], have been derived and as such, they represent refinements of the classical principles of thermodynamics.These theoretical relations have been tested in various experiments in small systems [15,16,17,18,19,20,21], which confirmed their validity under rather general non-equilibrium conditions.More recently, stochastic thermodynamics has also proved extremely helpful in studying energy fluxes and entropy production in more complex settings, such as in active matter [22,23,24,25,26], and systems with anomalous diffusion [27,28], thus providing a strong theoretical framework to describe small-scale thermodynamic processes arbitrarily away from equilibrium.
For mesoscopic systems in contact with reservoirs at constant temperature, e. g. colloids, vesicles, and macromolecules in aqueous solution, energy is continuously transferred between the system and the bath due to thermal collisions with the surrounding molecules even in the absence of external driving forces [29].Then, in many situations it is not enough to only know the second-law bounds imposed on the energy flows or the entropy production by the fluctuations relations [4,5,6,7,8] but also the detailed shape of their corresponding probability distributions.This is of special interest for processes taking place during short time intervals over which fluctuations are expected to largely exceed the mean values of the thermodynamic quantity of interest.Along these lines, a number of investigations using stochastic thermodynamics have been carried out over the past years in order to have a grasp of the statistics of the energy exchanged as heat between a Brownian system and its environment in absence of applied work.For instance, Fokker-Planck equations for the probability density function of the heat have been derived in presence of arbitrary confining potentials, whose long-time asymptotic solutions were experimentally verified for colloidal particles trapped in water by optical tweezers [30].Analytical expressions for the probability distribution of the heat transferred during an arbitrary time interval for the same system have also been obtained using path integrals [31,32].Moreover, the heat probability distribution has been determined for overdamped Brownian particles in non-stationary states relaxing toward thermal equilibrium after a temperature quench [33,34].In addition, the statistics of the steady heat fluxes for systems in contact with two thermostats at different temperatures have been investigated for quantum harmonic oscillators [35,36], RC electric circuits [37,38], pairs of hydrodynamically-coupled colloidal particles [39], and harmonic networks [40,41,42,43].Other effects on the heat distribution for Brownian particles, e.g.nonlinear potentials [44,45,46], inertia [47,48,49,50,51], relativistic motion [52], and non-isothermal transformations [53], have also been addressed theoretically.Furthermore, heat fluctuations have been studied for active matter systems, such as active chains in viscous heat baths [54], Brownian particles embedded in active media [55,56], and activity-driven harmonic chains [57].
It is worth mentioning that in most of the investigations on the heat probability distribution for thermostatted Brownian systems, there is a clear separation between the time scales of the system and those of the bath, which results in an effective Markovian description of the slow degrees of freedom of the system, for which the derivation of the heat distribution is straightforward.Nevertheless, the Markovian property is not possessed by a broad diversity of soft matter systems ranging from dense colloidal suspensions to polymeric fluids, which exhibit slow relaxations due to their crowded macromolecular microstructure [58].This leads to memory effects on the dynamics of Brownian-particle systems dispersed in such media, which are commonly described at equilibrium by the generalized Langevin equation, where the environment acts as a non-Markovian bath [59].Surprisingly, to the best of our knowledge, there is only a single theoretical work where the calculation of the first two moments of the distribution of the energy exchanged between a Brownian particle in a harmonic trap and its viscoelasticfluid surroundings is carried out based on the generalized Langevin equation [60].This leaves open the question of what is the full shape of the probability distribution under such non-Markovian conditions for an arbitrary time interval over which the stochastic energy transfer takes place.We point out that finding such a probability distribution for a Brownian system in equilibrium with a single non-Markovian bath is a first step towards the detailed knowledge on the statistics of the energy exchange in more intricate non-equilibrium situations, e.g.under time-dependent driving forces, reservoirs at different temperatures, and gradients.In turn, such problems are of great relevance in many applications at mesoscopic scales, e.g, the micromanipulation of colloidal probes in soft materials [61,62], the operation of Brownian heat engines in complex media [63,64,65,66], and the controlled microswimming in viscoelastic environments [67,68,69,70].
In view of the above considerations, the main purpose of this paper is to get insights into the effect of the non-Markovianity of a heat bath on the stochastic energetics for small thermostatted systems.To this end, we focus on a model system, namely, a micron-sized bead trapped by stationary optical tweezers in a viscoelastic fluid kept at constant temperature, whose motion is governed by the generalized Langevin equation.Based on the characteristic functional of the stochastic thermal force acting on the particle, we are able to derive an analytical expression in the overdamped limit for the energy transferred between the particle and the surrounding fluid over an arbitrary time interval.Our theoretical results are verified by experimental measurements in viscoelastic fluids such as aqueous worm-like micellar and polymer solutions.
The paper is organized as follows.In section 2 we present the model for the motion of a Brownian particle subject to a harmonic potential in a viscoelastic fluid, from which we find explicit formulae for some statistical quantities that describe the stochastic dynamics at thermal equillibrium.Then, in section 3 we derive an analytical expressions for the probability density function of the energy transferred from the fluid to the particle, and discuss its connection with previously obtained expressions for similar systems in the case of a Markovian heat bath.In section 4 we describe the experimental setup that we use to analyze the stochastic energetics of a colloidal bead harmonically trapped in some complex fluids with different viscoelastic responses, and make a comparison with our analytical results.Finally, in section 5, we summarize our main results and make some further physical remarks.

Model
We consider a spherical Brownian particle of mass m and radius a, which is embedded in an incompressible viscoelastic fluid medium with stress relaxation modulus G(t) and mass density ρ at constant temperature T .The viscoelastic fluid might be composed of macromolecules, e.g.polymer chains, micelles, or colloidal nanoparticles, dispersed in a viscous solvent, whose relaxation modulus is generally a slowly decaying function of t such that G(t → ∞) = 0, i.e. it behaves as a liquid in the long-time limit [58].We focus on the motion of a single coordinate of the center of mass of the particle, x, which is subject to a stationary harmonic potential of constant stiffness κ, i.e.U(x) = 1 2 κx 2 , as commonly implemented in experiments by means of optical tweezers.At thermal equilibrium, x evolves stochastically in time according to the generalized Langevin equation [59,71] The left-hand side of (1) corresponds to the inertial force on the particle, where x(t) is its instantaneous position at time t ≥ 0 starting from the initial conditions x 0 = x(0) and v 0 = ẋ(0) at time t = 0.This inertial term involves the effective mass m eff = m + 2 3 ρπa 3 that includes, in addition to the particle mass m, half of the mass of the fluid displaced by the particle [72].Besides, the first term on the right-hand side of (1) represents the coarse-grained drag force exerted on the trapped particle at time t by all the surrounding fluid particles.Here, Γ(t − t ′ ) is a memory kernel that encodes the delayed effect of the fluid at times 0 ≤ t ′ ≤ t, with Γ(t − t ′ ) = 0 for t ′ > t by causality.Furthermore, −∂ x U(x)| x=x(t) = −κx(t) corresponds to the value of the conservative force at time t that derives from the harmonic potential U(x).Moreover, ζ(t) is a Gaussian colored noise that accounts for the thermal collisions of the fluid particles, whose mean and autocorrelation function satisfy at thermal equilibrium respectively.The Laplace transform of the memory kernel, Γ(s) = ∞ 0 dt e −st Γ(t) is related to the fluid viscoelasticity by the expression [73] Γ(s) = 6πaη(s) + 6πa 2 ρsη(s). ( The first term on the right-hand side of (3) corresponds to conventional Stokes law for the drag force on a sphere moving in a fluid with a viscosity η(s) that is dependent on the Laplace frequency, s, which is determined by the Laplace transform of the fluid's relaxation modulus, η(s) ≡ G(s) = ∞ 0 dt e −st G(t), whereas the second term represents the Basset-Boussinesq force that originates from the motion of the displaced fluid [74,75].
We focus on the overdamped limit of (1), i.e. m eff , ρ → 0, which is a very good approximation in many colloidal experiments where a ∼ 10 −6 m, in such a way that the accessible frequencies to detect the single-particle motion are commonly |s| ≪ | Γ(s)|/m eff ∼ 10 6 rad s −1 and |s| ≪ |η(s)|/(ρa 2 ) ∼ 10 6 rad s −1 .In such a case, inertial terms can be totally neglected in the particle dynamics and, consequently, Eq. ( 1) becomes where the Laplace transform of the memory kernel Γ(t) can be approximated as Note that in the case of a Newtonian fluid of frequency-independent viscosity η(s) = η, the relaxation modulus is G(t) = 2ηδ(t), where δ(t) is the Dirac delta function.This yields the well-known overdamped Langevin equation for a Brownian particle trapped in a viscous medium where γ = 6πaη is the friction coefficient and ζ(t) reduces to a zero-mean deltacorrelated Gaussian noise, i.e.
By taking the Laplace transform of each term in (4), for the initial condition x(0) = x 0 the solution of the Laplace transform of x(t), i.e. x(s) = ∞ 0 dt e −st x(t), reads x(s) = χ(s)x 0 + ς(s) ζ(s), (7) where the functions χ(s) and ς(s) are defined in Laplace domain by the expressions respectively.Therefore, the solution in time domain, x(t), for a given stochastic realization of the thermal noise ζ(t ′ ) during the interval 0 ≤ t ′ ≤ t, can be obtained by Laplace inversion of ( 7) where χ(t) and ς(t) are the inverse Laplace transforms of the functions χ(s) and ς(s) defined in (8) and (9), respectively.Since we are interested in the situation where the particle is in equilibrium with the viscoelastic bath, at all times t ≥ 0 the stationary probability density function of the particle position x must be the canonical one Equation ( 11) is consistent with the equipartition relation for the particle position, x(t), at any time t ≥ 0, namely, 1 2 κ x(t) 2 eq = 1 2 k B T , where . . .eq denotes an average with respect to the stationary distribution (11).Therefore, from (10) and (11) the following expression can be readily derived for the autocorrelation function of x(t) computed between times t ′ and t ′ + τ where t ′ , τ ≥ 0. Note that (12) shows that the function χ(t) is proportional to the autocorrelation function of the particle position at thermal equilibrium.Moreover, regardless of the specific form of the memory kernel, χ(t) generally exhibits the two limit values In the second equality of ( 13), we have used the fact that lim s→0 s Γ(s) = lim t→∞ Γ(t) = 0 because the stress relaxation modulus of a viscoelastic fluid satisfies G(t → ∞) = 0. Furthermore, by computing the average of the square of each side of (10) over an infinite number of realizations of the thermal noise ζ(t ′ ) for a fixed value of x 0 , then with respect to the stationary distribution P eq (x 0 ) given in (11), and finally multiplying by κ/2, the use of the equipartition relation leads to the equality where . . .denotes the ensemble average over ζ(t ′ ) during the time interval 0 < t ′ ≤ t.
We now proceed to determine the probability density of finding the particle at position x at time t > 0 provided that it was located at position x 0 at time t = 0, that we denote as P (x, t|x 0 , 0).As described in section 3, this is necessary for the calculation of the probability density function of the variation of the harmonic potential energy U(x) at two different times.By definition, this conditional probability is where . . .0 denotes an ensemble average over an infinite number of realizations of the thermal noise ζ(t ′ ) during the time interval 0 < t ′ ≤ t for a fixed initial condition x(0) = x 0 .In the second line of (15), we have used the Fourier representation of the Dirac delta function, , as well as the expression of the time-domain solution x(t) with initial condition x(0) = x 0 given in (10).Note that the term Φ in the second line of (15) can be identified as the characteristic functional of the noise ζ(t ′ ).Since ζ(t ′ ) is assumed to be Gaussian, Φ[−kς] can be simply expressed as [76] Φ where we have made us of ( 14) in the second line of (16).A direct substitution of ( 16) into ( 15) and a straightforward calculation of the corresponding inverse Fourier transform leads the final expression of the conditional probability density P (x, t|x 0 , 0) We point out that the same expression for P (x, t|x 0 , 0) as in (17) can also be obtained by solving a non-Markovian Smoluchowski equation equivalent to the full generalized Langevin equation ( 1) with inertia provided that the system starts from an initial equilibrium state, i.e. with Boltzmann distribution P eq (x for the initial conditions (x 0 , v 0 ) [77,78], and then taking the overdamped limit m eff → 0. In either case, the assumption of initial equilibrium at t = 0 guarantees the stationarity of the probability density of x at all times t > 0, i.e. ∞ −∞ dx 0 P (x, t|x 0 , 0)P eq (x 0 ) = P eq (x), as can be verified by direct use of ( 17) with both P eq (x 0 ) and P eq (x) given by (11).

Probability distribution of the energy exchanged with the bath
In this section, we derive an analytical expression for the probability density function of the energy transferred during a given time interval from the viscoelastic fluid, which plays the role of a non-Markovian bath at constant temperature, to the mesoscopic system of interest, which is in this case Brownian particle trapped by the harmonic potential. .In situations like this, during a short time interval of duration of at least dt ∼ 10 −6 s, the variation of the bare Hamiltonian of the system, H(x, ẋ) = 1  2 m ẋ2 + U(x) = 1 2 m ẋ2 + 1 2 κx 2 , i.e. its energy without taking into account the interaction energy with the bath particles, can be approximated as dH = (mẍ + κx) ẋdt ≈ κx ẋdt = dU.This is because for such values of dt, the overdamped limit is valid, in such a way that the inertial term mẍ can be completely neglected with respect to the harmonic force, −κx.Furthermore, in absence of external forces performing work on the system, we can interpret the variation dH = dU as the energy randomly exchanged with the viscoelastic fluid during dt.If dH > 0, the systems absorbs energy from the surroundings, whereas it releases energy if dH < 0, thus increasing or decreasing its potential energy, respectively.Therefore, the energy exchanged between the system and the bath during a time interval [0, τ ] at the level of a single stochastic trajectory x(t), with 0 ≤ t ≤ τ , starting at x(0) = x 0 and ending at x(τ ) = x τ , which we denote as Q τ , can be determined from It should be noted that (18) has the same form as the first-law-like relation for the energy balance along a stochastic trajectory of a Brownian particle harmonically trapped in thermal equilibrium in a viscous solvent.In such a case, the left-hand side represents the heat stochastically exchanged with the bath via the fast thermal collisions with the solvent molecules and by viscous friction, and exactly amounts to minus the total energy variation of the environment.However, when the particle is trapped in a viscoelastic fluid, Q τ cannot be simply identified as heat since the stochastic energetics of the system also involves the interaction energy between the trapped particle and the surrounding macromolecules suspended in the solvent.The reason of this is that such interactions are not negligible and give rise to temporal correlations that do not relax as quickly as in the case of a purely viscous bath, as manifested by the presence of slowly decaying memory kernels and colored noises in the coarse-grained description of the particle dynamics provided by the generalized Langevin equation (1).In this situation, a more detailed description of the energy balance at the system trajectory level could be achieved by introducing the concept of Hamiltonian of mean force, which includes the effect of the interaction energy averaged over the degrees of freedom of the bath particles [79,80,81].
As discussed in [81,82], this approach could eventually allow to properly define the exchanged heat for mesoscopic systems that strongly interact with a bath such as the viscoelastic fluid considered here, even though there is also some criticism against the unambiguous identification of stochastic energy-like quantities other than the applied work under such conditions [83,84].Therefore, in the following, we will not address these ongoing issues on stochastic thermodynamics for systems that are strongly coupled to a heat bath and will only focus on the statistics of Q τ .In absence of work externally applied to the system, this quantity can be univocally identified as the energy absorbed from or released to the surrounding viscoelastic fluid by the particle, which allows it to explore the harmonic potential beyond the minimum at x = 0. Indeed, by multiplying the overdamped version of the generalized Langevin equation (1) by a small particle displacement dx = ẋ(t)dt performed between times t ≥ 0 and t+dt, and upon integration over [0, τ ] using (18), Q τ can be expressed as where the first term on the right-hand side corresponds to the energy supplied to the system by the fluctuating thermal force of the bath, ζ(t), during the stochastic realization of x(t) from t = 0 to t = τ , whereas the second term represents the energy exchanged by the action of the hydrodynamic drag exerted by the fluid, ), both terms including the effect of the interaction of the trapped particle with the macromolecules of the fluid.Equation ( 19) correctly reduces to the widely accepted definition of stochastic heat in the case of a purely viscous bath of constant viscosity η [2], in which case Γ(t) = 12πaηδ(t), thus yielding where γ = 6πaη is the friction coefficient of the particle.Note that the second terms on the right-hand sides of ( 20) and (19) highlight the subtle difference between the properties of the heat exchange that is properly defined for a particle weakly coupled to a viscous bath with respect to the general energy transfer when it is strongly coupled to a viscoelastic one.While in the former case viscous friction always produces energy dissipation into the bath, which corresponds to the negative definite term of (20), in the latter the viscoelastic drag could give rise to either a positive or a negative contribution to the energy transfer because of the long-lived particle interactions comprised in the memory kernel in (19).
Once the physical interpretation of Q τ has been clarified, the calculation of its probability distribution that we denote as P (Q, τ ), can be carried out in a straightforward manner.We realize that, according to (18), Q τ is only determined by the difference of the two stochastic variables x 2 0 and x 2 τ , which are in general correlated for finite values values of τ .Therefore, P (Q, τ ) can be expressed as where P eq (x 0 ) is the probability density function of the initial position x 0 , which is given by ( 11) since we assume that the system is at all times t ≥ 0 in thermal equilibrium with the viscoelastic bath, In addition, the function P (x τ , τ |x 0 , 0) in ( 21) corresponds to the probability density of finding the particle at position x τ at time τ > 0 provided that it was located at x 0 at time t = 0, whose explicit expression is given in (17) for any arbitrary memory kernel Γ(t − t ′ ) governing the stochastic dynamics of x(t).Therefore, by taking into account that the Dirac delta function in the integrand of ( 21) can be written as and by substituting ( 22) into (21) in order to integrate over −∞ < x τ < ∞, a more compact expression for P (Q, τ ) can be obtained where the functions P + (Q, τ ) and P − (Q, τ ), defined as only involve a single integration over −∞ < x 0 < ∞.Moreover, by means of the change of variables and by introducing the functions the integrals in (24) can be recast as Furthermore, by use of the identity where K 0 (w) = ∞ 0 du cos (w sinh u) denotes the zeroth order modified Bessel function of the second kind.Equation (28), which is the main result of the present paper, provides a surprisingly simple expression valid in the overdamped regime for the probability distribution of the energy exchanged in thermal equilibrium during a time interval of duration τ 10 −6 s between the harmonically trapped particle and the surrounding viscoelastic bath.Note that (28) reveals that all the information on the energetic coupling between the particle and the viscoelastic bath is encoded in a single function, namely, ϕ + (τ )ϕ − (τ ) = 1−χ(τ ) 2 , whose dependence on τ is determined by the strength of the harmonic potential and by the specific form of the relaxation modulus of the bath, as will be shown in section 4 for two specific examples of viscoelastic fluids.Moreover, since the function χ(t) is proportional to the autocorrelation function of the particle position at thermal equilibrium according to (12), the function ϕ + (τ )ϕ − (τ ) contains information of the temporal correlations induced by the energetic coupling of the trapped particle with its non-Markovian environment.
From the explicit expression given by ( 28) and using the definite-integral identity for the µ−th order modified Bessel function of the second kind ‡, all the moments of order n ≥ 0 of the distribution of P (Q, τ ), i.e.
can be calculated.Indeed, owing to the symmetry of P (Q, τ ) with respect to Q = 0 In particular, from (30) we can easily check that P (Q, τ ) is normalized, i.e.
∞ −∞ dQ P (Q, τ ) = 1.Additionally, the mean energy transferred from the bath to the particle is which simply reflects the fact that detailed balance guarantees the absence of net energy fluxes on average at thermal equilibrium even if the bath is non-Markovian, like the viscoelastic fluid at constant temperature considered here.Moreover, from (30), ‡ Here, Γ (z) = ∞ 0 du u z−1 e −u represents the well-known gamma function, which must not be confused with the memory kernel Γ(t − t ′ ) describing the Brownian particle dynamics in (1).
the variance of the exchanged energy during a time interval of duration τ , ∆Q whose asymptotic value for τ → ∞ is ∆Q 2 τ →∞ = (k B T ) 2 , according to the limit values of the function χ(t) given in (13) for an arbitrary relaxation modulus of the viscoelastic fluid.This result shows that, for sufficiently long intervals over which all temporal correlations induced either by the trapping potential or by the interactions with the medium vanish, the characteristic exchanged energy is ( ∆Q regardless of the specific features of the memory kernel and the value of the trap stiffness. It should be noted that expressions for Q τ and ∆Q 2 τ that are similar to those given in (31) and (32) were derived in [60] for the same model investigated here based on the general solution of (1), [x(t), ẋ(t)], for given initial conditions (x 0 , v 0 ) at t = 0.However, the calculation of higher-order moments following this approach becomes too involved, thus being feasible in practice only for the determination of the mean and the variance.Here, we overcome this limitation by directly deriving an expression for the probability distribution of Q τ , thereby having full knowledge of its statistical properties.For instance, from (30), the skewness of P (Q, τ ) is 0, whereas the kurtosis is independent of τ and has the constant value Some additional remarks can be made about (28).For instance, it can be easily checked that it correctly reduces to the well know expression for the probability distribution of the heat exchanged between a spherical particle of radius a confined by a harmonic potential and a viscous heat bath of constant viscosity η.Indeed, by taking the instantaneous memory kernel Γ(t − t ′ ) = 2γδ(t − t ′ ), which corresponds to a timeindependent friction coefficient γ = 6πaη, the functions in ( 5) and (8) become Γ(s) = γ, and χ(s , respectively, where τ γ = γ/κ is the characteristic relaxation time of the particle that originates from the balance between viscous friction and the restoring force exerted by the harmonic potential.In this case, the inverse Laplace transform of χ(s) is χ(t) = exp − t τγ , so that ϕ + (τ )ϕ − (τ ) = 1 − e −2τ /τγ = 2e −τ /τγ sinh(τ /τ γ ), thus leading to the following expression for the probability density function of the exchanged heat which is the same formula as previously derived in [31] using path integrals for the Markovian dynamics corresponding to the Langevin equation (6).Another interesting property of ( 28) is that, irrespective of the specific details of the fluid viscoelastic properties and of the strength of the trapping potential, in the long-time limit τ → ∞ it tends to where we have used the asymptotic value χ(t → ∞) = 0 of the function χ(t) given in (13).It is noteworthy that ( 35) is identical to the formula derived in [30] for long-time behavior of the probability distribution of the heat exchanged between a Brownian particle trapped by a stationary harmonic potential and a viscous fluid bath.
In addition, asymptotic expressions similar to (35) have also been derived for the heat probability distribution of harmonically trapped particles in non-stationary [33,34] and active media [55], for particles confined by non-linear potentials at low temperature [44], as well as for the fluctuations of the energy current through harmonic chains connected to two active reservoirs [57].This observation suggests that the zeroth order modified Bessel function of the second kind captures the asymptotic behavior of the energy exchange process with the environment of a broad variety of overdamped Brownian systems provided that the bath exhibits fluid-like behavior in the long-time limit.

Experimental results
In this section, we experimentally verify the main results that were previously derived in section 3. To this end, we use as viscoelastic baths two kinds of complex fluids with distinct rheological properties.The first fluid is a worm-like micellar solution, which is composed of the surfactant cetyltrimethylammonium bromide (CTAB) and the salt sodium salycilate (NaSal) at equimolar concentration of 2.5 mM in ultrapure water (resistivity 18.2 MΩ cm at 25 • C).At this concentration, which is above the critical micelle concentration of pure CTAB (0.83 mM), the surfactant molecules self-organize in flexible cylindrical micelles that have a radius of ∼ 2 nm, an average contour length of ∼ 200 nm, and a persistence length of ∼ 50 nm [86].The second fluid consist of a solution of the water-soluble polymer polyethylene oxide (PEO, molecular weight M w = 4 × 10 6 Da) dissolved in ultrapure water at a concentration of c = 0.1 %wt, which is smaller than the critical overlap concentration (c * ≈ 0.65 %wt) and corresponds to the semi-diluted particle regime in which polymer chains are non-entangled but inter-chain interactions are not negligible [87].For the sake of simplicity, in the following these fluids will be simply referred to as micellar solution and polymer solution, respectively.A complete homogenization of both solutions is achieved by continuous stirring over 24 h, thus becoming transparent to visible light and endowed with viscoelasticity due to the presence of the two types of macromolecules that were previously described.
Once the homogeneous micellar and polymer solutions are prepared, a small amount of spherical silica beads (radius a = 1 µm) are dispersed in them, thus resulting in viscoelastic fluid suspensions at very low concentration (less than 1 colloidal bead in 1 nl of solution).Sample cells made of a microscope glass slide and a coverslip stuck together parallel to each other by double-sided adhesive tape (thickness ∼ 100 µm), are filled with the viscoelastic suspension of interest, and sealed with epoxy glue in order to prevent leakage and evaporation of the fluid.Then, a green Gaussian laser beam (wavelenght 532 nm) is tightly focused inside the sample cell using an oil-immersion objective (100×, numerical aperture NA = 1.3), which allows one to trap a single bead at h ≈ 20 µm above the lower solid wall of the cell, as sketched in figure 1(a).Such an optical trap confines the stochastic particle motion within a harmonic potential, whose stiffness κ is kept constant for each fluid by fixing the laser power.Videos of the trapped particle are recorded using a CMOS camera at a sampling rate of f s = 2000 frames per second during approximately t max = 25 min.The position (x, y) of the center of mass of two-dimensional projection in the x-y plane of the particle is detected with a spatial resolution of 5 nm using standard particle-tracking routines, as illustrated by the snapshot shown in the upper panel of figure 1(b).The measurements are performed at temperatures of T = 23.5 ± 0.2 • C and T = 21.0 ± 0.2 • C for the micellar and polymer solutions, respectively.
In the lower panel of figure 1(b) we show an example of the stochastic time evolution of a single coordinate, e.g.x(t), of a particle trapped in the micellar solution, which is recorded over approximately 25 minutes.The trap stiffness can be computed by, e.g. a quadratic-polynomial fit of the potential U(x) = 1  2 κx 2 that is obtained from the Boltzmann distribution of the particle position given by ( 11), as shown in figure 1(c).An alternative method based on the equipartition relation for x is described in Appendix A, and leads to consistent values for κ.The same procedure is followed to characterize the trapping potential in the case of the polymer solution.Also, as described in Appendix A, by means of the numerical calculation of the Laplace transform of the mean-squared displacement of x(t), we can determine the Laplace-frequency dependent viscosity of the fluid, η(s), where we only consider positive real frequencies, s > 0, for the sake of clarity.In figure 1(d) we plot the dependence of η(s) on s for both the the micellar and the polymer solution.We note that in both cases, such dependencies can be fitted to the function whose functional form is consistent with the so-called Carreau-Yasuda equation that is commonly used as a phenomenological model for the shear-rate dependent viscosity of non-Newtonian fluids [88].In (36), η 0 and η ∞ correspond to the values of the fluid viscosity at vanishing and infinite frequencies, respectively, λ is the characteristic relaxation time of the fluid, whereas α ≥ 0 and n ≤ 1 are exponents that depend on the specific microscopic features of the fluid and describe a power-law-like decay of η(s) −η ∞ ∝ s −(1−n) at intermediate frequencies.The values of the parameters η 0 , η ∞ , λ, α and n for each fluid along with the values of κ for the trapping harmonic potentials in each experiment are listed in table 1.Note that in the case of the polymer solution, we have restricted the outer exponent in (36)  the literature [87].On the other hand, both n and α have been employed as independent free fitting parameters for the micellar solution.It should also be noted in figure 1(d) that the dependence of the viscosity on frequency is much more pronounced for the micellar solution than for the polymer solution, where the characteristic relaxation time λ of the former is three orders of magnitude larger than that of the latter.Such a marked difference in the rheological properties of these fluids allows us to analyze the stochastic energy exchange of the trapped particle with two types of non-Markovian baths, one with rather strong viscoelastic behavior (the micellar solution) and the other Table 1.Parameters characterizing the fluid viscoelasticity and the optical harmonic potential for each experiment.with much weaker viscoelasticity (the polymer solution).
From the experimental particle trajectories, we now proceed to compute the probability density function of the energy exchanged between the trapped particle and its corresponding viscoelastic surroundings, as well as its respective standard deviation and kurtosis for time intervals of different duration τ , typically f −1 s ≤ τ 10 s, where f −1 s = 5 × 10 −4 s.For this purpose, for a fixed value of τ we directly use (18) to compute stochastic values of Q τ , whose normalized histogram allow us to determine the probability density function P (Q, τ ).Since the system is at thermal equilibrium, its stationarity allows us to choose the initial time in (18) as any time t ′ along the recorded trajectory such that 0 In this way, the probability density function P (Q, τ ) and its corresponding moments are computed over approximately 3 × 10 6 data points of the random variable Q τ for each experiment.
In figure 2(a) we show the results for the micro-bead trapped in the micellar solution, where we verify that the probability density function of Q τ is accurately described by the zeroth order modified Bessel function of the second kind given in (28) over six orders of magnitude of P (Q, τ ) for all the explored values of τ spanning 5 order of magnitude (10 −3 −10 s).To this end, for each value of the duration τ the value of χ(τ ) is computed by a numerical Laplace inversion of the expression given in (8) by means of the Talbot method [90], where we use the experimental value of the trap stiffness as well as the functional form for η(s) given by (36), with the corresponding values of its parameters: η 0 = 0.0391 Pa s, η ∞ = 0.0009 Pa s, λ = 11.910s, α = 2.125, and n = 0.541.We also verify that, as τ increases, the experimental shape of P (Q, τ ) becomes independent of τ and tend to the function given in (35), as can be checked in figure 2(a) for τ ≥ 10 s, at which the curve of P (Q, τ ) is practically indistinguishable from P (Q, τ → ∞).In figure 2(b) we plot the dependence of the function χ(t) on time t, which was numerically determined as previously described using the Talbot method.For comparison, we also plot the behavior of χ(t) in the case of Newtonian fluids with viscosities η ∞ and η 0 for the same value of the trap stiffness.We observe that, unlike in a Newtonian fluid, where χ(t) is merely an exponential decay with characteristic time given by the ratio between the corresponding friction coefficient of the particle and the trap stiffness, in the micellar solution χ(t) exhibits an intricate non-monoexponential behavior as a function of t.For instance, for the specific value of κ used in the experiment, in the micellar solution χ(t) decays slower than in the Newtonian fluid with viscosity η ∞ but faster than in the case η 0 up to t ≈ 4 s.However, for .Solid lines represent the analytical expression given in (28) for the corresponding values of τ , whereas the black dashed line depicts the limiting expression for τ → ∞ given by (35).(b) Dependence of χ(t) on time t for the same system (solid line), which was obtained by numerical inversion of (8) using the functional form (36) for the frequency-dependent viscosity of the micellar solution shown in figure 1 t 4 s, the non-Markovian nature of the micellar fluid bath gives rise to a very slow decay of χ(t) that becomes even slower than the exponential decay e −κt/(6πaη 0 ) , then finally approaching χ(t) → 0 as expected in any fluid as t → ∞.To highlight the role of the viscoelasticity of the micellar solution on the energy exchange process with the particle, in figure 2(c) we plot as symbols the experimentally-determined standard deviation of Q τ , ∆Q 2 τ , as a function of the time interval duration τ .We verify that the square root of the theoretical expression of the variance given in (32), which is computed using the numerical curve of χ(t) plotted in figure 2(b), agrees very well with the experimental data, as shown by the solid line in figure 2(c) for all the values of τ investigated here.In the same plot, we can also check experimentally that the standard deviation of the energy exchange converges to the limit value k B T as τ → ∞, in accord with the expected asymptotic properties of the energy exchange fluctuations derived at thermal equilibrium from the generalized Langevin model (1).Once again, for comparison purposes, we also plot in figure 2(c) the theoretical predictions of (32) for the standard deviation of Q τ as a function of τ for the heat exchange in Newtonian fluids of viscosities η 0 (dashed line) and η ∞ (dotted-dashed line).We find that, for 0 < τ 4 s, the behavior of ∆Q 2 τ in the micellar solution is comprised between those in Newtonian fluids of viscosities η 0 and η ∞ , while for τ 4 s, the convergence to the limit value k B T becomes slower in the viscoelastic micellar fluid than in both Newtonian cases.Interestingly, we observe in the inset of figure 2(c) that even though none of these Newtonian curves describes correctly the full dependence of the standard deviation of Q τ on τ in the micellar solution, at small values of τ its behavior is very close that in a Newtonian fluid with viscosity η ∞ , i.e.
∆Q 2 τ ≈ k B T 1 − e −2κτ /(6πaη∞) .This observation can be explained by the fact that, during sufficiently short time intervals, the energy exchange process between the particle and the viscoelastic bath is dominated by the random heat transfer due to the fast collisions of the solvent molecules, which is characterized by the high-frequency friction coefficient 6πaη ∞ .For longer time intervals, interactions between the trapped bead and the surrounding wormlike micelles become important, thus leading to a transient energy storage in the micelles, which translates into systematic deviations of the standard deviation of Q τ from the Newtonian curve with viscosity η ∞ .As τ becomes much larger than the characteristic correlation times encoded in the behavior of χ(t), the stochastic energy exchange becomes independent of both the trapping potential and the non-Markovian properties of the viscoelastic micellar bath, thereby resulting in the convergence of the standard deviation of Q τ to the characteristic scale of thermal energy, k B T .We point out that, even though we verify experimentally that the standard deviation of Q τ depends on τ through the function χ(τ ), the kustosis remains independent of τ , very close to the constant value K τ = 9 that is predicted by (33), as shown in figure 2(d).This confirms that for all the values of τ analyzed here, the experimental shape of P (Q, τ ) is precisely described by the zeroth order modified Bessel function of the second kind given in (28).
The results for the probability density function of the energy exchange in the polymer solution are shown in figure 3(a).In the case of this weakly viscoelastic fluid, The solid line depicts the constant value K τ = 9 predicted by (33).
we also corroborate that the theoretical expression for this probability distribution given in ( 28) is in excellent agreement with the experimental results over six orders in magnitude for the values of P (Q, τ ) and at least 4 orders of magnitude in τ .We find that for τ ≥ 5 s, P (Q, τ ) has already converged to the τ -independent function given in (35), see the black dashed line in figure 3(a), in which case all the temporal correlations in the energy exchange process induced by the trapping potential and the bath vanish.In figure 3(b) we plot the behavior of the function χ(t) for the particle trapped in the polymer solution (solid line), and the hypothetical cases of the same particle trapped in Newtonian fluids with viscosities η 0 = 0.0090 Pa s (dashed line), and η ∞ = 0.0017 Pa s (dotted-dashed line).We observe that, because of the weak viscoelastic behavior of the polymer solution manifested through the very smooth frequency-dependence of its viscosity η(s) plotted in figure 1(d), the dependence of χ(t) on t in the viscoelastic case is rather close to e −κ/(6πaη 0 ) .However, the non-Markovianity of the bath induced by presence of the suspended polymers gives rise to deviations from such a mono-exponential decay, which can be seen upon closer inspection of this plot and of the semi-log representation shown in the inset of figure 3(b).Indeed, it turns out that the behavior of χ(t) in the polymer solution is qualitatively similar to that in the micellar fluid.Specifically, at short times t 2.6 s, the values of χ(t) in the polymer solution lie between those in Newtonian fluids of viscosities η ∞ and η 0 , whereas they exceed both Newtonian cases for t 2.6 s, thus revealing the presence of slowly decaying correlations in the dynamics of the trapped particle even if the fluid bath is only weakly viscoelastic.Then, using the numerical curve χ(t) shown in figure 3(b), we can calculate the dependence of the standard deviation of Q τ , ∆Q 2 τ , on the time interval duration τ given by (32), which is plotted as a solid line in figure 3(c).Again, we observe a conspicuous agreement between this prediction and the standard deviation of the energy exchange computed from the experimental data points for values of τ spanning five orders of magnitude (5 × 10 −4 − 6 s).Although the viscoelasticity of the polymer solution is not as pronounced as that of the micellar one, we are able to detect similar qualitative features of the dependence of ∆Q 2 τ on τ .In particular, we find in the inset of figure 3(c) that at sufficiently small values of τ , ∆Q 2 τ ≈ k B T 1 − e −2κτ /(6πaη∞) due to the dominance of the heat exchange with the solvent molecules, followed by an intermediate regime where deviations from the behavior in a Newtonian fluid occur due to the slow temporal correlations induced by interactions with the polymer chains, whereas for sufficiently large τ we observe the convergence to the asymptotic value ∆Q 2 τ →∞ = k B T .Finally, in figure 3(c) we verify that the values of the kurtosis of the experimentally determined distributions for the energy exchange in the polymer solution are very close to the theoretical one K τ = 9.

Summary and concluding remarks
In this paper, we have investigated theoretically and experimentally the statistics of the stochastic energetics of an overdamped Brownian particle trapped at constant temperature in a viscoelastic fluid acting as a non-Markovian thermal reservoir.By analyzing the simplest nontrivial situation where the trapping potential is harmonic and no work is applied to the system but the fluid has an arbitrary relaxation modulus, we were able to derive an analytical expression for the probability density function of the energy exchanged between the particle and its surroundings during a finite time interval.We have shown that this probability distribution, which can be written as a modified Bessel function of the second kind with even symmetry, as well as all its even moments, can be expressed in terms of a function that depends on the timeinterval duration and is proportional to the squared autocorrelation function of the particle position at thermal equilibrium.Such expressions include the effect of temporal correlations induced by both the confining potential as well as the interaction with the bath particles, thereby representing an extension of previously derived formulae for the heat exchange of Brownian systems in Markovian baths, e.g.viscous fluids.Our results are clearly illustrated in experiments using colloidal particles optically trapped in two types of viscoelastic fluids with distinc rheological and structural properties, namely, a wormlike micellar solution and a semidilute polymer solution, which are of widespread interest in soft matter systems.
In recent years, there has been a growing interest in extending the theoretical framework of stochastic thermodynamics beyond the usual assumption of weak coupling with a heat bath for small systems both in the classical and quantum realms, e.g.mesoscopic systems that interact with a non-Markovian bath like the viscoelastic one investigated here.Although the thermodynamic description of such systems subject to external driving forces are a focus of attention in various theoretical investigations [81,82,84,91,92,93,94], very few experimental studies have been conducted in a controlled fashion even close to thermal equilibrium [95,96,97].Therefore, the results presented in this paper represent a first step towards a quantitative characterization of stochastic energy exchanges of Brownian systems in non-Markovian environments and their connection with the slowly-decaying correlations that emerge due to the interactions with the bath particles.Further experimental and theoretical efforts based on this model system could help to investigate other relevant aspects on the stochastic energetics at strong coupling under more intricate conditions, e.g. in thermally-activated escape processes over energetic barriers with memory friction [98], the efficiency of cyclic Brownian engines operating in non-Newtonian fluids [63,65,66], and work fluctuations of colloids periodically or stochastically driven in viscoelastic media [99,100].

Figure 1 .
Figure1.(a) Schematic representation of a spherical bead of radius a, trapped at a distance h above the lower wall of the sample cell containing a viscoelastic fluid, which is achieved by means of optical tweezers realized by a tightly focused laser beam (green area).The corresponding coordinate system used in the data analysis is also sketched.(b) Upper panel: snapshot in the x-y plane of a bead trapped in the micellar solution, with the location of its center of mass (x, y) marked as an asterisk.The horizontal scale bar represents 1 µm.Lower panel: example of the stochastic time evolution of the coordinate x of the same particle trapped in the micellar solution.(c) Probability density function of the particle position, x, in thermal equilibrium with the micellar solution (vertical bars), experimental profile of the optical trapping potential (•), and corresponding quadratic fit (solid line), from which the trap stiffness, κ, is calculated.(d) Frequency-dependent viscosities of the polymer (•) and the micellar ( ) solutions, which are determined experimentally from the particle motion.Solid lines represent the best non-linear fitting of the experimental data to the model given in(36).

Figure 2 .
Figure 2. (a) Probability density function of the energy exchanged between a colloidal bead trapped by a harmonic optical potential (κ = 3.08 × 10 −7 N m −1) and the surrounding micellar solution, over time intervals of different duration: τ = 0.001 s (•), τ = 0.01 s ( ), τ = 0.1 s (▽), τ = 1 s (△), and τ = 10 s (⋄).Solid lines represent the analytical expression given in(28) for the corresponding values of τ , whereas the black dashed line depicts the limiting expression for τ → ∞ given by(35).(b) Dependence of χ(t) on time t for the same system (solid line), which was obtained by numerical inversion of (8) using the functional form(36) for the frequency-dependent viscosity of the micellar solution shown in figure1(d).The dashed and dotted-dashed lines represent χ(t) for a particle in Newtonian fluids of viscosities η 0 = 0.0391 Pa s and η ∞ = 0.0009 Pa s, respectively.Inset: semi-log representation of the main plot.(c) Variance of the experimentally-determined energy exchange as a function of the duration of the time interval over which it takes place (•).The solid line represents the theoretical prediction of (32) using the numerical curve plotted as a solid line in figure 2(b), whereas the dashed and dotted-dashed lines correspond to the variances for a particle in Newtonian fluids of viscosities η 0 = 0.0391 Pa s and η ∞ = 0.0009 Pa s, respectively.Inset: double-log representation of the main plot.(d) Kurtosis of the energy exchange distribution as a function of the duration of the time interval (•).The solid line depicts the constant value K τ = 9 predicted by(33).
Figure 2. (a) Probability density function of the energy exchanged between a colloidal bead trapped by a harmonic optical potential (κ = 3.08 × 10 −7 N m −1) and the surrounding micellar solution, over time intervals of different duration: τ = 0.001 s (•), τ = 0.01 s ( ), τ = 0.1 s (▽), τ = 1 s (△), and τ = 10 s (⋄).Solid lines represent the analytical expression given in(28) for the corresponding values of τ , whereas the black dashed line depicts the limiting expression for τ → ∞ given by(35).(b) Dependence of χ(t) on time t for the same system (solid line), which was obtained by numerical inversion of (8) using the functional form(36) for the frequency-dependent viscosity of the micellar solution shown in figure1(d).The dashed and dotted-dashed lines represent χ(t) for a particle in Newtonian fluids of viscosities η 0 = 0.0391 Pa s and η ∞ = 0.0009 Pa s, respectively.Inset: semi-log representation of the main plot.(c) Variance of the experimentally-determined energy exchange as a function of the duration of the time interval over which it takes place (•).The solid line represents the theoretical prediction of (32) using the numerical curve plotted as a solid line in figure 2(b), whereas the dashed and dotted-dashed lines correspond to the variances for a particle in Newtonian fluids of viscosities η 0 = 0.0391 Pa s and η ∞ = 0.0009 Pa s, respectively.Inset: double-log representation of the main plot.(d) Kurtosis of the energy exchange distribution as a function of the duration of the time interval (•).The solid line depicts the constant value K τ = 9 predicted by(33).

Figure 3 .
Figure 3. (a) Probability density function of the energy exchanged between a colloidal bead trapped by a harmonic optical potential (κ = 1.43 × 10 −7 N m −1) and the surrounding polymer solution, over time intervals of different duration: τ = 0.005 s (•), τ = 0.05 s ( ), τ = 0.5 s (▽), and τ = 5 s (△).Solid lines represent the analytical expression given in(28) for the corresponding values of τ , whereas the dashed line depicts the limiting expression for τ → ∞ given by(35).(b) Dependence of χ(t) on time t for the same system (solid line), which was obtained by numerical inversion of (8) using the functional form(36) for the frequency-dependent viscosity of the polymer solution shown in figure1(d).The dashed and dotted-dashed lines represent χ(t) for a particle in Newtonian fluids of viscosities η 0 = 0.0090 Pa s and η ∞ = 0.0017 Pa s, respectively.Inset: semi-log representation of the main plot.(c) Variance of the experimentally-determined energy exchange as a function of the duration of the time interval over which it takes place (•).The solid line represents the theoretical prediction of (32) using the numerical curve plotted as a solid line in figure3(b), whereas the dashed and dotted-dashed lines correspond to the variances for a particle in Newtonian fluids of viscosities η 0 = 0.0090 Pa s and η ∞ = 0.0017 Pa s, respectively.Inset: double-log representation of the main plot.(d) Kurtosis of the energy exchange distribution as a function of the duration of the time interval (•).The solid line depicts the constant value K τ = 9 predicted by(33).