Thermodynamic optimization subsumed in stability phenomena

In the present paper the possibility of an energetic self-optimization as a consequence of thermodynamic stability is addressed. This feature is analyzed in a low dissipation refrigerator working in an optimized trade-off regime (the so-called Omega function). The relaxation after a perturbation around the stable point indicates that stability is linked to trajectories in which the thermodynamic performance is improved. Furthermore, a limited control over the system is analyzed through consecutive external random perturbations. The statistics over many cycles corroborates the preference for a better thermodynamic performance. Endoreversible and irreversible behaviors play a relevant role in the relaxation trajectories (as well as in the statistical performance of many cycles experiencing random perturbations). A multi-objective optimization reveals that the well-known endoreversible limit works as an attractor of the system evolution coinciding with the Pareto front, which represents the best energetic compromise among efficiency, entropy generation, cooling power, input power and the Omega function. Meanwhile, near the stable state, performance and stability are dominated by an irreversible behavior.


the model and the regime
The LD model is a first-order irreversible deviation from a Carnot cycle, where irreversibilities are introduced only in the coupling between the working fluid and the hot and cold reservoirs at temperatures T c and T h , respectively 9 . The adiabatic processes are considered as instantaneous. The time dependent character of the model is introduced by means of the so-called irreversible coefficients c and h for each isothermal process, divided by the time these processes last, t c and t h , respectively. The fact that the baseline cycle is the Carnot one (whose efficiency does not depend on the working substance) when extending the reversible cycle to the irreversible framework allows for a structure independent description, where the specific properties of the working fluid are not relevant, and one can focus on the operation variables to optimize.
The input and output heats for a ref rigerator are Q c = T c �S[1 − � c /(�S t c )] and Q h = T h �S[−1 − � h /(�S t h )] , where S is the entropy change at the cold isotherm of the baseline reversible Carnot cycle. These expressions allow to decompose the heat exchange in a reversible component, T S , and an additional irreversible contribution proportional to the inverse of time. The mentioned irreversible coefficients c and h depend on the internal dynamics of the working fluid and thus, they are attached to intrinsic properties that for the purposes of optimization of the operation regime can be considered as constant. The total entropy change in the thermodynamic universe is given by �S tot = � c /t c + � h /t h , so that {t c , t h } → ∞ (or { c , h } → 0 ) the reversible performance is recovered. The main energetic quantities are the refrigerator coefficient of performance, COP ε = −Q c /(Q h + Q c ) , the cooling power, R = Q c /(t c + t h ) and the power input P in = W/(t c + t h ) = −(Q h + Q c )/(t c + t h ).
The function was proposed for a generic heat device as a trade-off figure of merit between the maximum useful energy and the unavoidable useful energy losses for a given energy input. For a RE, it is defined as the compromise between maximum cooling rate, R = Q c /t , and lost cooling rate for a fixed power input, P in . Then, the function is 14 where ε C = τ/(1 − τ ) is the COP of the Carnot cycle and τ ≡ T c /T h . Upper and lower bounds are obtained in the limits c → 0 and h → 0 , respectively In the symmetric case c = h 14 , The optimization of this function can be achieved equivalently using two different sets of time variables. In one case the total time and one partial contact time are used, allowing to reproduce behaviors typical from endoreversible and irreversible devices. In the other case the two partial contact times are considered. Different but complementary dynamical equations for the stability of this operation regime are obtained in each case. In the next section the first option is analyzed.

