# Work Extraction and Performance of Colloidal Heat Engines in Viscoelastic Baths

- Departamento de Sistemas Complejos, Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de México, Mexico

A colloidal particle embedded in a fluid can be used as a microscopic heat engine by means of a sequence of cyclic transformations imposed by an optical trap. We investigate a model for the operation of such kind of Brownian engines when the surrounding medium is viscoelastic, which endows the particle dynamics with memory friction. We analyze the effect of the relaxation time of the fluid on the performance of the colloidal engine under finite-time Stirling cycles. We find that, due to the frequency-dependence of the friction in viscoelastic fluids, the mean power delivered by the engine and its efficiency can be highly enhanced as compared to those in a viscous environment with the same zero-shear viscosity. In addition, with increasing fluid relaxation time the interval of cycle times at which positive power output can be delivered by the engine broadens. Our results reveal the importance of the transient behavior of the friction experienced by a Brownian heat engine in a complex fluid, which cannot be neglected when driven by thermodynamic cycles of finite duration.

## 1. Introduction

Historically, the study of heat engines has played a fundamental role in the general understanding of energy exchanges in macroscopic systems. For instance, the conception of the well-known Carnot cycle almost two centuries ago was motivated by the design of efficient engines capable of performing mechanical work by extracting energy from a hot reservoir and transferring heat to a cold reservoir, which finally led to the formulation of the second law of thermodynamics. Carnot theorem imposes a universal bound for the maximum efficiency that can be ideally achieved by any heat engine working reversibly in the quasi-static limit. Since then, further theoretical results on the efficiency of irreversible heat engines under finite-time thermodynamic cycles with non-zero power output have been obtained [1–6], which turn out to be important for practical applications.

In more recent years, advances in miniaturization technologies have allowed researchers in both basic and applied science to conceive the design of micron- and submicron-sized machines with the ability to perform specific tasks in the mesoscopic realm, e.g., controlled cargo transport through microchannels and nanopores, *in situ* cell manipulation, assembly of functional microstructures, micropumping, microflow rectification, micromixing of fluids, and bio-inspired artificial locomotion [7–9]. This has triggered an increasing interest in investigating the energetics and performance of mesoscopic heat engines, which, similar to their macroscopic counterparts, must be able to convert in an efficient manner the energy absorbed from their environment into useful work [10, 11]. An important issue that arises in the theoretical description and implementation of such devices is that they must operate under highly non-equilibrium conditions with pronounced thermal fluctuations, which poses important conceptual and practical challenges [12]. A significant progress in the theoretical analysis of mesoscopic heat engines has been made in the last two decades with the advent of stochastic thermodynamics, which extends concepts of classical thermodynamics such as heat, work and entropy production to the level of single stochastic trajectories for both equilibrium and driven systems [13–16]. Within this theoretical framework, it is possible to carry out a comprehensive analysis of the performance of stochastic heat engines based on Brownian particles subject to periodically time-dependent potentials and temperatures [17–21]. Along the same lines, optical micromanipulation techniques have facilitated during the last decade the experimental realization of simple colloidal heat engines, which are composed of a single colloidal particle as a working substance, embedded in water as a heat reservoir, undergoing thermodynamic cycles controlled by a harmonic optical potential [22–26]. In such colloidal systems, expansions and compressions during Stirling- and Carnot-like cycles are achieved by decreasing and increasing the trap stiffness, respectively, while a hot reservoir is realized either by an actual increase of the local temperature of the around the particle or by addition of synthetic noise of non-thermal origin. These experiments have paved the way for the investigation of stochastic models of colloidal heat engines in more intricate and realistic situations, such as passive Brownian engines operating in contact with active baths [27–32], Brownian engines with a self-propelled particle as working substance in contact with a viscous fluid [32–35] or in a suspension of passive Brownian particles [36] as a heat bath, as well as the realization of a colloidal Stirling engine in bacterial baths with tunable activity [37].

It must be pointed out that, in most of the situations envisaged for biological and technological applications, the fluid environment of a colloidal heat engine is not perfectly Newtonian with a constant viscosity, but possesses a complex viscoelastic microstructure because of the presence of macromolecules, e.g., biomolecular chains, polymers and wormlike micelles, or colloids suspended in a solvent, thus exhibiting time-dependent flow properties [38]. Therefore, the motion of a colloidal particle in such materials lacks a clear-cut separation from timescales of the surroundings, which results in memory effects with large relaxation times. All these features give rise to a wealth of intriguing transient effects that markedly manifest themselves when time-dependent driving forces are exerted on an embedded particle [39–44], and are absent in the case of purely viscous fluids. Although all these conditions are met by a colloidal heat engine operating in a complex fluid, to the best of our knowledge they have never been examined in the context of stochastic thermodynamic cycles. Therefore, it is of paramount importance to assess the role of viscoelasticity in the performance of this kind of engines, since the resulting frequency-dependent friction experienced by a colloidal particle can significantly impact the rate at which energy is dissipated into a viscoelastic bath [45–48].

Here, we investigate a model based on the generalized Langevin equation for the operation of a stochastic Stirling engine composed of a Brownian particle embedded in a viscoelastic fluid bath, which includes a memory kernel and colored noise to account for retarded friction effects and thermal fluctuations of the medium on the particle motion. By numerically solving the corresponding non-Markovian equation of motion, we analyze the effect of the characteristic relaxation time of the fluid on the performance of the engine under finite-time Stirling cycles, and compare our results with those found in the case of Brownian particle in a Markovian bath. We uncover a significant increase in the power output and the efficiency of the engine operating in a viscoelastic environment with respect to the corresponding values in a viscous bath at a given cycle time. Moreover, with increasing relaxation time of the fluid, the convergence to the quasi-static Stirling efficiency is shifted to monotonically decreasing values of the cycle period, thereby expanding the interval at which the engine is able to efficiently deliver positive power.

## 2. Model

We consider a stochastic heat engine consisting of a Brownian particle embedded in a viscoelastic fluid as a heat bath, whose motion is confined by a harmonic potential. Both the curvature of the confining potential and the temperature of the system can be varied in time according to a well-specified periodic protocol that mimics a macroscopic thermodynamic cycle. Therefore, a stochastic model of the particle dynamics that allows for temporal variations of the temperature is needed. Based on Zwanzig's pioneering work [49], Brey et al. [50] and Romero-Salazar et al. [51] derived the simplest equations of motion of a Brownian particle coupled to a heat bath with temperature changing in time. Their approach incorporates linear dissipative terms in the equations of motion of the surrounding bath particles, which account for continuous cooling or heating of the system controlled by some external mechanism in such a way that the bath particles are always in a canonical equilibrium at a well-behaved temperature dependent on time. In particular, in one dimension the generalized Langevin equation for the position *x*(*t*) at time *t* > 0 of the Brownian particle subject to a potential *U*(*x*(*t*), *t*), reads [50, 51].

where *m* is the mass of the particle, *T*(*s*) is the temperature of the system at time 0 ≤ *s* ≤ *t*, and *K*(*t* − *s*) is a memory kernel that weights the effect of the previous history of the particle motion at time *s* on its current drag force at time *t* due to the temporal correlations induced by the surrounding medium. In addition, in Equation (1), ζ(*t*) is a Gaussian stochastic force which accounts for thermal fluctuations in the system and satisfies

Extensions of Equation (1) to the three dimensional case, *r* = (*x, y, z*), which are relevant in many experimental situations using optical trapping techniques [52], are possible by a proper choice of the potential *U*(*r, t*) and a tensorial form of the memory kernel for particles of arbitrary shape [53]. Here, for the sake of simplicity we focus on the dynamics of a single coordinate of a spherical particle of radius *a*, which is confined by a harmonic potential $U(x,t)=\frac{1}{2}\kappa (t)x{(t)}^{2}$, where κ(*t*) is the stiffness at time *t* of the corresponding restoring force. Moreover, we assume that the fluid bath is incompressible and the time-dependent variation of κ(*t*) and *T*(*t*) are such that its rheological properties remain in the linear viscoelastic regime, which is completely characterized by the stress relaxation modulus *G*(*t*), or equivalently, by the complex dynamic shear modulus at frequency ω > 0, *G*^{*}(ω) = *iωη*^{*}(ω), where $i=\sqrt{-1}$ and η^{*}(ω) is the complex viscosity given by the Fourier transform of *G*(*t*), i.e., ${\eta}^{*}(\omega )={\int}_{-\infty}^{\infty}dt{e}^{-\text{i}\omega t}G(t)$ [54]. In general, *G*(*t*) is a function that decays to zero over a finite time-scale whose value is many orders of magnitude greater than those of simple viscous fluids [38]. For *a* larger than the characteristic length-scales of the fluid microstructure, the Fourier transform of the memory kernel, $\widehat{K}(\omega )={\int}_{-\infty}^{\infty}dt{e}^{-\text{i}\omega t}K(t)$, is related to η^{*}(ω) by the generalized Stokes relation [55, 56].

