Work extraction and performance of colloidal heat engines in viscoelastic baths

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.


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 transfering 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,2,3,4,5,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,8,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,14,15,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,18,19,20,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,23,24,25]. In such colloidal systems, expansions and compresions 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 [26,27,28,29], Brownian engines with a self-propelled particle as working substance in contact with a viscous fluid [30,31,32,29] or in a suspension of passive Brownian particles [33] as a heat bath, as well as the realization of a colloidal Stirling engine in bacterial baths with tunable activity [34].
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 contant 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 [35]. 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 timedependent driving forces are exerted on an embedded particle [36,37,38,39,40,41], 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 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 [42,43,44,45].
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 correspoding 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.

Model
We consider a stochastic heat engine consisting of a Brownian particle in a viscoelastic fluid, whose motion is confined by a harmonic potential. Both the stiffness of the restoring force exerted by trapping potential and the temperature of the system can be varied in time according to a well-specified periodic protocol. For the sake of simplicity, we focus on the dynamics of a single coordinate of the particle, and assume that the rheological properties of the fluid remain in the linear viscoelastic regime. In such a case, the one-dimensional particle motion is described by the following generalized Langevin equation [46,47] where m is the mass of the particle, x(t) is its position at time t, whereas κ(t) and T (t) are the trap stiffness of the harmonic potential U(x, t) = 1 2 κ(t)x(t) 2 and the temperature at time t, respectively. Moreover, K(t − s) is a memory kernel that weights the effect of the previous history of the particle motion at time s ≤ t on its current drag force at time t due to the temporal correlations induced by the surrounding medium. The function K(t) is related to the relaxation modulus of the fluid, G(t), which quantifies its linear viscoelastic behavior in response to external deformations. For instance, for a spherical particle whose radius r is larger than the characteristic length-scales of the fluid microstructure, K(s) = 6πrG(s) [48]. In general, for a viscoelastic fluid G(t) is a function that decays to zero over a finite time-scale which is many orders of magnitude greater than those of simple viscous fluids. In addition, in Equation (1), ζ(t) is a Gaussian stochastic force which accounts for the thermal fluctuations in the fluid and satisfies an extended version of the fluctuation-dissipation theorem of the second kind It should be noted that, for an instantaneous memory Kernel, K(t) = 2γδ(t), and in the overdamped limit, i.e., for time-scales ≫ m γ , Equation (1) reduces to where the thermal noise ζ(t) simply satisfies Equation (3) describes the motion of a particle subject to a restoring force −κ(t)x(t) and to viscous dissipation by means of the instantaneous frictional force −γ dx(t) dt , with friction coefficient γ. This situation corresponds to a purely viscous coupling with the heat bath, which was explicitly considered in many of the models of single-particle heat engines reported in the literature [18,20,26,27,28,29,30,31,32].
In the following, we focus on the overdamped motion of a spherical particle of radius r embedded in a fluid with relaxation modulus consisting of a Dirac delta function plus an exponential decay which models the rheological response of several viscoelastic fluids, such as intracellular fluids [49], polymer solutions [50,51], wormlike micelles [52,39], and λ-phage DNA [53,39], 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 in Equation (1) is where the complex conjugate of its Fourier transform,K(ω) = ∞ −∞ dt e −iωt K(t), 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 (1) 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 [54,55], i.e., ∆x(t) 2 ≈ 2k B T γ 0 t, whose dynamics is simply described by Equation (3) with γ = γ 0 . 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., η = γ 6πr = η 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 non-Markovian Langevin Equation (8) corresponds to the Markovian one (3) with γ = γ 0 . The operation of the Brownian engine during a Stirling cycle of duration τ is depicted in Figure 1, 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. For 0 ≤ t < τ /2, the colloidal engine undergoes an isothermal expansion at high themperature T h by linearly decreasing the trap stiffness from κ M to κ m .
2. 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.
4. At t = τ , the temperature is suddenly raised to T c , while keeping the trap stiffness at κ(t = τ ) = κ M , i.e. an isochoric-like process, thus completing the full cycle.
Then, the cycle is repetead until the system reaches a time-periodic non-equilibrium steady state, which becomes independent of the choice of the initial condition Figure 1. Schematic representation of a Stirling cycle of period τ perfomed 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) = 1 2 κ(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 . At t = τ /2, the temperature is suddenly decreased to T c < T h , and kept at that value for τ /2 < t < τ , while linearly increasing the trap stiffness from κ m to κ M . The cycle is completed at t = τ , at which the temperature is again increased to T h . 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).
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 (14), 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 transfered 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 (13) and (14) are the ones needed for the calculation of the efficiency of the Stirling heat engine. They involve the variance of the particle position at an arbitrary time t, x(t) 2 , computed over an ensemble of independent realizations of the colored nosie ζ(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 [45]. Therefore, to address the problem of the performance of a Brownian Stirling heat engine described by Equations (2), (8), (11) and (12), we opt for numerical simulations of the corresponding stochastic dynamics.

Numerical solution
In order to compute the probability distributions of the work and the heat defined in Equations (13) and (14), 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 (1) for x(t) can be written as a linear system of two coupled Markovian Langevin equations with α defined in Equation (10). In Equation (18), D ∞ (t) is a short-time diffusion coefficient associated to the infinite-frequency friction coefficient γ ∞ 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 + τ 2 , T (t) remains constant. Accordingly, the rate of change of the time-dependent temperature in Equation (18) vanishes during each half a Stirling cycle, i.e., d dt [ln T (t)] = 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 [56,53,57,39]: η 0 = 0.040 Pa s and η ∞ = 0.004 Pa s, which correspond to α = 9. The diameter of the colloidal particle is set to r = 0.5 µm, while the maximum and minimum values of the trap stiffness during the Stirling cycle are chosen as κ M = 5 pN µm −1 and κ m = 1 pN µm −1 , respectively, which are easily accessible with optical tweezers. 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 s−100 s, which also covers characteristic values in actual experimental systems. We solve Equations (18) 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 (3) with η = η 0 = 0.040 Pa s and the rest of the involved parameters, namely, κ m , κ M , r, 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.1 s ≤ τ ≤ 10 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, 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 [31]. Therefore, in both the viscous and viscoelastic baths analyzed here, all the timescales are normalized by τ κ , whereas energies are normalized by k B T c . 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 Figure 2(A) and (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  2(A). In addition, large possitive 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 2(B). In this case, clear exponential tails remain for values of τ as large as τ = 16τ κ , as revealed in the inset of Figure 2(B). with a higher probability of occurence of negative heat fluctuations. Hence, regardless of the 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 .

Results and discussion
In Figures 2(C) and (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, but the width of the distribution decreases 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 2(C). On the contrary, the heat distribution, which also has asymmetric exponential tails, is rather narrow in the case of the viscous bath, and broadens with increasing fluid relaxation time in the case of a viscoelastic bath, thus converging to a τ 0 -independent curve for large τ 0 , see inset of Figure 2(D). A similar behavior of ̺(W τ ) and ̺(Q τ /2 ) is observed for values of the cycle duration τ larger than τ κ , as shown in Figures 2(E) and (F) for τ = 10τ κ . The fact that for a given cycle time τ , the work distribution of the Brownian engine is narrower in a viscoelastic bath as compared to that in a viscous bath with the same zeroshear viscosity, 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 ocurring in a viscous fluid, while for sufficiently large τ 0 significant deviations must take place, as verified in Figures 2(C) and (E) for τ 0 = 0.265τ κ and τ 0 = 265τ κ , respectively. To investigate the efficiency of a Brownian engine operating in a viscoelastic bath, we compute the mean power produced during a cycle whose dependence on the cycle time τ is plotted as solid lines in Figure 3(A) for some exemplary values of the fluid relaxation time τ 0 . Besides, we also compute the mean power for a Brownian engine in a Newtonian fluid bath (τ 0 = 0) with the same zero-shear viscosity (η = η 0 ), which is plotted as a thick dotted line in 3(A). It is important to note that, for all values of τ 0 , P τ exhibits a non-monotonic behavior as a function of τ . In particular, it 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., the power output vanishes (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 3(A), 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 location 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), we find τ * = 6.15τ κ , while for a viscoelastic fluid (τ 0 > 0), τ * is smaller and decreases with increasing τ 0 . In the inset of Figure 3(A) 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 (τ * = 6.15τ κ ), and monotonically decreases with increasing τ 0 . In this 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 3(A) for τ 0 = 0.084τ κ , 0.265τ κ , 0.84τ κ . Around τ 0 = τ κ , a conspicuous change in the dependence on τ 0 of the operation of the engine happens. Indeed, for τ 0 > τ κ , the stall time converges to the constant value τ * = 0.52τ κ as τ 0 increases, as verified in the inset of Figure 3(A) 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, for values of τ 0 > τ κ , which are consistent with increasingly pronounced viscoelastic behavior of the bath, the power output curves converge to a limiting curve as τ 0 /τ κ → ∞, as shown in Figure 3(A) for τ 0 = 2.65τ κ , 8.4τ κ , 26.5τ κ , 84τ κ , 265τ κ . 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 instaneously in response to the particle motion under the variations of the trap stifness. 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 energy dissipation 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 coefficientK * (ω → 0) = γ 0 = 6πrη 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 ω ∼ τ −1 0 the imaginary part ofK * (ω) is not negligible. Therefore, the value τ 0 ≈ τ κ marks a qualitative change in the energy exchanges 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 tô K * (ω → ∞) = γ ∞ for τ 0 ≫ τ κ according to Equation (7). As a consequence, in this limit the relavant dissipation timescale is γ ∞ /κ m , which is in general smaller than τ κ because η ∞ = η 0 (1 + α) −1 ≤ η 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 .
We confirm the previously described mechanism by computing the mean power output P τ 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, as a function of the cycle time τ . The resulting power curve is shown as a thick dashed line in Figure 3(A). Remarkably, we verifiy that this curve actually corresponds to the limiting curve to which the power output of a Brownian engine in a viscoelastic bath converges as τ 0 ≫ τ κ increases. 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 η ∞ , τ * = 0.52τ κ , as shown in the inset of Figure 3(A), see the horizontal solid line. Furthermore, in Figure 3(A) we also 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 , and is bounded by value of P τ in a Newtonian bath with constant 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.
In Figure 3(B), we plot the mean rate of heat absorbed by the engine during the isothermal expansion at temperature T h as a function of the cycle duration τ for some representative values of τ 0 . In this case, unlike the strong dependence of P τ on the fluid relaxation time, the behavior of J τ is only weakly dependent on τ 0 . We find that J τ is always possitive, i.e., heat is absorbed, and exhibits a typical dependence J τ ∼ τ −1 on the Stirling cycle time, which becomes more apparent at large cycle times. Indeed, in the inset of Figure 3(B) we show that the mean heat absorbed by the particle during the hot step of the cycle, − Q τ /2 , remains rather unaffected by the duration of the cycle, and approaches the quasi-static value − Q (τ →∞)/2 = 1.203k B T c (see horizontal thin solid line), which is explicitly given by and can be readily derived from Equation (14) [31]. For comparison, we also plot as a thick dashed line the corresponding curves for a Brownian engine in a fluid with constant viscosity η ∞ , to which both J τ and − Q τ /2 converge as τ 0 increases. 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. 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 4(a), with some efficiency curves plotted in Figure 4(A) as a function of τ for exemplary values of τ 0 . In Figure 4(a) we show that, as a consequence of the energy exchange in a viscoelastic bath discussed in the previos paragraphs, a marked qualitative change in the performance of the engine emerges around τ 0 ≈ τ κ . For the investigated values of the cycle time τ , the interval at which the engine has a positive efficiency is rather narrow in the region τ 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 for which ǫ > 0. Because in the model (1) 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 energy that affect the performance of the Brownian engine. Therefore, it is expected that the quasi-static Stirling efficiency is never exceeded at finite τ regardless of the relaxation time of the viscoelastic fluid. In Equation (26), ǫ C = 1 − Tc 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 4(B) 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, we find that the quasi-static Stirling engine can be reached in a viscoelastic bath for typical experimental values of the parameters of the system, as shown in Figure 4(B) for τ 0 > τ κ .

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 environement 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 would like to 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 [58] and quantum Otto refrigerators [59].
To the best of our knowledge, our work represents the first investigation on the effect of memory friction in the perfomance 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 thermodyamic cycles with finite period, as those described by power law kernels and fractional Brownian noise [60,61,62,63], which describe the mechanical response of diverse soft materials encountered in the biological realm [64,65]. 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 [66]. Similar ideas could be extended to active Brownian heat engines [29] functioning in complex fluids, which could be implemented in practice by. e.g. light-activated colloids in non-Newtonian liquids [67,70,68,55,69] and hot Brownian particles [71,72,73].

Author Contributions
J.R.G.-S., concieved the model, carried out the numerical simulations, analysed the results, and wrote the manuscript.

Funding
This work was supported by UNAM-PAPIIT IA103320.