on time constraints, multi-objective optimization and stability
A fruitful representation of the system is achieved by using the operation total time and one partial contact time. This choice of optimization variables allows to recover endoreversible and irreversible behaviors which are tightly related with the multi-objective optimization of the RE, as will be shown later and it happens in the case of LD-HE's 45,46 . With this purpose the dimensionless variables, � c ≡ � c /� T , α ≡ t c /(t c + t h ) , t ≡ �S (t c + t h )/� T , with T = h + c that account for the system size are considered. From these definitions it is possible to work with dimensionless input and output (cooling power) heat fluxes 11 , A total entropy production, scaled with the size of the baseline Carnot cycle, can be obtained σ ≡ �S tot /( t �S) = [ � c /α + (1 − � c )/(1 − α)]/ t 2 . And from Eq. (4), the input power, P in = −˙ Q c −˙ Q h ; the COP, ε = R/ P in and � ≡ (2ε − ε C ) P in can be calculated. The optimization variables are α and t , meanwhile τ and c are fixed parameters, referring to external operation conditions and intrinsic material properties, respectively. The values of α and t that maximize are: A common feature in a variety of irreversible models is the appearance of loop-like power-efficiency parametric curves in heat engines and χ-ε curves for refrigerators. This has been considered as a characteristic signature of irreversibility. However, the so-called endoreversible model 54 consisting of a reversible Carnot engine irreversibly coupled to external heat reservoirs exhibits parabolic behaviors on those curves. In previous works it was shown that by fixing the value of α or t in the low dissipation model, it is possible to recover those types of behaviors. In the first case α = α M� while t can take values from 0 to ∞ (allowing to recover the reversible limit when t → ∞ ). For the latter t = t M and α takes values from 0 to 1 (this fixes the irreversibility of the system, since the reversible limit cannot be achieved). These two behaviors are discussed in detail in 13 , as part of a unified phenomenology for HE's and RE's and will become relevant in the present analysis of the stability-optimization relation. Some insights on the role of irreversible and endoreversible behaviors in the stability and optimization are also discussed for HE's in 45 .
Since the relaxation dynamics will be linked to an optimization process, a multi-objective optimization will be presented below. the best energetic compromise. It would be highly desirable to maximize the coefficient of performance, cooling power and simultaneously minimize entropy production and power input. However, no configurations can be found that fulfill all these requirements. It is common to search, instead, for the so-called Pareto front, which gives the best performance when one is looking to optimize simultaneously several objective functions 55 .
In this treatment two complementary outputs will be pursued: the location of the Pareto front, and second, the location of these points in the phase space (the Pareto optimal set) to compare them with the time constraints and the stability nullcline which will be discussed in in "Consecutive random perturbations" section.
The usual concept of dominance is used: in order to build the Pareto set a vector v = (v 1 , . . . , v n ) dominates another one w = (w 1 , . . . , w n ) if and only if v i ≥ w i ∀ i ∈ {1, . . . , n} (if one in looking for a maximum, ≤ for a minimum) and there is at least one j such that v j > w j . In other words, the improvement of any function will 1. In the phase space ( α, t ), the region of physical relevance is defined ( ε > 0 and R > 0).