with ρ the density of the fluid. Furthermore, when *a* is much smaller than the so-called viscoelastic penetration depth, $\sqrt{\frac{\left|{\eta}^{*}(\omega )\right|}{\rho \omega}}$, as typically occurs for micron-sized particles suspended in most viscoelastic fluids, inertial flow effects are negligible [57]. In such a case, Equation (3) can be approximated to $\widehat{K}(\omega )=6\pi a{\eta}^{*}(\omega )$ [58], which yields the simple relation *K*(*t*) = 6π*aG*(*t*) by Fourier inversion. This leads to the following Langevin equation for the position of the Brownian heat engine in the overdamped limit

In the following, we focus on a fluid relaxation modulus consisting of a Dirac delta function plus an exponential decay

which models the rheological response of several viscoelastic fluids, such as wormlike micelles [42, 59, 60], some polymer solutions [61, 62], and to a great extent, the linear viscoelasticity over certain time intervals of intracellular fluids [63, 64], block copolymers [65], and λ-phage DNA [42, 66], where τ_{0} is the relaxation time of their elastic microstructure, whereas η_{0} and η_{∞} represent the zero-shear viscosity and the background solvent viscosity, respectively. Therefore, the corresponding friction memory kernel is

where the complex conjugate of its Fourier transform, ${\widehat{K}}^{*}(\omega )=6\pi a\eta (\omega )$, represents a frequency-dependent friction

In Equations (6) and (7), γ_{∞} = 6π*rη*_{∞} and γ_{0} = 6π*rη*_{0} ≥ γ_{∞} are friction coefficients characterizing dissipation at short and long timescales, respectively, whereas elastic effects are quantified by (γ_{0} − γ_{∞})τ_{0}. Hence, in this case, Equation (4) takes the form

It is noteworthy that, at constant temperature *T* and in absence of a trapping potential, the mean square displacement of a particle whose motion is described by Equation (8), is

which implies that in the long-time limit, *t* ≫ γ_{∞}τ_{0}/γ_{0}, it would perform free diffusion like in a Newtonian fluid with constant viscosity η_{0} [67–69], i.e., $\langle \Delta x{(t)}^{2}\rangle \approx \frac{2{k}_{B}T}{{\gamma}_{0}}t$. This provides a clear criterion for a direct comparison of the performance of a Brownian engine in a viscoelastic fluid bath with that in a viscous medium of the same zero-shear viscosity, i.e., η = η_{0}, under identical time-dependent variations of κ(*t*) and *T*(*t*). Furthermore, we introduce the dimensionless parameter

in such a way that, for either α = 0 or τ_{0} → 0, the memory kernel becomes *K*(*t*) = 2γδ(*t*), with constant friction coefficient γ = γ_{0} = γ_{∞}. Consequently, in these cases Equation (4) reduces to

where the thermal noise ζ(*t*) simply satisfies [50].

Equation (11) describes the motion of a Brownian particle coupled to a viscous heat bath with time dependent temperature *T*(*t*) through the frictional force $-\gamma \frac{dx(t)}{dt}$ and the thermal stochastic force, subject to a restoring force −κ(*t*)*x*(*t*). It should be noted that this situation was explicitly considered in many of the models of single-particle heat engines reported in the literature [18, 20, 27, 28, 31–35].

We point out that the rheological properties of viscoelastic fluids are generally dependent on their temperature, which under a thermodynamic cycle would also become time-dependent. The inclusion of such thermal effects in the minimal Langevin model (8) is not trivial and even a phenomenological description through additional rheological parameters and time-scales would render it little useful for a clear interpretation of the memory effects of a frequency-dependent friction in the performance of the Brownian engine. Therefore, similar to the simplifications made in most single-particle models of heat engines working in purely viscous fluids, as a first approximation we assume that η_{0}, η_{∞} and τ_{0} remain constant over time. The effect of the temperature dependence of these parameters is out of the scope of the present paper and will be the subject of further work.

The operation of the Brownian engine during a Stirling cycle of duration τ is depicted in Figure 1A, where the trap stiffness and the temperature are varied in time *t* according to the following protocols

and

respectively, where δκ = κ_{M} − κ_{m} > 0 and *T*_{h} > *T*_{c}. More specifically, a full cycle consists of a sequence of four steps:

1 → 2: For 0 ≤ *t* < τ/2, the colloidal engine undergoes an isothermal expansion at high temperature *T*_{h} by linearly decreasing the trap stiffness from κ_{M} to κ_{m}.

2 → 3:At *t* = τ/2, the temperature is suddenly decreased to *T*_{c}, while keeping the trap stiffness at κ(*t* = τ/2) = κ_{m}, thus corresponding to a isochoric-like process.

3 → 4:For τ/2 < *t* < τ, the engine undergoes an isothermal compression at low temperature *T*_{c} by linearly increasing the trap stiffness from κ_{m} to κ_{M}.

4 → 1:At *t* = τ, the temperature is suddenly raised to *T*_{h}, while keeping the trap stiffness at κ(*t* = τ) = κ_{M}, i.e., an isochoric-like process, thus completing the full cycle.

**Figure 1. (A)** Schematic representation of a Stirling cycle of period τ performed by a colloidal heat engine, embedded in a viscoelastic fluid, by means of the temporal variation of the stiffness κ(*t*) of a trapping harmonic potential $U(x,t)=\frac{1}{2}\kappa (t)x{(t)}^{2}$ and of the bath temperature *T*(*t*). At time 0 ≤ *t* < τ/2, the trap stiffness is linearly decreased from κ_{M} to κ_{m} < κ_{M} while keeping the temperature of the surroundings at high temperature *T*_{h} (step 1 → 2). At *t* = τ/2, the temperature is suddenly decreased to *T*_{c} < *T*_{h} (step 2 → 3), and kept at that value for τ/2 < *t* < τ, while linearly increasing the trap stiffness from κ_{m} to κ_{M} (step 3 → 4). The cycle is completed at *t* = τ, at which the temperature is again increased to *T*_{h} (step 4 → 1). For a given realization of the cycle, the particle position *x*(*t*) encodes the information of the stochastic energy exchange between the particle and the surrounding fluid, as depicted by the noisy trajectory obtained by numerical simulations of Equation (8). **(B)** Schematic representation of the Brownian Stirling cycle in a 〈*x*^{2}〉-κ^{−1} diagram, similar to the pressure-volume diagram of a gas.

Then, the cycle is repeated until the system reaches a time-periodic steady state, which becomes independent of the choice of the initial condition *x*(*t* = 0) = *x*_{0}. Note that, by analogy with a macroscopic Stirling cycle of a gas as a working substance, here the inverse of the trap stiffness and the variance of the particle position play the role of the volume and pressure, respectively, as depicted in Figure 1B.

According to stochastic thermodynamics [14], the work done on the system by the time variation of the optical trap over a single stochastic realization of the (*n*+1)−th cycle starting at *t*_{n} = *nτ*, with *n* = 0, 1, 2, …, is

whereas the heat dissipated into the bath during the first half period of the cycle is given by

In Equation (16), *W*_{τ/2} is the work done during the first half of the cycle, and Δ*U*_{τ/2} is the corresponding variation of the potential energy in the harmonic trap, *U*(*x, t*), in accordance with the stochastic extension of the first law of thermodynamics. Positive and negative values of *W*_{τ} correspond to work done on the particle and work performed by the particle, respectively, whereas positive and negative values of *Q*_{τ/2} represent heat transferred from the particle to the bath and heat absorbed by the particle, respectively. It must be noted that the mean steady-state values of the two stochastic variables given by Equations (15) and (16), which will be denoted as 〈*W*_{τ}〉 and 〈*Q*_{τ/2}〉, respectively, are the ones needed for the calculation of the efficiency of the Stirling heat engine [17]. They involve the variance of the particle position at an arbitrary time *t* ≥ 0, 〈*x*(*t*)^{2}〉, with *t* = 0 the time defining the initial condition, computed over an ensemble of independent realizations of the colored noise ζ(*t*) defined by Equations (2). An analytical treatment of this problem requires the explicit solution of the generalized Langevin Equation (8), which is not trivial even in the simpler case of a constant trap stiffness and constant temperature [48]. Therefore, to address the problem of the performance of a Brownian Stirling heat engine described by Equations (2), (8), (13) and (14), we opt for numerical simulations of the corresponding stochastic dynamics.

### 2.1. Numerical Solution

In order to compute the probability distributions of the work and the heat defined in Equations (15) and (16), as well as their corresponding mean values, the non-Markovian Langevin Equation (8) must be numerically solved. To this end, we express it in an equivalent Markovian form by introducing an auxiliary stochastic variable, *z*(*t*), defined as

where

represents a diffusion coefficient associated to the effective friction γ_{0} − γ_{∞}, which depends on the instantaneous value of the temperature at time *s*, *T*(*s*), and ξ_{z}(*s*) is a Gaussian noise satisfying

Consequently, the non-Markovian Langevin Equation (4) for *x*(*t*) can be written as a linear system of two coupled Markovian Langevin equations

with α defined in Equation (10). In Equation (20), *D*_{∞}(*t*) is a short-time diffusion coefficient associated to the infinite-frequency friction coefficient ${lim}_{\omega \to \infty}{\widehat{K}}^{*}(\omega )={\gamma}_{\infty}$, see Equation (7), at temperature *T*(*t*), and is given by

whereas ξ_{x}(*t*) is a Gaussian noise which satisfies

Note that, apart from the step-like changes at *t* = *t*_{n} and $t={t}_{n}+\frac{\tau}{2}$, *T*(*t*) remains constant. Accordingly, the rate of change of the time-dependent temperature in Equation (20) vanishes during each half a Stirling cycle, i.e., $\frac{d}{dt}\left[ln\text{}T(t)\right]=0$.

To compute the probability distributions of *W*_{τ} and *Q*_{τ/2}, we carry out numerical simulations of the stochastic process [*x*(*t*), *z*(*t*)] starting from the initial condition [*x*(*t* = 0) = 0, *z*(*t* = 0) = 0] with a total length of 2 × 10^{4} times the period τ. To ensure that the system is always in a time-periodic non-equilibrium steady state independent of the choice of the initial condition, the first 10^{4} cycles are left out and the origin of time is shifted to the beginning of the (10^{4} + 1)−st cycle. Furthermore, without loss of generality we choose constant values of the low and high-frequencies viscosities that are typical of viscoelastic fluids prepared in aqueous solution in semidilute regimes [42, 62, 66, 70, 71]: η_{0} = 0.040 Pa s and η_{∞} = 0.004 Pa s, which correspond to α = 9. The diameter of the colloidal particle is set to *a* = 0.5μm, while the maximum and minimum values of the trap stiffness during the Stirling cycle are chosen as ${\kappa}_{M}=5\text{pN}\mu {\text{m}}^{-1}$ and ${\kappa}_{m}=1\text{pN}\mu {\text{m}}^{-1}$, respectively, which are easily accessible with optical tweezers [52]. The temperatures of the reservoir during the hot and cold part of the cycle are *T*_{c} = 5°C and *T*_{h} = 90°C, which are selected in such a way that they are within the temperature range in which water, which is a common solvent component of many viscoelastic fluids, remains liquid. On the other hand, to study the influence of the fluid relaxation time on the performance of the colloidal Stirling engine, τ_{0} is varied in the range of 0.01–100 s, which also covers characteristic values in actual experimental systems. We solve Equations (20) by means of an Euler–Cromer scheme with time step δ*t* = 10^{−4} s, which is about 75 times smaller than the shortest relaxation time of the system, γ_{∞}/κ_{M}. In the case of the Stirling heat engine in a Newtonian viscous fluid, we solve numerically Equation (11) with constant friction coefficient γ = 6π*aη*, where η = η_{0} = 0.040 Pa s and the rest of the involved parameters, namely, κ_{m}, κ_{M}, *a*, *T*_{c}, *T*_{h}, and δ*t*, are selected with the same values as described before for the viscoelastic case for a direct comparison between both systems. We also explore different values of the cycle period, 0.01 s ≤ τ ≤ 50 s, which allows us to examine the approach of the computed quantities to the quasi-static values τ → ∞. We note that τ_{κ} ≡ γ_{0}/κ_{m} represents the slowest dissipation time-scale of the system [26], and appears explicitly in the analytical expressions for the variance of a Brownian particle undergoing a finite-time Stirling cycle in contact with a viscous heat bath [34]. Therefore, in both cases of the viscous and viscoelastic baths analyzed here, all the timescales are normalized by τ_{κ}, whereas energies are normalized by *k*_{B}*T*_{c}.

## 3. Results and Discussion

Since *W*_{τ} and *Q*_{τ/2} are stochastic variables, we first present the results for their probability distributions, ϱ(*W*_{τ}) and ϱ(*Q*_{τ/2}), respectively, for different values of the time-scales τ and τ_{0}. In Figures 2A,B, we plot such distributions for a value of the fluid relaxation time that is comparable to the largest dissipation time-scale of the system: τ_{0} = 2.65τ_{κ}, at which memory effects due to the frequency-dependent friction must be important. In such a case, we observe that for fast Stirling cycles with period τ smaller or comparable to τ_{κ} the work distribution is asymmetric with respect to its maximum and exhibits pronounced exponential tails, as illustrated in the inset of Figure 2A. In addition, large positive work fluctuations occur for small τ, which indicates the existence of rare events where work is done on the particle during a cycle, thus effectively consuming energy as a heat pump. As τ increases, the exponential tails and their asymmetry vanish, thus giving rise to a narrower Gaussian-like shape for τ ≫ τ_{k}. This shows that the probability of finding positive work fluctuations decreases by increasing τ, i.e., the Brownian particle behaves more and more like a macroscopic Stirling engine, which on average is able to convert the heat absorbed from the viscoelastic bath into work. On the contrary, the heat distribution does not significantly change with the cycle time τ, as shown in Figure 2B. In this case, clear exponential tails remain even for large values of τ, as revealed in the inset of Figure 2B. where the probability of occurrence of negative heat fluctuations is higher than that of positive ones. Hence, regardless of the cycle period τ, it is more likely that heat is absorbed by the particle than dissipated into the bath during the isothermal expansion at temperature *T*_{h}.

**Figure 2. (A)** Probability density function of the work *W*_{τ}, and **(B)** the heat *Q*_{τ/2}, for a Brownian Stirling engine in contact with a viscoelastic fluid bath with relaxation time τ_{0} = 2.65τ_{κ}, for different values of the cycle time τ. **(C)** Probability density function of the work *W*_{τ}, and **(D)** the heat *Q*_{τ/2}, for a Brownian Stirling engine during a cycle of duration τ = τ_{κ}, in contact with viscoelastic fluid baths with the same zero-shear viscosity η_{0} = 0.040 Pa a, and distinct relaxation times τ_{0} spanning 5 orders of magnitude (solid lines). **(E)** Probability density function of the work *W*_{τ}, and **(F)** the heat *Q*_{τ/2}, for a Brownian Stirling engine during a cycle of duration τ = 10τ_{κ}, in contact with viscoelastic fluid baths with the same zero-shear viscosity η_{0} = 0.040 Pa s and distinct relaxation times τ_{0} spanning 5 orders of magnitude (solid lines). In **(C–F)**, the dotted lines represents the corresponding curves for a Brownian engine in a Newtonian fluid (τ_{0} = 0) with constant viscosity η_{0} = 0.040 Pa s. The insets are semilogarithmic representations of the main plots.

In Figures 2C,D, we analyze the dependence on the fluid relaxation time τ_{0} of the work and heat distributions, respectively, for Stirling cycles of period τ = τ_{κ}, i.e., similar to the largest viscous dissipation time-scale of the system. For comparison, we also plot as dotted lines the corresponding probability distributions for a colloidal engine in a fluid with constant viscosity η = η_{0}, for which τ_{0} = 0. Remarkably, we find that the fluid viscoelasticity, through the parameter τ_{0}, has a strong influence on the resulting shape of the distributions. For a viscous bath, the work has large exponential tails with a highly asymmetric shape. A similar shape is observed for a viscoelastic bath at sufficiently small τ_{0}, but the width and the asymmetry of the distribution gradually decrease as τ_{0} increases, then converging to a single limiting curve with a rather symmetric profile for sufficiently large values of the fluid relaxation time τ_{0} ≳ τ_{κ}, as shown in the inset of Figure 2C. In addition, the heat distribution also has exponential tails with a width that does not strongly depend on the fluid relaxation time τ_{0}, but the location of the maximum is slightly shifted to more and more negative values of *Q*_{τ/2} with increasing τ_{0}, as shown in Figure 2D. Finally, for values of the cycle duration τ larger than τ_{κ}, the shape of ϱ(*W*_{τ}) changes from a rather symmetric exponentially-tailed distribution to a limiting Gaussian curve with increasing τ_{0}, whereas ϱ(*Q*_{τ/2}) exhibits a symmetric profile with exponential tails peaked at a negative value of *Q*_{τ/2}, which remains unaffected by the τ_{0}, as respectively shown in Figures 2E,F for τ = 10τ_{κ}. It is important to realize that for τ > τ_{κ}, the work distribution of the Brownian engine is narrower in a viscoelastic bath as compared to that in a viscous bath with the same zero-shear viscosity. This can be attributed the elastic response in the former case, which prevents large instantaneous heat losses into the bath by viscous dissipation, thus resulting in a more efficient conversion into work of the energy extracted from the surroundings. This observation underlines the importance of the friction memory kernel of the particle motion in the viscoelastic fluid, which becomes strongly dependent on the frequency imposed by the Stirling cycle. Thus, for sufficiently small τ_{0} < τ_{κ} the energy exchanges between the Brownian particle and the viscoelastic bath must not be that different from those occurring in a viscous fluid, while for sufficiently large τ_{0} > τ_{κ} significant deviations must take place, as verified in Figures 2C,E for τ_{0} = 0.0265τ_{κ} and τ_{0} = 265τ_{κ}, respectively.