2.
A random set of points in the phase space is obtained and the thermodynamic functions are evaluated (energetic space). 3. The set of non-dominated points in the energetic space is obtained, giving a provisional Pareto front. 4. From the corresponding Pareto optimal set (phase space) a convex region is defined and extended in order to cover a larger region for searching new points in the Pareto front. Details on the definition of the extended region are given below. 5. From the new region a new set of random points is proposed and a new set of non-dominated points in the energetic space is obtained.
In the present analysis, the Kullback-Leibler divergence (KLD) 56 is calculated between the distribution of points of the ith and the i − 1 th iterations. The radii to extend the search region in the phase space decreases with the KLD value. When this relative entropy is very small, there is no information gain in iterating more times, then, the search for new points in the Pareto optimal set stops. In Fig. 1a the Pareto front is shown using as objective function σ (minimization), ε , and R (maximization). The consideration of more functions, such as P in does not contribute to obtain new points. For a LD-HE it was found that the Pareto front was remarkably close to the endoreversible limit. In this case the match is higher. Thus, this limit offers the best energetic compromise. The corresponding points in the phase space are depicted in Fig. 1b, along with the boundary of the physical region of interest (where the efficiency and the cooling power are positive).
Stability dynamics. The problem of stability is frequently tackled through a first order equations system. By means of a linearization about the stable point, in many cases the resulting dynamics has the form f (y) = −dy/dt ≡ −ẏ , being y the dynamical variable and t the dynamical time. This is the kind of dynamical equation describe the stability of RC electrical circuit, overdamped spring, or any stable point in a potential well, V(y), near its minimum ( y * = 0 ) and where ẏ = −dV /dy. which is the well-known approximation to a harmonic oscillator. This is the kind of dynamics appearing in colloidal particles undergoing micrometric Carnot and Stirling heat engines through optical traps, that experimentally validated the low dissipation model 34,35 .
Below, it will be shown that if the operation regime of a time-periodic heat device has a steady state (an equilibrium point specified through the operation variables), for sufficiently small perturbations, it is expected that the (probably yet unknown) dynamics of the system can be expressed as a first order system. By means of a Taylor expansion of energetic functions (in this case Q c and ) around the steady state, a generic first order system can be transformed to a dynamical equation linking the operation variables and the energetic functions, allowing to effectively associate variations on the operation regime to fluctuations on the heat fluxes between the system and the heat reservoirs. In this way external perturbations can be linked to variations of the operation www.nature.com/scientificreports/ times (affecting energy fluxes). Notice that these perturbations come from external sources, and are not linked to the internal dynamics, already accounted by the coefficients c and h . As mentioned above, let us focus only on external perturbations on the operation regime (and for extension, to the variables that define such regime). The starting point is to assume that the system has an equilibrium point at the operation regime. With no further information regarding the specific energy transport, t and α are assumed to follow, within the first order scheme 57 , a typical (and the simplest) relation for an autonomous system in one dimension given by a dynamical equation The dynamical time, t has no dimensions and has a characteristic timescale that will be chosen later and A ∈ M 2×2 .
The natural energy flux associated to the time variable α would be Q c (variations of α will produce changes in the input heat). Meanwhile , a global energetic function will be linked to variations on the total time t . Analogous results are obtained if another global function is chosen, such as the so-called χ function (the equivalent to power output for heat engines 13 ).
From a first order expansion of ˙ Q c and around the steady state ( Y = 0 ) one obtains with J the Jacobian matrix C and D are positive constants determining the response speed to perturbations from the steady state that we will refer as the restitution strength 45 . Their values may depend on multiple characteristics, but usually the system size is the most important. Because large systems are more likely to respond slowly to perturbations on the control variables than small systems, the larger the system the smaller the values of C and D. From a dynamical perspective, their inverse values set a characteristic time scale, so that large values of C and D correspond to large restitution strength and short characteristic times. In the forthcoming analyses, results are referred to this time scale. Notice that the dynamical Eqs. (9) and (10) are not known a priori, but they are a mathematical requirement from stability relying, of course, on the assumption that the operation regime has a steady state. There are however, plausible arguments to think this is the case for several systems, since natural energy converters have displayed robustness, the capability to maintain and change operation regimes under stress in the operation conditions 62 .
In the linear approximation, the local stability of the above-given steady state is determined by the Jacobian matrix in Eq. (8) which by definition, in the M regime The determinant of J is zero as well as one eigenvalue; the other one is given by The relaxation time, t 1 (which refers to the variable α ), and the operation time, t M , are linked, being the first indication that the two phenomena are related. In the stability of a cyclic process one would require that t 1 ≤ t M , leading to the constraint Beyond the linear approximation presented above, the system of equations given by Eqs. (9) and (10) can be numerically solved. This dynamic, as in the case of HE's, produce stability basins, see Fig. 2. Trajectories inside www.nature.com/scientificreports/ the stability basin will arrive to the stationary state, the rest will diverge to a non-desirable state of infinite operation time and null cooling power. Even when the constraint provided by Eq. (12) does not involve D, it affects the stability basin shape which depends additionally on τ and c . Figure 2 shows the influence of the dissipation coefficient ( c ) on the basin of attraction for the complete solution obtained by solving numerically Eqs. (9) and (10) for given initial conditions/states (state after a perturbation). In Fig. 2a some representative trajectories for three dissipation symmetries, c = {0.05, 0.5, 0.95} , are depicted in the phase space. All the trajectories evolve in such a way that the variable t does not decrease. This is intuitively understood as this variable is related to the irreversibility of the system. In Fig. 2b, some trajectories are depicted over the -ε-σ surface. Notice that for symmetric dissipation, � c = 1/2 , the basin of attraction is less constraint and accepts larger perturbations. Figure 3 shows the influence of the relation between the magnitudes of C and D in the basin of attraction. Trajectories in the phase space after a perturbation for three cases: D = {1, C, 3/2C} are presented. In the first row of Fig. 3 the line integral convolution plot of Eqs. (9) and (10)  As can be seen in Fig. 4 for the representative case D = C (for C = D the results are similar) the endoreversible and irreversible limits have meaningful information regarding the performance of the refrigerator. The COP, cooling power, and entropy generation behaviors show that the trajectories tend to approach the endoreversible limit, which represents the better energetic compromise. After arriving to the endoreversible limit, the trajectories orbit around the stable state, however, those inside the stability basin display the smallest drawback in both εσ and those outside have the largest decay in the energetic efficiency and increment on entropy.
A characteristic of trajectories inside the basin of attraction is that the approaching to the endoreversible limit occurs with t < t M� ; then, in the arrival to the stable state their orbits are bounded by the irreversible limit (yellow curve) which is tangent to the basin of attraction in the steady state. Similar features were reported in the case of HE's. From this analysis it can be said that the endoreversible limit represents an attractor involved in an energetic self-improvement of the system and the irreversible limit bounds the basin of attraction. Both behaviors are relevant in the stability and thermodynamic improvement for both HE's and RE's. consecutive random perturbations. It is interesting to analyze the case where there is a lack of control in the optimization variables. Fluctuations on these time-variables will lead to variations on the heat fluxes and thus fluctuations around the steady state are expected. Since the internal dynamical nature of the system is already accounted through the macroscopic description given by the low-dissipation model, the source of these fluctuations for the problem at hand comes from outside the system.
For the analysis of these fluctuations the system at M conditions will undergo consecutive random perturbations along one cycle. To this purpose a cycle time will be divided by N equal sub-intervals of length t . The final state after one cycle is computed by solving the stochastic differential equation based on the proposed dynamic Eqs. (9) and (10), using a normally-distributed random variable as an additive white noise. Here, two independent stochastic variables ξ 1 and ξ 2 in the αt directions follow a 2-dimensional Gaussian distribution: Figure 2. Some representative trajectories given by solving numerically Eqs. (9) and (10). The symmetry in the dissipation coefficient c affects the shape of the basin of attraction. In (a) three cases are depicted: c = {0.05, 0.5, 0.95} , the symmetric case exhibits the largest stability basin. In (b) some trajectories are mapped into the energetic surface with τ = 3/5 . The red curves indicate trajectories inside the basin of attraction; dashed (red and blue online) curves denote those trajectories with initial conditions α → 0 and in continuous lines (red and blue) those with initial conditions α → 1. www.nature.com/scientificreports/ A representative case, γ = 2/�t 2 ≈ 2 × 10 −8 is considered; thus, the standard deviations are By using the Euler-Maruyama method 58 the points in the phase space are calculated iterating where (α 1 , t 1 ) = (α M� , t M� ) ; t is 10 −4 th of a cycle period t M , the response case C = D is studied. The energetic functions and the position in the phase space are averaged in each cycle. Then, it will be assumed that after each cycle-time the system evolves or it is driven to its time-periodic steady state, without any information from past cycles. This procedure is repeated for 10 5 trajectories. We have checked that this number of trajectories gives a good convergent statistics according to the relative entropy measured by Kullback-Leibler divergence 56 (see (15)   www.nature.com/scientificreports/ "A.1"), meaning that adding more trajectories do not offer further information and that the statistical description is conclusive. Figure 5a shows the final states in the phase space after each cycle: in green the points that ended inside the basin of attraction and in blue those that ended outside. For this kind of dynamics there exists a nullcline given by the curve t(α) for which α = 0 (see the returning point in the α direction in trajectories evolving from large values of α in Fig. 5). The consequences of this is that α evolves slowly around this curve. Note in Fig. 5a the relevant influence of the nullcline on the final point locations, which are agglomerated around it. The number of points inside the attraction basin is slightly larger than outside. After each cycle is completed, the arithmetic mean of each thermodynamic function is calculated. Points inside the basin of attraction tend to have shorter contact times with the cold reservoir and also shorter total operation times, while the opposite occurs for the points outside the basin of attraction. The Pareto front is depicted to compare the locus of the final states with the best thermodynamic compromise. The mean values obtained for all trajectories are depicted in Fig. 5b as horizontal lines. In Fig. 5c the corresponding probability distribution functions for , ε , and σ show the distinctive behavior for points inside and outside the stability region.
In Fig. 6a, the final states after each cycle are displayed in the -*ε-σ space. The endoreversible limit establishes an upper bound for all these trajectories and below the energetic configurations are located around the  www.nature.com/scientificreports/ irreversible limit. In Fig. 6b the averaged states over each stochastic trajectory/cycle are presented. Finally, in Fig. 6c the label "A" indicates the averaged energetic states for the mean values, and the label "B" the average of the final states, both points show that the performance of many cycles is very closed to that of the irreversible limit. This part of the analysis points out to a key role of time constraints in the obtaining of endoreversible and irreversible behaviors and their implications on the overall thermodynamic performance when the operation regime is stochastically perturbed. Many works have been reported on the endoreversible limit and its validity range. The stochastic results here obtained from induced perturbations on a steady state give new insights about the status and true nature of this endoreversible limit: it behaves as an attractor of the overall dynamics and describes the best compromise between the most relevant energetic functions, while the irreversible limit distinguishes trajectories inside or outside the stability basin and under perturbations it describes the statistical irreversible behavior induced over the system. These results are reproduced for different random variables distributions and C, D combinations.