To investigate the performance of a Brownian engine operating in a viscoelastic bath, in Figure 3A we plot the mean work *done by the Brownian engine* during a cycle, i.e., −〈*W*_{τ}〉. In a Newtonian fluid, −〈*W*_{τ}〉 is positive at sufficiently large τ and monotonically saturates to a constant positive value in the quasi-static limit τ → ∞ [34].

whereas it becomes negative at small values of τ and tends to zero as τ → 0 according to Equation (15), thus implying that it has a minimum at a certain value of τ. This is verified in Figure 3A, where we plot as dotted and dashed lines the curves corresponding to the work done by a particle in viscous fluids with constant viscosities η = η_{0} = 0.040 Pa s and η = η_{∞} = 0.004 Pa s, respectively, i.e., equal to the viscosities characterizing the long-time and short-time dissipation of the viscoelastic fluid. The location of the minimum, which is depicted by arrows, depends on the specific value of η, but the general shape of the curve in a linear-logarithmic representation is the same, as observed in Figure 3A. Interestingly, in the case of viscoelastic fluids with non-zero values of τ_{0}, the work done by the particle exhibits an intermediate behavior between these two curves. For instance, for τ_{0} = 0.0265τ_{κ} ≪ τ_{κ}, the dependence of −〈*W*_{τ}〉 on τ is very similar to that in a Newtonian fluid with viscosity η = η_{0}, with a single minimum at the same location (τ ≈ τ_{κ}) and only small deviations of the respective values along the vertical axis. Nevertheless, as τ_{0} increases, a second local minimum emerges at τ ≈ 0.1τ_{κ}, i.e., at the location of the minimum of the curve corresponding to the Newtonian fluid of viscosity η = η_{∞}, as observed in Figure 3A for τ_{0} = 0.084τ_{κ}. Such a second minimum becomes more and more apparent with increasing τ_{0}, whereas the first minimum at τ ≈ τ_{κ} becomes less and less dominant, as seen for τ_{0} ≥ 0.265τ_{κ}. Unexpectedly, for τ_{0} ≫ τ_{κ}, the curves for the viscoelastic case converge to that for a Newtonian fluid with a viscosity η = η_{∞}. These observations suggest that, depending of the specific values of the fluid relaxation time and the cycle time with respect to τ_{κ}, different dissipation mechanisms take place in order for the particle to convert the energy taken from the bath into work by means of the applied thermodynamic cycle.

**Figure 3. (A)** Mean work done by the Brownian Stirling engine during a cycle, −〈*W*_{τ}〉, as a function of the cycle time τ, for different values of the fluid relaxation time τ_{0} (solid lines). The dotted and dashed lines represent the mean power output of a Brownian Stirling engine in Newtonian fluids with constant viscosities η = η_{0} and $\eta ={(1+\alpha )}^{-1}{\eta}_{0}={\eta}_{\infty}$, respectively. The arrows depict the location of the corresponding minima. **(B)** Mean power output per cycle of the colloidal Stirling engine, *P*_{τ}, as a function of the cycle duration τ, for different values of the fluid relaxation time τ_{0} (solid lines). The dotted and dashed lines represent the mean power output of a Brownian Stirling engine in Newtonian fluids with constant viscosities η = η_{0} and $\eta ={(1+\alpha )}^{-1}{\eta}_{0}={\eta}_{\infty}$, respectively. Same color code as in **(A)**. Inset: stall time of the engine defined in Equation (26), τ^{*}, as a function of the fluid relaxation time, τ_{0} (thick line). The thin line represents the value of τ^{*} for an engine operating in a Newtonian fluid with constant viscosity η = η_{∞}.

Next, we compute the mean power produced by the engine during a cycle

whose dependence on the cycle time τ is plotted as solid lines in Figure 3B for some exemplary values of the fluid relaxation time τ_{0}. Besides, we also plot in Figure 3B as a dotted line the mean power for a Brownian engine in a Newtonian fluid bath with viscosity η = η_{0}. It is important to note that, for all values of τ_{0}, *P*_{τ} exhibits a non-monotonic behavior as a function of τ, which gradually deviates from the behavior in a Newtonian fluid with viscosity η = η_{0} as τ_{0} increases. This is the result of the pronounced non-monotonic dependence of −〈*W*_{τ}〉 on τ shown in Figure 3A. In particular, *P*_{τ} has a maximum that originates from the trade-off between high energy dissipation at small cycle times τ (high frequency operation) and large τ (slow operation), at which net work is produced by the engine with low dissipation. Additionally, the general shape of all power curves displays three different operation regimes. For sufficiently slow Stirling cycles (large τ), the engine is able to deliver net power on average (*P*_{τ} > 0), where the irreversible energy dissipation into the bath becomes negligible. On the other hand, there is a specific value of the cycle time at which the engine stalls, i.e., both the mean work and the power output vanish: 〈*W*_{τ}〉 = 0, *P*_{τ} = 0 [17]. Finally, for sufficiently fast cycles (small τ), the engine absorbs energy (*P*_{τ} < 0) rather than delivering it, thus behaving like a heat pump. This regime is the consequence of the large amount of energy irreversibly dissipated when the particle is quickly driven by the periodic variation of κ(*t*) and *T*(*t*). Interestingly, in Figure 3B, we show that the value of the fluid relaxation time τ_{0} has a considerable impact on the mean power output, and in particular, on the value of the cycle time at which the Brownian engine stalls, which we denote as τ^{*}

For instance, in the case of the Newtonian fluid (τ_{0} = 0) with η = η_{0}, we find ${\tau}^{*}=6.15{\tau}_{\kappa}$, while for a viscoelastic fluid (τ_{0} > 0), τ^{*} is smaller and decreases with increasing τ_{0}. In the inset of Figure 3B, we plot the dependence of τ^{*} on τ_{0}, where we can see that for sufficiently short fluid relaxation times, the stall time is close to that for a Newtonian fluid bath (${\tau}^{*}=6.15{\tau}_{\kappa}$), and monotonically decreases with increasing τ_{0}. In this short-τ_{0} regime, the performance of the engine is very sensitive to the specific value of τ_{0}, as shown by the strong variation of the shape of the power curves plotted in Figure 3B for τ_{0} = 0.0265τ_{κ}, 0.084τ_{κ}, 0.265τ_{κ}, 0.84τ_{κ}. Around τ_{0} = τ_{κ}, a conspicuous change in the dependence on τ_{0} of the operation of the engine happens. Indeed, as τ_{0} increases the stall time converges to the constant value ${\tau}^{*}=0.52{\tau}_{\kappa}$, as verified in the inset of Figure 3A for τ_{0} > τ_{κ}. The monotonic decrease of τ^{*} implies that the interval of cycle times at which the engine is able to efficiently deliver positive power output is expanded with increasingly larger τ_{0}. Moreover, with increasing fluid relaxation times τ_{0} > τ_{κ}, which is consistent with increasingly pronounced viscoelastic behavior of the bath, the power output curves converge to a limiting curve, as shown in Figure 3B for τ_{0} = 2.65τ_{κ}, 8.4τ_{κ}, 26.5τ_{κ}, 84τ_{κ}, 265τ_{κ}. Remarkably, we find that such a limiting curve corresponds to the power curve of a Brownian Stirling engine in a Newtonian bath with viscosity equal to high frequency value η = η_{∞} = 0.004 Pa s, i.e., the viscosity of the solvent component in the viscoelastic fluid, which is represented as a dashed line in Figure 3B. As a consequence, the limit of the stall time of an engine working in a viscoelastic fluid with increasing τ_{0} corresponds to the stall time of a Brownian engine operating in a Newtonian one with constant viscosity η = η_{∞}, ${\tau}^{*}=0.52{\tau}_{\kappa}$, as verified in the inset of Figure 3B, see the horizontal solid line. Furthermore, in Figure 3B, we check that, for a given cycle of finite duration τ, the mean power output of the engine operating in a viscoelastic fluid is enhanced with increasing values of τ_{0} with respect to the power output in a Newtonian fluid of the same zero-shear viscosity. We also find that the location of the global maximum of each power output is shifted to smaller and smaller values of τ with increasing τ_{0}, whereas the value of *P*_{τ} at the maximum increases with increasing τ_{0} because of the decreasing irreversible dissipation taking place in a fluid with pronounced viscoelastic behavior.

These findings allows us to uncover the underlying mechanism behind the influence of fluid viscoelasticity on the performance of the engine. In a Newtonian fluid with constant viscosity η_{0}, the largest time-scale associated to viscous dissipation due to temporal changes in the trap stiffness is precisely τ_{κ}, which is proportional to η_{0}, and represents the largest relaxation time in the system. In this case, the viscous bath simply acts as a mechanically inert element of the engine which equilibrates instantaneously in response to the particle motion under the variations of the trap stiffness. On the other hand, when the bath is a viscoelastic fluid, the hidden degrees of freedom of its elastic microstructure, e.g., entangled micelles, polymers, interacting colloids, etc., also come into play in the dynamics and mechanically respond within a characteristic time τ_{0} > 0 to the temporal changes periodically imposed on the particle. Therefore, the interplay between τ_{κ} and τ_{0} determines the resulting energetic behavior of the system:

• If τ_{0} ≪ τ_{κ}, the fluid microstructure fully relaxes before the energy dissipation into the bath takes place on a time-scale τ_{κ}. In such circumstances, the Brownian particle has enough time to probe the long-time (low frequency) properties of the fluid environment with friction coefficient ${\widehat{K}}^{*}(\omega \to 0)={\gamma}_{0}=6\pi a{\eta}_{0}$, see Equation (7), thereby leading to a stochastic energetic behavior similar to that in a Newtonian fluid with constant viscosity η = η_{0}.

• If τ_{0} ≲ τ_{κ}, excessive irreversible energy losses by viscous dissipation are counterbalanced by the transient energy storage in the elastic structure of the bath, because at frequencies $\omega ~{\tau}_{0}^{-1}$ the imaginary part of ${\widehat{K}}^{*}(\omega )$ is not negligible. Therefore, the value τ_{0} ≈ τ_{κ} marks a qualitative change in the energy exchange between the particle and bath.

• If τ_{0} > τ_{κ}, the elastic fluid microstructure does not have enough time to mechanically relax to the temporal changes of the cycle, thus preventing the particle from undergoing the long-time friction characterized by the coefficient γ_{0}. Therefore, the particle can only probe the short-time response of the surrounding fluid through the high-frequency components of the friction, which correspond to ${\widehat{K}}^{*}(\omega \to \infty )={\gamma}_{\infty}=6\pi a{\eta}_{\infty}$ for τ_{0} ≫ τ_{κ} according to Equation (7). As a consequence, in this limit the relevant dissipation timescale is γ_{∞}/κ_{m}, which is in general smaller than τ_{κ} because ${\eta}_{\infty}={\eta}_{0}{(1+\alpha )}^{-1}\le {\eta}_{0}$. For instance, for the numerical values chosen in the simulations presented here, γ_{∞}/κ_{m} = 0.1τ_{κ}. Accordingly, less irreversible dissipation must take place in the viscoelastic fluid under finite-time Stirling cycles, thus enhancing the net power output of the engine at a given cycle time τ as compared to that in a Newtonian fluid with the same zero-shear viscosity η_{0}.

To confirm the previously described mechanism of energy storage and dissipation during the Stirling cycle, in Figure 4A, we plot the mean heat absorbed by the particle during the hot step of the cycle, −〈*Q*_{τ/2}〉, as a function of the total duration τ of a full cycle. We find that, for all values of τ and of the fluid relaxation time τ_{0}, −〈*Q*_{τ/2}〉 ≥ 0, which means that the particle absorbs heat on average during the first half of the cycle. In particular, for a given τ_{0} the mean absorbed heat increases monotonically from the value −〈*Q*_{τ=0}〉 = 0, and saturates to a constant value corresponding to a quasi-static process as τ → ∞. For comparison, in Figure 4A, we also plot as a dotted line the mean heat absorbed by the Brownian engine when operating in a Newtonian fluid with viscosity η = η_{0}. In such a case, it can be readily demonstrated from Equation (16) that −〈*Q*_{τ/2}〉 actually approaches a quasi-static value, which is explicitly given by [34]

For the numerical values of the parameters investigated here, −〈*Q*_{(τ → ∞)/2}〉 = 1.203*k*_{B}*T*_{c}, see horizontal thin dotted line in Figure 4A. We observe that, regardless of τ_{0}, all heat curves converge to such a value for τ ≫ τ_{κ}, but depending on the specific value of the fluid relaxation time, different behaviors occur at short and intermediate cycle durations. Once again, we find that with increasing τ_{0}, the mean-heat curves gradually deviate from the behavior in a Newtonian fluid with viscosity η = η_{0}, and for τ_{0} ≫ τ_{κ} they converge to that in a Newtonian fluid with η = η_{∞}, see dashed line in Figure 4A. This provides another evidence that, as τ_{0} increases, the energy dissipation of an engine operating in a viscoelastic fluid is mainly determined by the friction with the solvent.

**Figure 4. (A)** Mean heat absorbed by the Brownian Stirling engine during the isothermal expansion at high temperature, −〈*Q*_{τ/2}〉, as a function of the cycle time τ, for different values of the fluid relaxation time τ_{0} (solid lines). The dotted and dashed lines corresponds to the mean heat absorbed by the colloidal engine in Newtonian fluids with constant viscosities η = η_{0} and $\eta ={(1+\alpha )}^{-1}{\eta}_{0}={\eta}_{\infty}$, respectively. The horizontal thin dotted line represents the quasi-static value given by Equation (27). **(B)** Mean rate of heat absorption by the colloidal engine during the isothermal expansion at high temperature, *J*_{τ}, as a function of the duration of the cycle τ, for different values of the fluid relaxation time τ_{0} (solid lines). Same color code as in **(A)**. The dotted and dashed lines correspond to mean rate of heat absorption in Newtonian fluids with constant viscosities η = η_{0} and η = η_{∞}, respectively. The arrows depict the location of the corresponding maxima. The dotted-dashed line depicts the behavior ~ τ^{−1}.

In Figure 4B, we plot as solids lines the mean rate of heat absorption by the engine from the bath during the isothermal expansion at temperature *T*_{h}

as a function of the cycle duration τ for some representative values of τ_{0} > 0. The corresponding curves for a particle in Newtonian fluids with η = η_{0} and η = η_{∞} are represented as dotted and dashed lines, respectively. In such cases, we find that *J*_{τ} exhibits a maximum, which corresponds approximately to the location of the minima in −〈*W*_{τ}〉 shown in Figure 3A. For τ ≲ τ_{κ}, a marked dependence on the fluid relaxation time is observed if τ_{0} ≲ τ_{κ}, while for τ ≫ τ_{κ} a dependence ${J}_{\tau}~{\tau}^{-1}$ on the Stirling cycle time emerges for all values of τ_{0}, thus indicating the onset of the quasi-static thermodynamic behavior.

The previous findings reveal that, unlike the performance of Brownian heat engines in a Newtonian environment with a single relevant time-scale γ_{0}/κ_{m} of energy dissipation, in a viscoelastic fluid bath the low-frequency and the high-frequency values of the friction, γ_{0} and γ_{∞}, give rise two meaningful dissipation time-scales, namely τ_{κ} = γ_{0}/κ_{m} and the apparently hidden time-scale ${(1+\alpha )}^{-1}{\tau}_{\kappa}={\gamma}_{\infty}/{\kappa}_{m}$ due to the friction of the particle with the solvent. When the Stirling cycle time τ is comparable to one of such time-scales, the corresponding channel of irreversible dissipation is strongly activated. This in turn leads to a large amount of energy absorbed by the particle from the heat bath at a very high rate, as manifested by the minima and maxima depicted by arrows in Figures 3A, 4B, respectively. For an arbitrary cycle time, the interplay between the two channels of irreversible dissipation along with the transient energy storage by the elastic microstructure of the fluid determine the resulting performance of the Brownian engine.

Finally, we determine the efficiency of the Brownian Stirling engine, defined as

as a function of the cycle time, τ, and the relaxation time of the viscoelastic fluid, τ_{0}. The results are represented as a 2D color map in Figure 5A, with some efficiency curves plotted in Figure 5B as a function of τ for exemplary values of τ_{0}. Additionally, in Figure 5A, we also plot the stall time τ^{*} defined in Equation (26) as a function of the fluid relaxation time. As a consequence of the energy exchange with a viscoelastic bath discussed in the previous paragraphs, τ^{*} divides the efficiency diagram into two regions. For τ < τ^{*} the Brownian particle behaves a heat pump, where ϵ_{τ} < 0 exhibits a rather intricate dependence on τ_{0} and τ due to the competition between the different energy storage and dissipation channels of the bath, which results on average in net energy absorption from the bath. On the other hand, for τ > τ^{*}, the efficiency is positive, ϵ_{τ} > 0, i.e., the Brownian particle behaves as a heat engine with positive power output. In this case, the efficiency is a monotonic increasing function of both τ and τ_{0}. Note that for the investigated values of the cycle time τ, the interval at which the engine has a positive efficiency is rather narrow for small fluid relaxation times τ_{0} < τ_{κ}, because for the large value of the zero-shear viscosity considered in the simulations (η_{0} = 0.040 Pa s, typical of biological fluids), there is a large amount of heat dissipation even at comparatively slow Stirling cycles. However, when the value of τ_{0} is similar or larger than τ_{κ}, the elastic response of the fluid takes effect, hence the decrease in energy dissipation with a subsequent broadening of the interval of cycle times by one order magnitude for which ϵ_{τ} > 0.