Relaxation velocities and self-optimization
A second pair of time-variables can be used to study the optimization and stability of the LD RE. The optimization analysis is analogous to that of the previous Section, but the dynamics involved in the stability is different and as it will be shown below, the operation regime is described through a global stable state. This will provide a complementary vision of the linking between optimization and stability.
By introducing the dimensionless variables t c = t c �S/� h , t h = t h �S/� h , the input and output heats are Q c = Q c /(T h �S) , Q h = Q h /(T h �S) , respectively, and � = � c /� h (taking under consideration the system size). The heat exchanges per cycle can be written as in this way, the total entropy change of the thermodynamic universe per cycle is �S tot = ( t −1 h + � t −1 c ) and the entropy production σ = �S tot /( t c + t h ) . The information related to the internal dynamics is accounted by . Together with τ , it will remain as a fixed parameter. Due to the normalizing definitions the fluxes R ≡ Q c /( t c + t h ) , P in , σ and differ from the functions appearing in the previous section by a factor (1 + �) . The optimization of is achieved through the partial contact times t c and t h : The same upper and lower bounds for ε M� appearing in Eqs. (2) and (3) are obtained for → 0 and → ∞ , respectively.
the best energetic compromise. The same method to obtain the Pareto front (see in "The best energetic compromise" section) is used for this set of variables. The resulting Pareto front is depicted in Fig. 7a and coincides with that appearing in Fig. 1a with a scale factor (1 + �) . The Pareto optimal set is shown in Fig. 7b. The endoreversible and irreversible behaviors stemming from the time constraints are plotted in this pair of variables, see Fig. 7b.
As mentioned before, the two pairs of time variables discussed in the previous section and in this one are equivalent, the Pareto front, as well as the Pareto optimal set can be mapped from one system into the other and www.nature.com/scientificreports/ the exact results are obtained. Now, the dynamic equations for the second set of variables will be introduced, but one must notice that in the previous section the quantities that is tried to be fixed are energy fluxes, which is specially relevant in irreversible stationary processes. In this case, on the other hand, the stability is described in terms of a fixed energy (see below). For that reason, both dynamics are not equivalent.
Stability dynamics. Now, time variables are associated to heat exchanges between the system and the surroundings. It is then plausible to model the external perturbations on the operation regime as variations on the input and output heat, affecting engine operation time. For this analysis τ and , which are intrinsic properties of the system, will remain fixed. Each heat only depends on the associated partial contact time, i. e., Q c = Q c ( t c ) and Q h = Q h ( t h ) [see Eq. (17)], then, variations on the contact times can be effectively linked to variations on the corresponding input/output heats. The matrix formulation given in the previous section can now be addressed as a first-order one-dimensional uncoupled system 57  In a linear approximation the local stability of this steady state is determined by the eigenvalues, 1 and 2 , and eigenvectors of the Jacobian matrix: and since both 1 and 2 are real and negative, the operation regime steady state is stable. Relaxation times are defined as t 1 ≡ −1/ 1 and t 2 ≡ −1/ 2 , and by using Eqs. (18) and (23) they are given by From the above equations it is easy to check that relaxation is dominated by t 1 (linked to t c ). This is more noticeable for small values of τ and large values of , that is, when the dissipation coefficient at the cold reservoir is larger. By defining the total operation time t M tot ≡ t M h + t M c , that from Eq. (18) is given by the operation time can be compared with the total relaxation time, defined by t relax ≡ t 1 + t 2 , Since it is desirable that the system returns to the steady state within a cycle period t relax ≤ t M tot , A and B should fulfill the following condition Beyond the linear approximation, the system given by Eq. (22) can be solved numerically. Figure 8 shows the stream plot of the dynamics, which exhibits a stable point. Level curves for the dynamical velocity v dyn ≡ (d t c /dt) 2 + (d t h /dt) 2 = const are depicted. Notice that the velocity of relaxation is different in each www.nature.com/scientificreports/ quadrant. In this case the constraint t 1 = t 2 is used. Shorter relaxations are achieved for t c < t M� c and t h < t M� h , for other constraints, such as A = B this qualitative behavior holds, but the asymmetry between relaxation times will affect these contours. On the other hand, the Pareto optimal set, which do not depend on the restitution dynamics is depicted as a reference. The points in the Pareto optimal set with larger contact times (slower operation time) are closer to the minimum entropy generation state, meanwhile those with smaller contact times (faster operation time) are closer to the maximum R regime. One remarkable feature is that it appears a region where transition between velocity contours are more separated, meaning that restitution takes longer compared with changes in operation time. On the other hand, for another region transitions are faster. The largest restitution times are exhibited when only one partial time is perturbed. This will have consequences in the statistics of consecutive perturbations, as will be shown in next subsection.
In Fig. 9a the line integral convolution plot of the dynamic equations (22) over a random distribution of initial conditions is depicted, simulating streamlines of fixed arc length; darker shaded regions indicates smaller velocities. Figure 9b shows some representative trajectories with initial conditions at the border of the depicted region (initial state after a perturbation). Trajectories in each quadrant, labeled as I-IV are represented with different colors to emphasize differences in the energetic planes plotted in Fig. 9c-f. Trajectories in each quadrant evolve in a slightly different way, but this fact yields noticeable energetic repercussions. Of special relevance are regions I and III (colored in green and blue, respectively). Figure 9c shows the trajectories in the S tot -ε-surface. It is observed that trajectories in quadrants II and IV present large variations in the COP and entropy change before arriving to the stable state. This last feature could be considered as a non-desirable behavior. On the other hand, a kind of "benefit" in quadrants I and III is obtained since trajectories in these regions produce the best compromise between a given cooling power and the COP (larger values of ε for given R ) and the entropy change (the less entropy changes for a given R ) at the same time. These features are detailed in the 2D-plots of Fig. 9d-f showing that the trajectories in quadrants I (green) and III (blue) evolve directly towards the stable point in a very narrow region, meanwhile trajectories in quadrants II (red) and IV (yellow) present sudden changes of direction before arriving to the steady state.
The consequences of the dynamics on the system energetic properties could imply the use of a disadvantage, such as limited control, as a self-optimization mechanism throughout a biased control focused mostly in perturbations towards quadrants II and IV, to favor perturbations with both t c > t M� c and t h > t M� h or t c < t M� c and t h < t M� h . One could be interested not only in one perturbation but a series of them affecting the performance of one cycle. The influence of continuous random perturbations over a cycle will be addressed below in order to get more insights about this issue.
consecutive random perturbations. To analyze the effect of the dynamics in random perturbations a cycle time will be divided by N equal sub-intervals of length t . The final state after continuous perturbations along one cycle is computed by solving the stochastic differential equation based on the proposed dynamic equations (22), using a normally-distributed random variable as an additive white noise. Here, the two independent stochastic variables ξ 1 and ξ 2 in the t c and t h directions follow a 2-dimensional Gaussian distribution. By using the Euler-Maruyama method 58 phase space points are calculated iterating for the N steps where γ = 100 , the standard deviation σ t c,h is t M� c,h /10 . The initial state is ( t c 1 , t h 1 ) = ( t M� c , t M� h ) ; t is t M� tot /10 4 (25), so that after 10 4 steps one cycle time is covered. The constants A = 4 × 10(2 − τ )/(τ (1 − τ )) and B = 4 × 10/(1 − τ ) are used, so t relax = t M� tot /10 [see Eq. (27)] along with the condition t 1 = t 2 , so the dynamics is almost the same in the t c,h directions. After one cycle has ended, the system attains the time-periodic steady state and starts another random trajectory without any information regarding previous cycles. This is repeated for 5 × 10 4 trajectories/cycles. The statistical convergence has been tested using the Kullback-Leibler divergence 56 of the system energetic distributions [see equation ("A.2")]. The results are shown in Fig. 10.
In Fig. 10a the final states for each trajectory are depicted; colors indicate the quadrant of the system average position. The geometric centers of the points in each quadrant are displayed (see stars). The Pareto optimal set is shown as well. Notice that average behavior of the system is displaced towards the region of slow relaxation, following the direction of the Pareto front at which entropy change is minimized. By looking at the number of trajectories averaged in each region a trend is found: around 30.066% of the trajectories stayed in region I, 20.374% in region III and 49.56% almost equally distributed in regions II and IV, as it is emphasized in the close caption of Fig. 10a.
The mean values for the thermodynamic functions in each cycle are depicted in Fig. 10b-c for the energetic space involving S tot , ε , and R . The color representation is maintained. Points in region I have the best performance involving entropy production and efficiency, meanwhile trajectories in region III (region of fast relaxations) maintain closer to the steady state. In regions II and IV the energetic functions fluctuate very close to the steady state. From the above discussion the stability dynamics favors an energetic behavior where efficiency and entropy changes are improved at the cost of decreasing cooling power. www.nature.com/scientificreports/