**Figure 5. (A)** 2D color map representation of the efficiency of the colloidal Stirling engine, ϵ_{τ}, as a function of the fluid relaxation time, τ_{0}, and the duration of a Stirling cycle, τ. The dotted-dashed line corresponds to the stall time τ^{*}, which separates the values of the parameters τ_{0} and τ for which the system operates as a heat engine (ϵ_{τ} > 0) from those for which it behaves as a heat pump (ϵ_{τ} < 0). **(B)** Examples of efficiency curves as function of the duration of a Stirling cycle, τ, for some particular values of the fluid relaxation time (solid lines). The dotted and dashed lines represent the efficiency curves of a Brownian Stirling engine in Newtonian fluids with constant viscosity η = η_{0} and η = η_{∞}, respectively. The dotted-dashed line depicts the quasi-static value of the Stirling efficiency, ϵ_{τ → ∞}, given by Equation (30).

Because in the model (4), we assume that the only source of stochasticity of the system is the thermal fluctuations of the fluid, apart from the driving potential of the harmonic trap there are no other sources of energy that affect the performance of the Brownian engine. Therefore, it is expected that the quasi-static Stirling efficiency

which can be determined from the ratio of Equations (24) and (27), is never exceeded at finite τ regardless of the relaxation time of the viscoelastic fluid. In Equation (30), ${\u03f5}_{C}=1-\frac{{T}_{c}}{{T}_{h}}$ corresponds to the efficiency of a Carnot engine operating quasi-statically between two reservoirs at temperatures *T*_{c} and *T*_{h}. For the numerical values of the parameters characterizing the Stirling cycle considered here, we find ϵ_{τ → ∞} = 0.2043. In Figure 5B, we demonstrate that, indeed, all the efficiency curves are bounded by such a value and approach it as the cycle time τ increases. The typical value of cycle period at which such efficiency is reached strongly depends on τ_{0}. While for a Brownian engine in a Newtonian fluid of viscosity η_{0} the convergence is very slow, the quasi-static Stirling efficiency can be reached in a viscoelastic bath for typical experimental values of the parameters of the system, as shown in Figure 5B for τ_{0} > τ_{κ}.

To compare the performance of a Stirling Brownian engine in a viscoelastic bath with other situations of practical interest, we first determine its efficiency at maximum power in a Newtonian fluid bath with the same zero-shear viscosity η = η_{0}. Although not as general as the Carnot efficiency, under some circumstances the so-called Curzon-Ahlborn efficiency [1, 2] represents a good approximation for the upper bound of the efficiency of stochastic heat engines working at maximum power [4, 6, 17, 19]. For the values of the parameters investigated in this work, we find that the power output *P*_{τ} in a purely viscous fluid reaches the maximum value ${P}_{{\tau}_{MP}}=0.00679{k}_{B}{T}_{c}{{\tau}_{\kappa}}^{-1}$ at a cycle time of τ_{MP} = 18.6τ_{κ}, at which the efficiency is η_{τMP} = 0.1218. This value compares well with the Curzon-Ahlborn efficiency, ${\u03f5}_{CA}=1-\sqrt{\frac{{T}_{c}}{{T}_{h}}}=0.1248$, and is approximately 60% the Carnot efficiency η_{C} = 0.2043. In Table 1, we list some exemplary values of the mean power output over a cycle, *P*_{τ =τMP}, and the corresponding efficiencies of the Brownian engine, ϵ_{τ =τMP}, operating at the same Stirling cycle time τ_{MP} = 18.6τ_{κ} in viscoelastic fluid baths with distinct values of their relaxation time τ_{0}. We verify that with increasing τ_{0}, both the absolute power delivered by engine and its efficiency are enhanced with respect to those in a Newtonian fluid. In particular, the efficiency at τ_{MP} = 18.6τ_{κ} converges to approximately 93% the Carnot efficiency for τ_{0} ≫ τ_{κ}.

**Table 1**. Mean power output produced by a Brownian Stirling engine during a cycle τ = τ_{MP} = 18.6τ_{κ} and corresponding efficiency for distinct values of the fluid relaxation time τ_{0}.

## 4. Summary and Final Remarks

In this work, we have investigated a stochastic model based on the generalized Langevin equation for a Brownian Stirling engine in contact with a viscoelastic fluid bath. The slow rheological behavior of the fluid is taken into account in the model by an exponentially decaying memory kernel, which captures the basic features of the linear viscoelastic behavior of many non-Newtonian fluids. Our findings demonstrate that the memory friction exerted by the surrounding fluid has a tremendous impact on the performance of the heat engine in comparison with its operation in a viscous environment with the same zero-shear viscosity. In particular, a pronounced enhancement of the power output and the efficiency of the engine occurs as a result of the frequency-dependent response of the fluid under finite-time Stirling cycles, thus converging to limiting curves determined by the high frequency component of the friction of the particle as the fluid relaxation time increases. Moreover, the minimum value of the duration of the Stirling cycle at which the Brownian engine can convert energy from the medium into work becomes monotonically shorter with increasing fluid relaxation time, which broadens the interval of possible values of the Stirling cycle duration over which the engine is able to efficiently deliver positive power. From a wider perspective, our results highlight the importance of the non-equilibrium transient nature of the particle friction under temporal cycles of finite duration. We point out that, although in a different context, qualitatively similar effects have been discussed in systems with frequency-dependent properties due to their coupling to non-Markovian baths, such as Brownian particles driven into periodic non-equilibrium steady states [72] and quantum Otto refrigerators [73]. Furthermore, the link between a frequency dependent friction and the noise correlations of the bath is in turn an important issue for the correct interpretation of the efficiency of stochastic heat engines operating in non-equilibrium baths, as recently examined in the case of underdamped active Brownian particles [74].

To the best of our knowledge, our work represents the first investigation on the effect of memory friction in the performance of a Brownian Stirling engine in contact with a viscoelastic fluid reservoir. Thus, we expect that the results presented in this paper will contribute to a better understanding and potential applications of efficient work extraction and heat dissipation in other types of mesoscopic engines operating in complex fluids. Further steps of our work aim at addressing long-term memory effects during stochastic thermodynamic cycles with finite period, as those described by stretched exponentials [75] and power law kernels and fractional Brownian noise [76–79], which describe the mechanical response of diverse soft matter systems such as glasses and biological materials [80, 81]. One further aspect that could be investigated in the future is the effect of temporal changes in the fluid parameters, as it is well-known that the rheological properties of viscoelastic fluids are dependent on their temperature, which under a thermodynamic cycle would become time-dependent. We would like to point out that, since the parameters characterizing the operation of the heat engine presented in this paper are representative of typical soft matter systems, we expect that this process can be realized in a straightforward manner by use of optical tweezers [52]. Similar ideas could be extended to Brownian particles in non-linear potentials [82], and active Brownian heat engines [32] functioning in complex fluids, which could be implemented by means of light-activated colloids in non-Newtonian liquids [69, 83–86] and hot Brownian particles [87–89].

## Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author Contributions

JG-S conceived the model, carried out the numerical simulations, analyzed the results, and wrote the manuscript.

## Funding

This work was supported by UNAM-PAPIIT IA103320.

## Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## References

1. Novikov II. The efficiency of atomic power stations (a review). *J Nuclear Energy*. (1958). 7:125–28. doi: 10.1016/0891-3919(58)90244-4

2. Curzon FL, Ahlborn B. Efficiency of a Carnot engine at maximum power output. *Am J Phys*. (1975) 43:22–4. doi: 10.1119/1.10023

3. Leff HS. Thermal efficiency at maximum work output: new results for old heat engines. *Am J Phys*. (1987) 55:602–10. doi: 10.1119/1.15071

4. Van den Broeck C. Thermodynamic efficiency at maximum power. *Phys Rev Lett*. (2005) 95:190602. doi: 10.1103/PhysRevLett.95.190602

5. Izumida Y, Okuda K. Molecular kinetic analysis of a finite-time Carnot cycle. *Europhys Lett*. (2008) 83:60003. doi: 10.1209/0295-5075/83/60003

6. Esposito M, Lindenberg K, Van den Broeck C. Universality of efficiency at maximum power. *Phys Rev Lett*. (2009) 102:130602. doi: 10.1103/PhysRevLett.102.130602

7. Ozin GA, Manners I, Fournier-Bidoz S, Arsenault A. Dream nanomachines. *Adv Mater*. (2005) 17:3011–8. doi: 10.1002/adma.200501767

8. Hänggi P, Marchesoni F. Artificial Brownian motors: controlling transport on the nanoscale. *Rev Mod Phys*. (2009) 81:387–442. doi: 10.1103/RevModPhys.81.387

9. Kim K, Guo J, Liang ZX, Zhu FQ, Fan DL. Man-made rotary nanomotors: a review of recent developments. *Nanoscale*. (2016). 8:10471–90. doi: 10.1039/C5NR08768F

10. Martinez IA, Roldan D, Dinis L, Rica RA. Colloidal heat engines: a review. *Soft Matter*. (2017). 13:22–36. doi: 10.1039/C6SM00923A

11. Pietzonka P, Fodor E, Lohrmann C, Cates ME, Seifert U. Autonomous engines driven by active matter: energetics and design principles. *Phys Rev X*. (2019) 9:041032. doi: 10.1103/PhysRevX.9.041032

12. Ciliberto S, Gomez-Solano R, Petrosyan A. Fluctuations, linear response, and currents in out-of-equilibrium systems. *Annu Rev Cond Matter Phys*. (2013) 4:235–61. doi: 10.1146/annurev-conmatphys-030212-184240

13. Sekimoto K. Langevin equation and thermodynamics. *Prog Theor Phys Suppl*. (1998) 130:17–27. doi: 10.1143/PTPS.130.17

14. Seifert U. Stochastic thermodynamics, fluctuation theorems and molecular machines. *Rep Prog Phys*. (2012) 75:126001. doi: 10.1088/0034-4885/75/12/126001

15. Speck T. Stochastic thermodynamics for active matter. *Europhys Lett*. (2016) 114:30006. doi: 10.1209/0295-5075/114/30006

16. Ciliberto S. Experiments in stochastic thermodynamics: short history and perspectives. *Phys Rev X*. (2017) 7:021051. doi: 10.1103/PhysRevX.7.021051

17. Schmiedl T, Seifert U. Efficiency at maximum power: an analytically solvable model for stochastic heat engines. *Europhys. Lett*. (2007) 81:20003. doi: 10.1209/0295-5075/81/20003

18. Rana S, Pal PS, Saha A, Jayannavar AM. Single-particle stochastic heat engine. *Phys Rev E*. (2014) 90:042146. doi: 10.1103/PhysRevE.90.042146

19. Holubec V. An exactly solvable model of a stochastic heat engine: optimization of power, power fluctuations and efficiency. *J Stat Mech Theory Exp*. (2014) 2014:P05022. doi: 10.1088/1742-5468/2014/05/P05022

20. Tu ZC. Stochastic heat engine with the consideration of inertial effects and shortcuts to adiabaticity. *Phys Rev E*. (2014) 89:052148. doi: 10.1103/PhysRevE.89.052148

21. Bauer M, Brandner K, Seifert U. Optimal performance of periodically driven, stochastic heat engines under limited control. *Phys Rev E*. (2016) 93:042112. doi: 10.1103/PhysRevE.93.042112

22. Blickle V, Bechinger C. Realization of a micrometre-sized stochastic heat engine. *Nat Phys*. (2012) 8:143–46. doi: 10.1038/nphys2163

23. Quinto-Su PA. A microscopic steam engine implemented in an optical tweezer. *Nat Commun*. (2014) 5:5889. doi: 10.1038/ncomms6889

24. Martinez IA, Roldan Dinis L, Petrov D, Parrondo JMR, Rica RA. Brownian Carnot engine. *Nat Phys*. (2016) 12:67–70. doi: 10.1038/nphys3518

25. Argun A, Soni J, Dabelow L, Bo S, Pesce G, Eichhorn R, et al. Experimental realization of a minimal microscopic heat engine. *Phys Rev E*. (2017) 96:052106. doi: 10.1103/PhysRevE.96.052106

26. Albay JAC, Zhou ZY, Chang CH, Jun Y. Shift a laser beam back and forth to exchange heat and work in thermodynamics. *Sci Rep*. (2021) 11:4394. doi: 10.1038/s41598-021-83824-7

27. Zakine R, Solon A, Gingrich T, Van Wijland F. Stochastic stirling engine operating in contact with active baths. *Entropy*. (2017) 19:193. doi: 10.3390/e19050193

28. Saha A, Marathe R, Pal PS, Jayannavar AM. Stochastic heat engine powered by active dissipation. *J Stat Mech Theory Exp*. (2018) 2018:113203. doi: 10.1088/1742-5468/aae84a

29. Chaki S, Chakrabarti R. Entropy production and work fluctuation relations for a single particle in active bath. *Phys A Stat Mech Appl*. (2018) 511:302–15. doi: 10.1016/j.physa.2018.07.055

30. Chaki S, Chakrabarti R. Effects of active fluctuations on energetics of a colloidal particle: Superdiffusion, dissipation and entropy production. *Phys A Stat Mech Appl*. (2019) 530:121574. doi: 10.1016/j.physa.2019.121574

31. Saha A, Marathe R. Stochastic work extraction in a colloidal heat engine in the presence of colored noise. *J Stat Mech Theory Exp*. (2019) 2019:094012. doi: 10.1088/1742-5468/ab39d4

32. Holubec V, Steffenoni S, Falasco G, Kroy K. Active Brownian heat engines. *Phys Rev Research*. (2020) 2:043262. doi: 10.1103/PhysRevResearch.2.043262

33. Ekeh T, Cates ME, Fodor E. Thermodynamic cycles with active matter. *Phys Rev E*. (2020) 102:010101. doi: 10.1103/PhysRevE.102.010101

34. Kumari A, Pal PS, Saha A, Lahiri S. Stochastic heat engine using an active particle. *Phys Rev E*. (2020) 101:032109. doi: 10.1103/PhysRevE.101.032109

35. Szamel G. Single active particle engine utilizing a nonreciprocal coupling between particle position and self-propulsion. *Phys Rev E*. (2020) 102:042605. doi: 10.1103/PhysRevE.102.042605

36. Martin D, Nardini C, Cates ME, Fodor É. Extracting maximum power from active colloidal heat engines. *Europhys Lett*. (2018) 121:60005. doi: 10.1209/0295-5075/121/60005

37. Krishnamurthy S, Ghosh S, Chatterji D, Ganapathy R, Sood AK. A micrometre-sized heat engine operating between bacterial reservoirs. *Nat Phys*. (2016) 12:1134–8. doi: 10.1038/nphys3870

38. Larson RG. *The Structure and Rheology of Complex Fluids*. New York, NY: Oxford University Press (1999).

39. Wilson LG, Harrison AW, Poon WCK, Puertas AM. Microrheology and the fluctuation theorem in dense colloids. *Europhys Lett*. (2011) 93:58007. doi: 10.1209/0295-5075/93/58007

40. Démery V, Bénichou O, Jacquin H. Generalized Langevin equations for a driven tracer in dense soft colloids: construction and applications. *N J Phys*. (2014) 16:053032. doi: 10.1088/1367-2630/16/5/053032

41. Gomez-Solano JR, Bechinger C. Probing linear and nonlinear microrheology of viscoelastic fluids. *Europhys Lett*. (2014) 108:54008. doi: 10.1209/0295-5075/108/54008

42. Gomez-Solano JR, Bechinger C. Transient dynamics of a colloidal particle driven through a viscoelastic fluid. *N J Phys*. (2015) 17:103032. doi: 10.1088/1367-2630/17/10/103032

43. Berner J, Müller B, Gomez-Solano JR, Krüger M, Bechinger C. Oscillating modes of driven colloids in overdamped systems. *Nat Commun*. (2018) 9:999. doi: 10.1038/s41467-018-03345-2

44. Mohanty RP, Zia RN. Transient nonlinear microrheology in hydrodynamically interacting colloidal dispersions: flow cessation. *J Fluid Mech*. (2020) 884:A14. doi: 10.1017/jfm.2019.912

45. Toyabe S, Sano M. Energy dissipation of a Brownian particle in a viscoelastic fluid. *Phys Rev E*. (2008) 77:041403. doi: 10.1103/PhysRevE.77.041403

46. Vishen AS. Heat dissipation rate in a nonequilibrium viscoelastic medium. *J Stat Mech Theory Exp*. (2020) 2020:063201. doi: 10.1088/1742-5468/ab7e2f

47. Di Terlizzi I, Baiesi M. A thermodynamic uncertainty relation for a system with memory. *J Phys A Math Theor*. (2020) 53:474002. doi: 10.1088/1751-8121/abbc7d

48. Di Terlizzi I, Ritort F, Baiesi M. Explicit solution of the generalised Langevin equation. *J Stat Phys*. (2020) 181:1609–35. doi: 10.1007/s10955-020-02639-4

49. Zwanzig R. Nonlinear generalized Langevin equations. *J Stat Phys*. (1973) 9:215–20. doi: 10.1007/BF01008729

50. Brey JJ, Casado J. Generalized Langevin equations with time-dependent temperature. *J Stat Phys*. (1990) 61:713–22. doi: 10.1007/BF01027298

51. Romero-Salazar L, Velasco RM. Generalized Fokker-Planck equation with time dependent temperature. *Revista Mexicana de Fisica*. (1995) 41:358–64.

52. Gieseler J, Gomez-Solano JR, Magazzù A, Castillo IP, García LP, Gironella-Torrent M, et al. Optical tweezers: a comprehensive tutorial from calibration to applications. *Appl Opt Photon*. (2021) 13:74–241. doi: 10.1364/AOP.394888

53. Squires TM, Mason TG. Tensorial generalized Stokes-Einstein relation for anisotropic probe microrheology. *Rheol Acta*. (2010) 49:1165–77. doi: 10.1007/s00397-010-0490-5