concluding remarks and perspectives
The local stability of the Maximum Omega regime has been analyzed for the low dissipation refrigerator under two complementary dynamics. One with perturbations on one heat flux and the Omega function and another one with perturbations on the input/output heats. For each one an equivalent pair of time variables, which suitable describe the system and provide analogous optimizations reveal different features linking optimization with self-optimization. The restitution forces are modeled by means of dynamic equations stemming from requiring a steady state and small deviations of the input/output heat and fluxes. On the one hand, a dynamic with a basin of attraction is obtained, along with a nullcline which outside the attraction basin, leads the system to a non-physical region. The analysis of relaxation trajectories indicates that relaxation paths evolve in such a way that the thermodynamic performance of the system is improved. The system moves directly to the endoreversible limit due to the restitution dynamics. This endoreversible limit coincides with the Pareto front, which represents the best energetic compromise between all the relevant thermodynamic functions. Thus, the restitution dynamics induces an optimization process. For the case of many perturbations along one cycle the stability dynamics constitutes an irreversible mechanism, characterizing points inside and outside the stability basin.
On the other hand, for the second dynamics (and 2nd pair of variables), a global steady state is obtained and an analysis of relaxation velocities is presented. Fast and slow relaxations are obtained and together with the Pareto front (and the endoreversible limit) have an influence in describing the case of many perturbations. When the system moves to a worse energetic state, it improves due to the restitution dynamics in a fast manner, meanwhile if the system is perturbed to a better state, it evolves slowly to the steady state. This will produce a small departure in the average behavior of the system to a more optimum state regarding entropy generation and larger efficiency (COP).
The improvement due to perturbations on the steady state is related with the phenomenon of antifragility, in which a system exhibits an improvement under stress 62 . This could be of particular interest in medicine and biological systems. The extension of the low dissipation model to chemical engines 10 could also provide a path to test properties such as antifragility in biochemical processes.
The analysis here presented for RE's has been reported for low dissipation heat engines as well 45,46 , reinforcing the idea of a universal energetic characterization of heat devices as they are perturbed externally. Remarkably, the analyses of refrigeration systems are nowadays increasing its relevance since their coupling with other energy converters from different nature are being under study as an strategy for solving the problem of energy storage and heat recovery 59,60 . The study on a setup with fixed cooling power would be of interest, since in many applications the objective of refrigeration is to maintain this quantity as steady as possible.
Ideally, cyclic processes should be able to operate in a time-periodic steady state, one can make the idealization that the starting and ending points are well defined and analyze this as a unitary process. In reality, the switching between, aerobial respiration (linked to energetic trade-off functions) and anaerobial respiration (maximum power 61 ), to name an example, as well as the transient states in between in the changing of regime are beyond the simplified description of the present model. But it is precisely in the connection with this unknown mechanisms, where adaptation and robustness play a major role. So far, the results presented here for the self-optimization correspond to an emergent phenomenon (theoretically) that should be tested in specific structure dependent models.
The results presented in the present paper could open the window to the joint analysis of stability and optimization of coupled energy conversion systems with different purposes (heat engines, refrigerators, and heat pumps) and with different size scales, focusing on improving the strategies of production (aiming for a sustainability) and www.nature.com/scientificreports/ seize on energy. Possible applications of these results might involve colloidal particles constrained in an optical trap, describing Carnot or Stirling cycles. The equivalence between low dissipation model, irreversible Curzon-Ahlborn engine and Otto and Brayton engines might also allow for applications on solar energy heat devices. Additionally, this study could open new perspectives to analyze possible self-optimization or organization mechanism in specific systems and could be useful to understand non-revealed properties of non-equilibrium systems and their energetic bounds.

A Statistical convergence: Kullback-Leibler divergence
A.1 Dynamics involving Q c and . The Kullback-Leibler divergence, D KL allows to test statistical convergence. The value of D KL gives a measure of how distant are two distributions. If D KL = 0 , the information stemming from both distributions is the same. This is a relevant issue to demonstrate that the obtained trend is not due to the lack of additional data.
From the averaged cooling power R is computed of 10 5 trajectories, the interval between the largest and smallest R values is divided by √ N (rounded to the upper next integer 63 ) equal intervals, or bins, in this way, the same partition is used for the to compute the discrete probability distributions of the first k-thousand trajectories, ρ k are obtained, and the D KL is calculated comparing ρ k−1 with ρ k . The entropy divergence, D KL,k is given by giving a measure of how much information is gained by adding more trajectories. In Fig. 11, it is is shown that from 4 × 10 4 trajectories the statistical behavior does not vary significantly and adding trajectories will not provide much further information.