54. Bird R, Armstrong R, Hassager O. *Dynamics of Polymeric Liquids, Volume 1: Fluid Mechanics*. 2nd ed. New York, NY: Wiley (1987).

55. Felderhof BU. Estimating the viscoelastic moduli of a complex fluid from observation of Brownian motion. *J Chem Phys*. (2009) 131:164904. doi: 10.1063/1.3258343

56. Indei T, Schieber JD, Córdoba A, Pilyugina E. Treating inertia in passive microbead rheology. *Phys Rev E*. (2012) 85:021504. doi: 10.1103/PhysRevE.85.021504

57. Xu K, Forest MG, Klapper I. On the correspondence between creeping flows of viscous and viscoelastic fluids. *J NonNewtonian Fluid Mech*. (2007). 145:150–72. doi: 10.1016/j.jnnfm.2007.06.003

58. Cordoba A, Indei T, Schieber JD. Elimination of inertia from a Generalized Langevin Equation: Applications to microbead rheology modeling and data analysis. *J Rheol*. (2012) 56:185–212. doi: 10.1122/1.3675625

59. Fischer P, Rehage H. Rheological master curves of viscoelastic surfactant solutions by varying the solvent viscosity and temperature. *Langmuir*. (1997) 13:7012–20. doi: 10.1021/la970571d

60. Ezrahi S, Tuval E, Aserin A. Properties, main applications and perspectives of worm micelles. *Adv Coll Interface Sci*. (2006) 128–30:77–102. doi: 10.1016/j.cis.2006.11.017

61. Paul S, Kundu A, Banerjee A. Active microrheology to determine viscoelastic parameters of Stokes-Oldroyd B fluids using optical tweezers. *J Phys Commun*. (2019) 3:035002. doi: 10.1088/2399-6528/ab0833

62. Paul S, Narinder N, Banerjee A, Nayak KR, Steindl J, Bechinger C. Bayesian inference of the viscoelastic properties of a Jeffrey's fluid using optical tweezers. *Sci Rep*. (2021) 11:2023. doi: 10.1038/s41598-021-81094-x

63. Wilhelm C, Gazeau F, Bacri JC. Rotational magnetic endosome microrheology: viscoelastic architecture inside living cells. *Phys Rev E*. (2003) 67:061908. doi: 10.1103/PhysRevE.67.061908

64. Vaippully R, Ramanujan V, Bajpai S, Roy B. Measurement of viscoelastic properties of the cellular cytoplasm using optically trapped Brownian probes. *J Phys Cond Matter*. (2020) 32:235101. doi: 10.1088/1361-648X/ab76ac

65. Raspaud E, Lairez D, Adam M, Carton JP. Triblock copolymers in a selective solvent. 2. Semidilute solutions. *Macromolecules*. (1996) 29:1269–77. doi: 10.1021/ma951172x

66. Zhu X, Kundukad B, van der Maarel JRC. Viscoelasticity of entangled lambda-phage DNA solutions. *J Chem Phys*. (2008) 129:185103. doi: 10.1063/1.3009249

67. Bellour M, Skouri M, Munch JP, Hébraud P. Brownian motion of particles embedded in a solution of giant micelles. *Eur Phys J E*. (2002) 8:431–6. doi: 10.1140/epje/i2002-10026-0

68. Grimm M, Jeney S, Franosch T. Brownian motion in a Maxwell fluid. *Soft Matter*. (2011) 7:2076–84. doi: 10.1039/c0sm00636j

69. Narinder N, Gomez-Solano JR, Bechinger C. Active particles in geometrically confined viscoelastic fluids. *N J Phys*. (2019) 21:093058. doi: 10.1088/1367-2630/ab40e0

70. Handzy NZ, Belmonte A. Oscillatory rise of bubbles in wormlike micellar fluids with different microstructures. *Phys Rev Lett*. (2004) 92:124501. doi: 10.1103/PhysRevLett.92.124501

71. Chapman CD, Robertson-Anderson RM. Nonlinear microrheology reveals entanglement-driven molecular-level viscoelasticity of concentrated DNA. *Phys Rev Lett*. (2014) 113:098303. doi: 10.1103/PhysRevLett.113.098303

72. Wulfert R, Oechsle M, Speck T, Seifert U. Driven Brownian particle as a paradigm for a nonequilibrium heat bath: effective temperature and cyclic work extraction. *Phys Rev E*. (2017) 95:050103. doi: 10.1103/PhysRevE.95.050103

73. Camati PA, Santos JFG, Serra RM. Employing non-Markovian effects to improve the performance of a quantum Otto refrigerator. *Phys Rev A*. (2020) 102:012217. doi: 10.1103/PhysRevA.102.012217

74. Holubec V, Marathe R. Underdamped active Brownian heat engine. *Phys Rev E*. (2020) 102:060101. doi: 10.1103/PhysRevE.102.060101

75. Cui B, Yang J, Qiao J, Jiang M, Dai L, Wang YJ, et al. Atomic theory of viscoelastic response and memory effects in metallic glasses. *Phys Rev B*. (2017) 96:094203. doi: 10.1103/PhysRevB.96.094203

76. Qian H. Fractional Brownian motion and fractional Gaussian noise. In: Rangarajan G, Ding M, editors. *Processes with Long-Range Correlations: Theory and Applications*. Berlin; Heidelberg: Springer-Verlag (2003) p. 22–33. doi: 10.1007/3-540-44832-2_2

77. Rodriguez RF, Fujioka J, Salinas-Rodriguez E. Fractional correlation functions in simple viscoelastic liquids. *Phys A Stat Mech Appl*. (2015) 427:326–40. doi: 10.1016/j.physa.2015.01.060

78. Sevilla FJ, Rodríguez RF, Gomez-Solano JR. Generalized Ornstein-Uhlenbeck model for active motion. *Phys Rev E*. (2019) 100:032123. doi: 10.1103/PhysRevE.100.032123

79. Gomez-Solano JR, Sevilla FJ. Active particles with fractional rotational Brownian motion. *J Stat Mech Theory Exp*. (2020) 2020:063213. doi: 10.1088/1742-5468/ab8553

80. Balland M, Desprat N, Icard D, Féréol S, Asnacios A, Browaeys J, et al. Power laws in microrheology experiments on living cells: comparative analysis and modeling. *Phys Rev E*. (2006) 74:021911. doi: 10.1103/PhysRevE.74.021911

81. Kobayashi Y, Tsukune M, Miyashita T, Fujie MG. Simple empirical model for identifying rheological properties of soft biological tissues. *Phys Rev E*. (2017) 95:022418. doi: 10.1103/PhysRevE.95.022418

82. Ferrer BR, Gomez-Solano JR, Arzola AV. Fluid viscoelasticity triggers fast transitions of a Brownian particle in a double well optical potential. *Phys Rev Lett*. (2021) 126:108001. doi: 10.1103/PhysRevLett.126.108001

83. Gomez-Solano JR, Samin S, Lozano C, Ruedas-Batuecas P, van Roij R, Bechinger C. Tuning the motility and directionality of self-propelled colloids. *Nat Commun*. (2017) 7:14891. doi: 10.1038/s41598-017-14126-0

84. Gomez-Solano JR, Roy S, Araki T, Dietrich S, Maciolek A. Transient coarsening and the motility of optically heated Janus colloids in a binary liquid mixture. *Soft Matter*. (2020) 16:8359–71. doi: 10.1039/D0SM00964D

85. Narinder N, Bechinger C, Gomez-Solano JR. Memory-induced transition from a persistent random walk to circular motion for achiral microswimmers. *Phys Rev Lett*. (2018) 121:078003. doi: 10.1103/PhysRevLett.121.078003

86. Lozano C, Gomez-Solano JR, Bechinger C. Active particles sense micromechanical properties of glasses. *Nat Mater*. (2019) 18:1118–23. doi: 10.1038/s41563-019-0446-9

87. Rings D, Schachoff R, Selmke M, Cichos F, Kroy K. Hot Brownian motion. *Phys Rev Lett*. (2010) 105:090604. doi: 10.1103/PhysRevLett.105.090604

88. Rings D, Chakraborty D, Kroy K. Rotational hot Brownian motion. *N J Phys*. (2012) 14:053012. doi: 10.1088/1367-2630/14/5/053012

Keywords: stochastic thermodynamics, heat engine, fluctuations, noise, viscoelasticity, memory effects, non-equilibrium processes, colloids

Citation: Gomez-Solano JR (2021) Work Extraction and Performance of Colloidal Heat Engines in Viscoelastic Baths. *Front. Phys.* 9:643333. doi: 10.3389/fphy.2021.643333

Received: 17 December 2020; Accepted: 04 March 2021;

Published: 29 March 2021.

Edited by:

Ramon Castañeda-Priego, University of Guanajuato, MexicoReviewed by:

Roberto Cerbino, University of Vienna, AustriaErick Sarmiento-Gomez, University of Guanajuato, Mexico

Copyright © 2021 Gomez-Solano. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Juan Ruben Gomez-Solano, r_gomez@fisica.unam.mx