Optimal Control of Uniformly Heated Granular Fluids in Linear Response

  Citation: Ruiz-Pino, N.; Prados, A. Optimal Control of Uniformly Heated Granular Fluids in Linear Response. Entropy 2022,24, 131. https://doi.org/10.3390/e24010131 Academic Editor: Antony N. Beris Received: 15 November 2021 Accepted: 12 January 2022 Published: 16 January 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). entropy Article Optimal Control of Uniformly Heated Granular Fluids in Linear Response Natalia Ruiz-Pino and Antonio Prados * Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain; [email protected] *Correspondence: [email protected] Abstract: We present a detailed analytical investigation of the optimal control of uniformly heated granular gases in the linear regime. The intensity of the stochastic driving is therefore assumed to be bounded between two values that are close, which limits the possible values of the granular temperature to a correspondingly small interval. Specifically, we are interested in minimising the connection time between the non-equilibrium steady states (NESSs) for two different values of the granular temperature by controlling the time dependence of the driving intensity. The closeness of the initial and target NESSs make it possible to linearise the evolution equations and rigorously—from a mathematical point of view—prove that the optimal controls are of bang-bang type, with only one switching in the first Sonine approximation. We also look into the dependence of the optimal connection time on the bounds of the driving intensity. Moreover, the limits of validity of the linear regime are investigated. Keywords: optimal control; granular fluids; linear response; Sonine approximation; bang-bang controls 1. Introduction The study of granular media, beyond its own theoretical interest, is particularly important for industrial applications such as improving their transport or storage. Granular materials are discrete clusters of macroscopic particles that exhibit two fundamental features. First, collisions between particles are inelastic, so that energy is not conserved: It monotonically decreases with time if there is no external mechanism that injects energy into the system. Second, thermal energy is many orders of magnitude lower than the characteristic potential energy, making thermal fluctuations largely irrelevant for the behaviour of granular systems [1]. In the simplest model for granular fluids, particles are d -dimensional smooth, hard spheres of mass m that undergo inelastic binary collisions. In each collision, the tangential component of the relative velocity is unchanged, whereas the normal component is reversed and shrunk by a factor α , 0 ≤α≤ 1, which is termed the restitution coefficient. Energy is only kinetic, and the energy dissipated in each collision is thus proportional to 1 −α2—the elastic limit corresponds to α= 1. In the undriven system, after a few collisions per particle, the so-called homogeneous cooling state (HCS) is reached [ 2 – 7 ], in which the system remains homogeneous and the granular temperature T —basically the average kinetic temperature—monotonically decreases following an algebraic decay, the Haff law [2]. In order to allow the system to reach a stationary state, an energy injection mechanism is needed. A simple but also relevant situation is the uniformly heated granular fluid [ 8 , 9 ] that we consider throughout this work. Therein, independent white noise forces act on the particles of the granular fluid, the intensity of which is characterised by a parameter χ≥ 0 related to the variance of the stochastic force. The granular fluid reaches a non-equilibrium steady state (NESS) in the long-time limit, in which the system remains homogeneous. Therein, the energy injected by the stochastic thermostat balances—in average—the energy Entropy 2022,24, 131. https://doi.org/10.3390/e24010131 https://www.mdpi.com/journal/entropy Entropy 2022,24, 131 2 of 29 loss in collisions and the value of the granular temperature depends on the intensity of the driving, whereas higher-order cumulants of the velocity are independent thereof. The steady state is always stable, at variance with the HCS, which is known to be unstable for perturbations of large enough wavelength [ 4 , 10 , 11 ] (On another note, velocity correlations stemming from the dissipative character of collisions are relevant in certain physical situations, e.g., for the understanding of the total energy fluctuations in the HCS [ 12 , 13 ]). The uniformly heated granular gas has been extensively studied, both its properties at the NESS [8,9,14,15] and its dynamical evolution [16–18]. Granular systems are intrinsically out-of-equilibrium systems. Their dissipative dynamics entail that their velocity distribution function (VDF) is non-Gaussian, even in the long-time limit in which a hydrodynamic state, independent of the initial condition, is reached. This is true for both the HCS in the undriven case and the NESS in the uniformly heated situation. The non-Gaussianities of the VDF are essential to understand the behaviour of granular fluids and are incorporated to the picture by implementing a Sonine expansion [ 3 ] of the Enskog–Fokker–Planck equation. This leads to an infinite hierarchy of equations for the cumulants, which is typically closed by introducing the so-called first Sonine approximation: Only the fourth cumulant or excess kurtosis a2 is retained—higher order cumulants are neglected. Therein, the granular temperature and the excess kurtosis obey a system of two coupled ordinary differential equations, the accuracy of which for describing the dynamical evolution of the granular fluid has been validated in many works, e.g., [ 8 , 9 , 15 , 17 – 21 ]. In this context, especially relevant are those analysing memory effects such as the Kovacs hump or the Mpemba crossing, in which non-Gaussianities are key to facilitate their emergence [18–21]. Only very recently has the possibility of controlling the dynamical evolution of granular systems been analysed [ 22 ]. This might be surprising at first sight, since the control of physical systems has been considered for some time in different physical contexts, such as quantum mechanics [ 23 – 26 ] and statistical mechanics [ 27 – 32 ]. A paradigmatic case of control of a mesoscopic system is that of an optically trapped colloidal particle [27,28,31,33–39] . When the confining potential is harmonic, the time dependence of the stiffness of the trap κ(t) can be externally controlled, and one aims at optimising the connection between two given equilibrium states, corresponding to different values of the stiffness of the trap—i.e. , the colloidal particle is being confined or deconfined. Here, optimising means that some relevant physical observable (irreversible work, entropy production, connection time, . ..) is minimised. The time-dependent stiffness κ(t) plays the role of the control function—sometimes together with the temperature of the bath, which can be changed in an effective way by adding a random force [ 40 , 41 ]. The control problem is greatly simplified by the following three features. First, the initial and target states are equilibrium states, so that their corresponding probability distribution functions (PDFs) are perfectly known. Second, the PDF is Gaussian for all times, so that it is completely characterised by its average and variance. Third, the evolution equations for the average and the variance are exactly solvable in closed form. The delay in posing the problem of controlling granular systems probably stems from the challenging character of the control problem in this case, both at the conceptual and mathematical level. None of the three simplifying features above, holding for the harmonically trapped Brownian particle, is present in granular fluids. First, the initial and target states are NESS, and their PDFs are only approximately known. Second, the PDF is non-Gaussian for all times. Third, the evolution equations are non-linear and thus not exactly solvable. It is interesting to compare the situation in the granular case described above with the one appearing in other paradigmatic system, the Brownian gyrator [ 42 – 44 ]. Although the initial and final states are also NESSs in that case, the PDF is Gaussian for all times, and the evolution equations for the relevant moments can be exactly solved. Thus, the control problem of this system is simpler, although only non-optimal connections have been worked out, to the best of our knowledge [45]. Entropy 2022,24, 131 3 of 29 One may thus pose the problem of connecting two NESSs of the granular fluid corresponding to different values of the driving intensity χ , χi and χf , i.e., to different values of the granular temperature Ti and Tf . The control function here is the intensity of the driving χ(t) . We are interested in the time optimisation problem, i.e., to find the protocol χ(t) —starting from (and ending at) the desired initial (and target) NESS—that minimises the connection time between the initial and final states. This kind of time optimisation problem is important from a fundamental point of view and also has relevance for applications. For the connection between equilibrium states, related problems emerge in the optimisation of irreversible heat engines [ 46 ], the analysis of the Mpemba effect [ 21 , 47 – 49 ], and the optimisation of the relaxation route to equilibrium [50–52]. The limiting situation in which all the power of the stochastic thermostat is available, i.e., 0 ≤χ<∞ , was investigated in Ref. [ 22 ] within the first Sonine approximation. Despite the challenges mentioned above, the unboundedness of the control makes it possible to give analytical predictions for the connecting time—the evolution equations are heavily simplified in the limiting cases χ= 0 and χ=∞ . In addition, the analytical predictions were compared with numerical simulations, and an excellent agreement was found. Indeed, this was expected because the accurateness of the first Sonine approximation is a wellestablished fact for the smooth hard-sphere granular gas, even for the study of subtle behaviours such as the Kovacs or Mpemba memory effects [ 18 – 21 ]. This is the reason why we follow a purely analytical approach in this paper. In this work, we analytically investigate the more realistic case in which the driving intensity is bounded between two values, χmin ≤χ≤χmax . In order to make analytical progress, we consider the linear response regime, in which χmin and χmax are close: This allows us to linearise the evolution equations and make exact—in the linear response limit—predictions for the optimal connecting time as a function of the bounds (χmin , χmax) . The linearisation of the equations also allows us to employ rigorous mathematical results of optimal control theory (OCT) and to check that the underlying hypotheses are fulfilled, a program that was unattainable in the non-linear case [ 22 ]. Moreover, we also explore the limits of validity of the linear response regime, by taking the double limit (χmin  1, χmax  1 ) and comparing the obtained behaviour with those for the non-linear case with unbounded driving [22]. The structure of this paper is as follows. In Section 2, we put forward the model, write the evolution equations for the temperature and the excess kurtosis, and linearise them around the final NESS. Section 3is devoted to the derivation of the optimal controls, in the sense of minimising the connection time. The trajectories of the temperature and the excess kurtosis—both as functions of time and in the phase plane—for the optimal controls are analysed in Section 4. The dependence of the minimum connection time on the bounds of the driving is the subject of study of Section 5. We investigate the limits of the validity of the linear response approximation as the bounds in the driving are loosened in Section 6. Finally, a discussion of the obtained results is presented in Section 7. The Appendices deal with some technicalities that are omitted in the main text. 2. The Model Our system is a granular fluid with number density n , comprising N d -dimensional ( d= 2, 3) hard-spheres of mass m and diameter σ (hard discs in d= 2). Specifically, we consider smooth inelastic hard spheres. Collisions between them are binary, and the post-collisional velocities (v0 1,v0 2)are given in terms of the pre-collisional ones (v1,v2)by: v0 1=v1−1+α 2(v12 ·b σ)b σ,v0 2=v2+1+α 2(v12 ·b σ)b σ, (1) where b σ is the unit vector along the direction joining the centre of the particles and α is the restitution coefficient, 0 ≤α≤ 1. In addition, the system is heated by a stochastic thermostat, i.e., a white-noise force Fi independently acts on every particle verifying hFi(t)i=0, Fi(t)Fj(t)=m2ξ2δijδ(t−t0),∀i,j=1, . . . , N, and ∀(t,t0). Entropy 2022,24, 131 4 of 29 In the first Sonine approximation that we employ throughout, the system is described by two variables, the granular temperature T and the excess kurtosis a2 . Their definitions in terms of moments of the velocity are: Dv2E=dT m,a2=d d+2v4 hv2i2−1. (2) As stated in the introduction, the system reaches an NESS in the long-time limit due to the balance—on average–of the energy input and dissipation. The stationary values of T and a2are given by: T3/2 s=mξ2 ζ0(1+3 16 as 2)≡χ,ζ0=2nσd−1(1−α2)πd−1 2 √mdΓ(d/2), (3) as 2=16(1−α)(1−2α2) 73 +56d−24dα−105α+30(1−α)α2, (4) Note that as 2 is independent of the thermostat intensity, as measured by χ , it only depends on (d , α) . From the kinetic equation, the following coupled system of ordinary differential equations (ODEs) are obtained—see, e.g., [17,19]: ˙ T=ζ0χ(1+3 16as 2)−T3/2(1+3 16a2), (5a) ˙ a2=2ζ0 Th(T3/2 −χ)a2+BT3/2(as 2−a2)i, (5b) where the parameter Bis given by [17,19] B=aHCS 2 aHCS 2−as 2 ,aHCS 2=16(1−α)(1−2α2) 25 +2α2(α−1) + 24d+α(8d−57); (6) aHCS 2is the value of the excess kurtosis in the HCS [8,9]. 3. Optimal Control in Linear Response Above, we have considered that the driving intensity of the thermostat χ is constant. In general, it may be time-dependent, a certain given function of time χ(t) that determines the externally enforced driving program. Looking at the evolution Equation (5a,b) for (T , a2) in the light of OCT, this means that χ(t) is the control function. In this paper, we consider the following control problem: the connection of two NESS, i.e., bringing the system from an initial state (Ti , a2i =as 2) , to a target, final one (Tf , a2f =as 2) , by engineering a suitable driving program χ(t) . Moreover, we would like to perform this connection in the shortest possible time. The case in which all the power of the thermostat is available has been considered in Ref. [ 22 ]. Here, we analyse the more realistic case in which the driving intensity is bounded between two limiting values, χmin ≤χ(t)≤χmax , with χmin ≥ 0 and χmax <∞. In order to solve the control problem analytically, we restrict ourselves to the linear response regime, i.e., χmin and χmax are close—and so are Ti and Tf . To look into the dynamics of the system, it is preferable to introduce scaled variables as follows: t∗=ζ0T1/2 ft,T∗=T Tf ,χ∗=χ T3/2 f ,A2=a2 as 2 . (7) In this way, we have defined dimensionless time t∗ , granular temperature T∗ , and driving χ∗ ; moreover, scaling the excess kurtosis with its steady value simplifies our analysis (Both aHCS 2 and as 2 change sign for α= 1 /√2 , so that a2 typically changes sign with the inelasticity. On the other hand, the scaled variable A2 always remains positive). For these scaled variables, we have the following evolution equations: Entropy 2022,24, 131 5 of 29 ˙ T=χ(1+3 16as 2)−T3/2(1+3 16as 2A2), (8a) ˙ A2=2 Th(T3/2 −χ)A2+BT3/2(1−A2)i. (8b) We have omitted the superscript ∗ in the dimensionless variables in order to simplify the notation since; from now on, these are the variables used. The term χ( 1 +3 16 as 2) on the right hand side (rhs) of (8a) represents the energy injection due to the action of the thermostat, while the term −T3/2( 1 +3 16 as 2A2) collects the energy losses due to the inelastic collisions. Of course, if χ is kept constant and equal to its target value, i.e., χ(t) = 1 ∀t≥ 0, the system reaches the NESS (Ts= 1, As 2= 1 ) in the long-time limit, consistently with our discussion in the previous section. In linear response, we thus write: T=1+δT,A2=1+δA2,χ=1+δχ, (9) with δT 1, δA2 1, and δχ  1. Note that, to be consistent, we must assume that δχmin , δχmax  1. This allows us to linearise the evolution equations of T and A2 as follows: d dt δT δA2=β −2δχ +−3 2β1−β 3−2BδT δA2, (10) where we have defined β≡1+3 16as 2. (11) Now it is δχ that plays the role of the control function, δχmin ≤δχ ≤δχmax. Once the evolution equations are linearised, the control problem is stated as follows: We would like to bring the system from the initial NESS corresponding to δT(t=0) = δTi,δA2(t=0) = 0, (12) to the target NESS δT(tf) = 0, δA2(tf) = 0, (13) in the minimum possible time tf . Moreover, the system remains stationary for t< 0 and t>tf : This means that for t< 0, we have prepared the system in the NESS with the initial value of the temperature by driving it with the corresponding intensity and that for t≥tf , the driving intensity for the target temperature is applied, i.e., δχ(t) = δχi=3 2δTi,t<0, δχ(t) = 0, t≥tf. (14) Equation (10) is linear in both the variables (δT , δA2) and the control function δχ , and therefore the rigorous theorems for linear control systems are applicable—see, for example, chapter III of Ref. [ 53 ]. For our specific situation in which δχmin ≤δχ ≤δχmax , these theorems ensure that the optimal protocol that minimises the connection time tf is of bang-bang type with at most one change. That is, δχ(t) is piece-wise continuous, taking either the value δχmax or δχmin and presenting, at most, one jump between these two values in the time window ( 0, tf) (More specifically, this result stems from Theorem 10 in Section 17 of Ref. [ 53 ], and we check that the hypotheses of this theorem are fulfilled in Appendix A). These kind of bang-bang optimal protocols arise in different physical situations [ 22 , 23 , 54 – 56 ]. In general, bang-bang protocols emerge as the optimal ones when Pontryagin’s Hamiltonian is linear in the controls—i.e., when the evolution equations are linear in the controls, although they may be non-linear in the relevant physical variables [30,57–59]. Entropy 2022,24, 131 6 of 29 To be able to determine the optimal protocol, we must distinguish two cases according to the initial temperature, a global cooling process ( Ti>Tf= 1, δTi> 0) and a global heating process (Ti<Tf=1, δTi<0): • For δTi> 0, CH protocol: In the time window [ 0, tJ) , the driving δχ is set to its minimum value δχmin (cooling), whereas in the time window [tJ , tf) , it is set to its maximum δχmax (heating): δχ(t) =        δχi,t<0, δχmin, 0 ≤t<tJ, δχmax,tJ≤t<tf, 0, t≥tf. (15) • For δTi< 0, HC protocol: In the time window [ 0, tJ) , the driving is set to its maximum value χmax, whereas in the time window [tJ,tf), it is set to its minimum χmin: χ(t) =        δχi,t<0, δχmax, 0 ≤t<tJ, δχmin,tJ≤t<tf, 0, t≥tf. (16) The switching time tJ , 0 ≤tJ≤tf , will be determined later as a function of the parameters of the problem, i.e., as a function of (δχi , δχmin , δχmax) . As already stated above, the values of δχ for t< 0 and for t≥tf ensure that the system starts from the NESS with T=Ti and, after the application of the bang-bang protocol, remains in the target NESS with T=Tf= 1, for both the CH and HC protocols. At first, there is no clear reason to assign the CH protocol to the case δTi> 0 and the HC protocol to the case δTi< 0 (aside from the analogy with the full-thermostat-power case analysed in Ref. [ 22 ]). In order to show that this is indeed the case, one needs to study the behaviour of the trajectories swept by the point in the phase plane (δA2 , δT) . We defer this analysis until Section 4. 3.1. Ti>Tf= 1 Cooling-Heating Bang-Bang In this section, we integrate the solution of the system in two time windows: the first one, [ 0, tJ) , when δχmin is applied, and a second one [tJ , tf) , when δχmax is applied. We also determine the time tJ as well as the value of the variables δTJ and δA2J at that time. The point (δA2J , δTJ) constitutes the set of initial conditions for the control system in the second window. Equation (10) is inhomogeneous, due to the term proportional to δχ on its rhs. (It is only homogeneous when δχ = 0, i.e., when the control is set to the constant value χf=T3/2 f corresponding to the final temperature.) Over each time window, δχ(t) = δχext , where the subscript “ ext ” includes both bangs, δχext =δχmin (first window) and δχext =δχmax (second window). The inhomogeneity can be thus understood as the system relaxing towards the NESS corresponding to δχext . Let us denote by Text = 1 +δText the temperature corresponding to the NESS reached when the system is driven with constant intensity χext = 1 +δχext . Since the steady value of the excess kurtosis does not depend on the driving intensity, we have only to subtract: δText = (χext)2/3 −1=2 3δχext +O(δχext)2(17) from δTto make the system homogeneous. Thus, we define: δ˜ T≡δT−δText. (18) The homogeneous system for δ˜ Tand δA2reads: Entropy 2022,24, 131 7 of 29 d dt δ˜ T δA2=−3 2β1−β 3−2Bδ˜ T δA2. (19) The eigenvalues (−λ1,−λ2)and eigenvectors (v1,v2)of this system are given by: λ1=1 23 2β+2B+k>0, v1=1 62λ2−3β 6, (20a) λ2=1 23 2β+2B−k>0, v2=1 62λ1−3β 6, (20b) where we have introduced the parameter: k≡λ1−λ2=s3 2β−2B2 +12(1−β)>0. (21) With the definitions above, both λ1and λ2are positive, and λ1>λ2. Now we can write the solution in both time windows, separately, because δ˜ T is different over each one. In the first step of the bang-bang, t∈[0, tJ), where δχ(t) = δχmin: δ˜ T δA2=δT−2 3δχmin δA2=C1v1e−λ1t+C2v2e−λ2t. (22) In the second step of the bang-bang, t∈[tJ,tf), where δχ(t) = δχmax: δ˜ T δA2=δT−2 3δχmax δA2=C3v1e−λ1t+C4v2e−λ2t. (23) The constants (C1 , C2 , C3 , C4) are obtained by imposing the initial conditions in each time window. For t=0, we have the initial condition (12), which determines C1and C2: C2=−C1=3δTi−δTmin k=2δχi−δχmin k. The point at the final time tJof the first time window is: δTJ=δ˜ TJ+δTmin =C1hv1(1)e−λ1tJ−v2(1)e−λ2tJi+2 3δχmin, (24a) δA2J=C1e−λ1tJ−e−λ2tJ. (24b) The initial conditions for Equation (23) are supplied by (δTJ , δA2J) . Therefore, we can obtain (C3,C4)as a function of the switching time tJ: δTJ−δTmax δA2J=C3v1e−λ1tJ+C4v2e−λ2tJ. (25) Note that (C3 , C4) also depend on the bounds of the driving (δχmin , δχmax) through δTmin and δTmax . By imposing that we have to reach the target state, i.e., Equation (13) , we write: −δTmax 0=C3v1e−λ1tf+C4v2e−λ2tf. (26) Equations (25) and (26) are four equations for the four unknowns (C3 , C4 , tJ , tf) , which thus provide us with the solution to the control problem. Solving for C3and C4, we obtain: C3=2δχmax k  δχtot δχi−δχmin eλ1tJ−1 δχtot δχi−δχmin eλ2tJ−1  λ1 k ,C4=−C3e−ktf. (27) Entropy 2022,24, 131 8 of 29 We have introduced the total amplitude of the allowed interval for the driving: δχtot ≡δχmax −δχmin =χmax −χmin >0. (28) In this way, the final time tfis given as a function of the switching time tJ: tf=1 kln  δχtot δχi−δχmin eλ1tJ−1 δχtot δχi−δχmin eλ2tJ−1 =tJ+1 kln  1−δχi−δχmin δχtot e−λ1tJ 1−δχi−δχmin δχtot e−λ2tJ , (29) which is in turn given by the solution of the implicit equation: δχtot1−δχi−δχmin δχtot e−λ2tJλ1 k=δχmax1−δχi−δχmin δχtot e−λ1tJλ2 k. (30) The set of Equations (29) and (30) provides an analytical solution for the minimum connection time tf in the CH protocol, which is valid in the linear approximation we are considering in this paper. 3.2. Ti<Tf= 1 Heating-Cooling Bang-Bang Let us now start from an initial state with δTi< 0. The analysis of this case is similar to that just carried out for δTi> 0, but the order of the bangs is reversed. In the first time window, [ 0, tJ) , the maximum driving δχmax is applied, whereas in the second time window, [tJ , tf) , the minimum driving δχmin is applied. Therefore, the homogenisation procedure for the temperature δTis also reversed. In the first window t∈[0, tJ)we have: δ˜ T δA2=δT−2 3δχmax δA2=C5v1e−λ1t+C6v2e−λ2t, (31) whereas in the second window t∈[tJ,tf)it is: δ˜ T δA2=δT−2 3δχmin δA2=C7v1e−λ1t+C8v2e−λ2t. (32) The initial conditions are given by Equation (12) . Inserting them into Equation (31) , we obtain: C5=−C6=3δTmax −δTi k=2δχmax −δχi k. (33) The evaluation of Equation (31) at time tJ gives the initial condition for the second time window. Taking into account the difference in the δ˜ T variables in Equations (31) and (32) , due to our switching the value of the driving intensity at t=tJ , one obtains (C7 , C8) in terms of tJ and also of the bounds (δχmin , δχmax ) in complete analogy with the CH protocol. In addition, at the final time tf , one must impose that the system reaches the target NESS, i.e., Equation (13) , which provides the two extra equations needed to determine the switching time tJand the connection time tfas functions of the system parameters. The result is: C7=2δχmin k  δχtot δχmax−δχieλ1tJ−1 δχtot δχmax−δχieλ2tJ−1  λ1 k ,C8=−C7e−ktf, (34) for (C7,C8)in terms of (tJ,tf): tf=1 kln  δχtot δχmax−δχieλ1tJ−1 δχtot δχmax−δχieλ2tJ−1 =tJ+1 kln  1−δχmax−δχi δχtot e−λ1tJ 1−δχmax−δχi δχtot e−λ2tJ , (35) Entropy 2022,24, 131 9 of 29 for the minimum connection time in terms of the switching time, and the following implicit equation: δχtot1−δχmax −δχi δχtot e−λ2tJλ1 k=−δχmin1−δχmax −δχi δχtot e−λ1tJλ2 k, (36) for tJ . In complete analogy with the CH case, the set of Equations (35) and (36) gives the minimum connection time for the HC protocol in the linear approximation. Note that exchanging δχmin ↔δχmax (which entails δχtot → −δχtot ) leads from Equations (35) and (36) to Equations (29) and (30)—and vice versa. Above, we have derived analytical expressions for different physical variables of interest, for both the CH and HC bang-bang protocols, in the linear response approximation. More specifically, we have (i) the complete description of the trajectory followed by the system in the phase plane, i.e., the time evolution of the point (δA2(t) , δT(t)) , and (ii) the switching time tJ and the minimum connection time tf . The linear response approximation has allowed us to obtain analytical predictions as functions of all the relevant physical parameters: not only of the initial temperature Ti , as measured by δχi= 3 δTi/ 2, but also of the bounds of the driving (χmin , χmax) , as measured by (δχmin , δχmax) . Therefore, it is interesting to inspect the behaviour of the obtained expressions as a function of (δχi , δχmin , δχmax) in order to understand the response of the system to the optimal control designed. 4. Trajectories for the Temperature and the Excess Kurtosis In this section, we look into the trajectories of the temperature and the excess kurtosis, to understand the need for a two-step bang-bang protocol on a physical basis. The time evolution of δT and δA2 is presented in Figure 1, for both the CH and the HC cases ( δT solid lines, δA2 dashed lines). First, let us analyse the CH protocol (left panel), i.e., δTi> 0 (also δχi> 0). Therein, δT relaxes to δTmin under the action of δχmin in the time window [ 0, tJ) (without reaching it, since the relaxation at constant driving lasts for an infinite time). Simultaneously, δA2 starts to increase from its steady value, equal to zero, because Equation (10) implies the following: d dt δA2t=0+ =−2δχmin +3δTi=2(δχi−δχmin), (37) which is non-negative because δχmin ≤ 0 ≤δχi . (Otherwise the connection of the two NESS would be impossible, as rigorously proven in the next section—and in agreement with physical intuition.) A decrease (an increase) in the granular temperature makes the VDF separate from (closer to) the Gaussian shape, i.e., the scaled excess kurtosis A2 correspondingly increases (decreases). Once the target temperature Tf= 1 is reached inside this first time window, i.e., the temperature curve crosses the horizontal axis δT= 0, the action of δχmin cannot be interrupted by setting the thermostat intensity to unity because δA2> 0 and the system is not in the target NESS. This is why we must let δT continue to drop to a value δTJ such that δTJ< 0, associated with a kurtosis value δA2J> 0. This point (δA2J , δTJ) is determined by the condition that, at the end of the subsequent relaxation with δχmax in the time window [tJ , tf) , δA2 and δT must simultaneously reach their target value (zero). Second, we analyse the HC case (right panel), the discussion is completely analogous and thus we summarise it in the following. In the first time window [ 0, tJ) with δχmax , the horizontal axis δT= 0 is crossed at some time smaller than tJ , but it is necessary to continue applying δχmax to overshoot it, since δA2< 1 for that time and the system has not reached the target NESS. Once more, the point (δA2J , δTJ) is determined by the condition that, at the end of the second time window with δχmin , both variables simultaneously vanish. The need of a two-step protocol, and the order of the bangs, can also be understood —maybe more clearly—by looking at the trajectories in phase space. The trajectories of Entropy 2022,24, 131 16 of 29 2 4 6 8 10 12 14 1 2 3 4 5 Figure 8. Convergence to the non-linear expression (38) as the bounds go to more extreme values for the HC case. Symbols code of the data shown are the same as in Figure 7. The breakage of the convergence “from above” to the non-linear result is clearly seen: For extreme enough values of the bounds, the linear time becomes smaller than the non-linear prediction for a full-strength thermostat. It is interesting to remark that as α decreases (i.e., as the inelasticity increases), less extreme values of χmax and χmin are needed to provoke the inversion. In other words, the linear approximation breaks down for less extreme bounds. For example, let us consider α= 0.1, which can be regarded as a high-inelasticity case—as opposed to the lowinelasticity case α= 0.9. Fixing χmin = 0, the connection times of the linear approximation become shorter than those given by Equation (38) for δχmax ≥ 9.2 in the range of δTi∈[−0.01, 0] (or, equivalently, δχi∈[− 0.015, 0 ] ), smaller than the value δχmax ≥ 9.8 for α= 0.9. This trend with α of the bounds leading to the inversion of the convergence “from above” can be understood by recalling that as 2is a decreasing function of α. Consequently, the importance of the heating term in the evolution equation of the temperature (8a) , χ1+3 16 as 2 , increases as α is lowered: A smaller value of χ is needed to obtain the same value of the heating term. 7. Discussion Our work improves the understanding of the optimal control of driven granular gases. The results obtained in this paper complement and enrich those obtained in Ref. [ 22 ] for a full-strength thermostat. The inclusion of bounds in the driving, χmin ≤χ≤χmax raises non-trivial questions that have been answered by our study, such as the range of initial and target temperatures that can be connected. Our investigation has been carried out in the linear response regime, i.e., the initial and target states are close enough—and so are the bounds of the driving χmin and χmax . This allows us to linearise the evolution equations around the final (target) NESS. The linear response approximation leads to a set of evolution equations that are linear both in the control function and the dynamical variables—more precisely, in their deviations from their target values. Therefore, we obtain a linear control problem that can be completely solved. A rigorous mathematical theorem ensures that the optimal control, minimising the connection time, is of bang-bang type with at most one switching: i.e., the optimal control comprises two time windows [ 0, tJ) and [tJ , tf) , with the control being equal to one of its limiting values, either χmin or χmax , in the first time window [ 0, tJ) and changing to the other limiting value at the switching time tJ . Therefore, two types of Entropy 2022,24, 131 17 of 29 bang-bang protocols arise, depending on the order of the bangs: χmin followed by χmax , which we have termed CH, or χmax followed by χmin , which we have termed HC. We have shown that the CH protocol is the optimal one when the initial temperature Ti is larger than the final one Tf , whereas the HC protocol is the optimal one in the reverse situation, Ti<Tf. We have investigated the behaviour of the connection time as a function of the bounds in the driving intensity. This study has allowed us to elucidate the range of initial and final temperatures that can be connected. We have shown that the final temperature has to lie between the temperatures Tmin and Tmax , where Tmin ( Tmax ) is the steady temperature corresponding to the constant driving χmin ( χmax ). On the other hand, the initial temperature may lie outside the interval [Tmin , Tmax] , and the connection is still possible: The minimum connection time is still finite when the upper bound χmax crosses the initial temperature Ti>Tf (or the lower bound crosses the initial temperature Ti<Tf ) (It could be argued that, still, the most relevant physical situation corresponds to the case Ti∈[Tmin , Tmax] because one needs to prepare the system in the initial NESS). In addition, we have explored the limits of the validity of the linear response approximation we have employed throughout. We have done this by loosening the restrictions on the bounds χmin and χmax . Specifically, we have analysed the behaviour of our (linear response) prediction for the minimum connection time, tf , as χmin is decreased to very small values, and χmax is increased to very large values. This behaviour has been compared with the minimum connecting time for the non-linear case tn` f , which was obtained when the thermostat has its full strength, χmin = 0 and χmax =∞ . Specifically, we have compared the linear time with the asymptotic expression for tn` f for small temperature jumps—in which the corresponding connection times are also very small [ 22 ]. For the CH case, we have found that tf tends to tn` f always “from above”, tf>tn` f . This is logical, since the largest set of controls—such as that of the full-power thermostat—should lead to the shortest connection times. However, for the HC case, we have found that the tendency from above towards tf is broken for large enough values of χmax . This marks a limit of validity for the linear response approximation in this case. The asymmetry between the CH and HC protocols can be understood on a physical basis: In the latter case, heating precedes cooling, and thus the temperature departs from the vicinity of Tf for high enough χmax —whereas in the former, cooling precedes heating, and the system remains closer to the target state even when χmax becomes large. In order to further look into the behaviour described in the previous paragraph, we have looked into the regime of short connecting times tf 1 within the linear response framework. Note that our linear response predictions for tf , as given by Equations (29) and (35) , contain all the powers of δχi/δχtot , δχmax/δχtot , δχmin/δχtot . Linear response assumes that both δχi 1, δχmax  1, and δχmin  1, but the ratios between one another are in principle of the order of unity. It is only when δχi is much smaller than δχmax and δχmin that the connection time becomes small. In this regime, we have obtained a simple approximate expression for tf valid to the lowest order in δχi , which also depends on the bounds in the driving. This approximate expression always gives connection times that are longer than that for the full-power thermostat, both for the CH and HC cases. This means that the inversion of the tendency “from above” towards tn` f comes from higher-order terms in the ratios δχi,max,min/δχtot. Our work also opens the door to finding new optimal controls for other non-equilibrium systems. For example, let us look at a colloidal particle moving in the vicinity of a minimum of the trapping potential—which can be thus considered to be harmonic. In that case, the temperature of the thermal bath in which the particle is immersed plays the role of the driving intensity. Interestingly, the temperature of the bath can be effectively increased by adding a random forcing that can be modelled as a Gaussian white noise [ 40 , 41 ]. In this way, the effective temperature changes from Tmin (room temperature) to Tmax (thousands of kelvins). The similitude of the mathematical framework, linear evolution equations and bounded control, makes it appealing to analyse the optimal connection—also in the sense Entropy 2022,24, 131 18 of 29 of minimising the connection time—in that case and compare the corresponding results with those derived here. It is also relevant to extend the optimal controls derived here—and also of those for a full-power thermostat in Ref. [ 22 ]—to more complex situations, such as rough granular gases [ 60 – 64 ]. Translational and rotational velocities are correlated, even for small roughness [ 60 , 61 ]. In addition, the number of relevant variables increases: The values of the translational and rotational temperatures are, in general, different, and additional cumulants emerge (See, for example, [ 62 – 64 ] for the general kinetic framework). The impact of these features on the optimal controls is thus an interesting prospect for future work. Author Contributions: Conceptualization, A.P.; data curation, N.R.-P.; formal analysis, N.R.-P.; funding acquisition, A.P.; investigation, N.R.-P.; methodology, A.P.; project administration, A.P.; software, N.R.-P.; supervision, A.P.; validation, N.R.-P. and A.P.; visualization, N.R.-P. and A.P.; writing—original draft, N.R.-P. and A.P.; writing—review and editing, N.R.-P. and A.P. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by project PGC2018-093998-B-I00 FEDER/Ministerio de Ciencia e Innovación–Agencia Estatal de Investigación (Spain). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data supporting the published results are available from the authors upon reasonable request. Acknowledgments: A.P. thanks Pablo Maynar for useful discussions. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Appendix A. Maximum Principle for Linear Systems: Verifying Hypothesis Let us consider the linear, both in the variables xand the controls u, control system: dx dt =Ax +Bu, (A1) in which x:R→Rn , u:R→U , where the control set U is a m -dimensional parallelepiped, and A and B are two matrices of suitable dimensions. Now we analyse the problem of bringing the system from xi to xf in the minimum possible time tf , which is known as the time optimisation problem. The columns of the matrix B are denoted by bj , and we introduce the assumption that the set of vectors {bj , Abj , A2bj , ... An−1bj} constitutes a basis of Rn for each j= 1, . . . , m . Under this controllability hypothesis, we can formulate the following theorem: Theorem A1. If all the eigenvalues of A are real, then the optimal controls are bang-bang, i.e., they take the most extreme values of their definition domain and present, at most, n −1switchings. Our system (10) perfectly fits into the framework given by Equation (A1) , with the identifications: x=δT δA2,u=δχ,A=−3 2β1−β 3−2B,B=β −2. (A2) Therefore, n= 2 and m= 1, with U≡[δχmin , δχmax] . We know that the eigenvalues of the matrix A are real, since they are given by (−λ1 , −λ2) in Equation (20) . Thus, the theorem above applies, and the controls are bang-bang with at most one switching if the vectors {b1 , Ab1} form a basis of R2 . The determinant of the matrix with columns b1 and Ab1 is ∆=4(βB−1)>0, since βB>1 for all α. In fact, as seen in Figure A1,∆increases with α. Entropy 2022,24, 131 19 of 29 0 0.2 0.4 0.6 0.8 0 50 100 150 Figure A1. Discriminant ∆ as a function of the restitution coefficient α . Both the d= 2 (solid line) and d= 3 (dashed) cases are shown. The discriminant remains positive for all values of α , guaranteeing that the optimal control is of the bang-bang type with at most one switching. Appendix B. Order for the Bangs In this Appendix B, we prove that the optimal protocol for the case Ti<Tf= 1 is of HC type. A completely analogous proof links the case Ti>Tf= 1 to the CH protocol. We proceed by showing that one can only reach NESSs with Ti<Tf= 1, i.e., with δTi< 0, making use of an HC protocol, it is impossible with a CH bang-bang. The idea of the proof is based on rigorously establishing that the qualitative behaviour of the motion of the system in the phase plane (δA2 , δT) is the one depicted in Figure A2. We start by analysing the shape of the cooling curve, which starts from a NESS ( 0, δTi) . For the cooling steps, δT(t) and δA2(t) are given by their respective expressions in Equation (24) , with the substitution tJ→t. Therefore, their time derivatives are: d dt δT=2δχi−δχmin kλ1v1(1)e−λ1t−λ2v2(1)e−λ2t, (A3) d dt A2=2δχi−δχmin khλ1e−λ1t−λ2e−λ2ti. (A4) Note that dA2/dt|t=0= 2 (δχi−δχmin)> 0. On the one hand, δT(t) monotonically decreases from δTi= 2 δχi/ 3 for t= 0 to δTmin = 2 δχmin/ 3 for t→∞ because d(δT)/dt does not vanish for t> 0. In fact, the possible extremum of δT occurs at a time kt1=log v1(1)λ1 v2(1)λ2 , which either does not exist (for α> 1 /√2 ) or is negative (for α< 1 /√2 ). On the other hand, d(δA2)/dt vanishes at a time t0given by: t0=1 klog λ1 λ2 >0, (A5) i.e., δA2 increases from zero to positive values in the interval [ 0, t0) , then reaches a maximum at t=t0and decreases back to zero for t0<t. Entropy 2022,24, 131 20 of 29 Figure A2. Qualitative picture of the heating and cooling trajectories in the phase plane. Curves for t> 0 are drawn with solid lines, curves for t< 0 with dashed lines; heating ones ( δχmax) are in red, cooling ones ( δχmin) in blue. The common tangent to the heating and cooling curves at the origin is represented by a black dotted line. The system starts cooling from an initial state δTi< 0. Let us assume that the bang-bang protocol is of CH type. In the first step of the bang, the system follows the blue cooling curve: If it is allowed to relax during an infinite time, it reaches the NESS ( 0, δTmin =2 3δχmin) over the vertical axis. The cooling must be interrupted at some time tJ> 0, where the driving is switched to δχmax : Since the system must reach the target NESS at the origin, it needs to move over the branch of the heating curve corresponding to t< 0 (red dashed line). However, this is impossible, since this heating curve is always above the tangent line and does not intersect the cooling curve. Therefore, it is not feasible to drive the system to the origin using a CH protocol for δTi<0. The above discussion entails that the motion of the point (δA2 , δT) along the cooling curve in the phase plane follows indeed the shape depicted by the blue solid line in Figure A2. An analogous study shows that the shape of the heating trajectories that starts from a NESS, i.e., from the vertical axis δA2= 0, must be like the red solid line in Figure A2: δT increases monotonically for t> 0, and δA2 starts decreasing, reaches a minimum, and afterwards increases back to zero. Once we know the qualitative behaviour of the curves (heating and cooling) that start from an NESS, let us apply a cooling process to an initial state ( 0, δTi< 0 ) . We know that all of these cooling trajectories must be contained in the region between the vertical axis and the cooling curve starting from the origin, since phase plane trajectories for the same value of δχ cannot intersect. The only way for the system to reach the point ( 0, 0 ) with a CH protocol is that the cooling trajectory beginning at ( 0, δTi< 0 ) intersects the heating trajectory ( δχmax ) that crosses the origin ( 0, 0 ) , i.e., to the branch of the red solid line corresponding to negative times (red dashed line). Our aim is to prove that this trajectory cannot enter the region described above, where all the cooling curves starting from ( 0, δTi< 0 ) are confined. To prove this, note that the cooling and heating trajectories going through the origin have a common tangent, since d dt δT d dt δA2 =−βδχ 2δχ =−β 2<0 (A6) is independent of δχ (black dotted line). We proceed to prove that the heating curve that goes through a point P over the tangent line is always above this tangent line and, therefore, cannot intersect any cooling curve. It is enough to show that the slope of the heating curve Entropy 2022,24, 131 21 of 29 that goes through an arbitrary point over the tangent line is always larger (lower in absolute value) than the slope of the tangent. To do this, we take the point Pover the tangent: δTP=−β 2δA2P, (A7) and making use of Equation (10), we have that: d dt δTP=βδχmax −3 2βδTP+ (1−β)δA2P=βδχmax +3 4β2+1−βδA2P, (A8) and d dt δA2P=−2δχmax +3δTP−2BδA2P=−2δχmax −3 2β+2BδA2P. (A9) We have that 3 4β2+1−β>0: 3 4β2+1−β>0=3 41+3 16as 22 −3 16as 2=3 41+1 8as 2+27 1024(as 2)2>0, (A10) and 3 2β+ 2 B> 0, ∀(α , d) . Therefore d dt δTP> 0 and d dt δA2P< 0 over the tangent line, and the slope of the heating curve on this point is negative: mP= d dt δTP d dt A2P =βδχmax + (3 4β2+1−β)δA2P −2δχmax −(3 2β+2B)δA2P <0. (A11) Now we compare it with the slope of the tangent: mP+β 2= 3 2β2+1−βδA2P−β3 2β+2BδA2P 2−2δχmax −(3 2β+2B)δA2P=[4βB−3(1−β)]δA2P 2δχmax + (3 2β+2B)δA2P≥0 (A12) for δA2P≥0 because 4βB−3(1−β)>0. Then: mP>−β 2. and the heating curve cannot cross any cooling curve (since all of them are under the tangent line). We conclude that it is not possible to drive the system from ( 0, δTi< 0 ) to the origin with a CH protocol. On the other hand, the HC protocol starting from ( 0, δTi< 0 ) and ending up at the origin is indeed feasible, as shown in Figure A3. The cooling curve that goes through the origin—specifically its dashed branch—divides the semi-plane δA2< 0 into two parts. The initial point ( 0, δTi< 0 ) and the NESS for the heating part of the protocol, ( 0, δTmax) , lie at different sides thereof. As a consequence, the heating curve for the first bang—over which δT monotonically increases and δA2 has only one minimum—intersects at only one point the cooling curve for the second bang, giving rise to the optimal connection. Entropy 2022,24, 131 22 of 29 Figure A3. Qualitative picture of the HC protocol in the phase plane. Curves for t> 0 are drawn with solid lines, curves for t< 0 with dashed lines; heating ones ( δχmax) are in red, cooling ones ( δχmin) in blue. The system starts heating from an initial state δTi< 0, corresponding to point I . In the first step of the bang, the system follows the red heating curve: If allowed to relax during an infinite time, it reaches the NESS ( 0, δTmax =2 3δχmax) over the vertical axis. The heating is interrupted at the point J over the cooling curve for t< 0 (blue dashed line), where the driving is switched to δχmin . The optimal connection thus comprises the arcs IJ and JF. Appendix C. Limit Values for the Intensity of the Thermostat In this Appendix C, we are interested in studying the range of values for the intensity of the thermostat that make it possible to connect the initial and target NESS (for both the CH and HC protocols). In particular, we would like to discern whether the connection is possible for any pair of values (χmin , χmax) or there appears a region in parameter space that make the connection impossible. We analyse the CH case (cooling, Ti> 1) in detail, since the analysis of the HC case follows completely analogous lines. Physically it seems clear that, in order to reach the final temperature Tf= 1 from an initial temperature Ti> 1, it is imperative that the minimum intensity of the driving verifies χmin < 1—the minimum value of the thermostat intensity has to be smaller than that corresponding to the final temperature. However, what about χmax ? Is the connection always possible as long as χmax > 1? Or, on the contrary, is there a lower bound that makes it impossible to connect the two NESS? The limit values (χmin , χmax) , beyond that the connection is no longer possible, are those that bring about a divergent minimum connection time tf. In the following, we show how the line of reasoning above gives answers to the questions posed: the physical intuition on the limit value of χmin is correct, tf diverges in the limit as δχmin → 0 − , and the limit value of χmax is also unity, tf diverges in the limit as δχmax → 0 + . In order to prove these statements, we will follow the following procedure: to elucidate the behaviour with δχmin ( δχmax ), we keep δχmax ( δχmin ) fixed and progressively increase δχmin (decrease δχmax ) from negative (positive) values until the minimum connection time diverges. From Equation (29) , which gives tf as a function of δχi , δχmin and δχmax in the CH case, we can infer the values of δχmin that make tf diverge. This divergence only comes about when either the numerator of the logarithm tends to ∞ or the denominator tends to 0. We analyse both possibilities in the following. For the numerator to diverge, either δχmin =δχi> 0 (recall that we are studying the CH protocol) or tJ→∞ . In the latter case, tf always diverges (Physically , it is evident that the connection time cannot be shorter than the switching time, tf≥tJ . Mathematically, the divergence of the numerator always wins because λ1>λ2 and tf∼tJ ; the time spent in the second part of the bang-bang becomes negligible as compared with tJ ). The switching time tJ is determined by Equation (30) , which tells us that when tJ→∞: Entropy 2022,24, 131 23 of 29 δχtot →δχmax. (A13) Therefore, tJ→∞ when δχmin → 0 − , in agreement with the physical intuition described above. This makes it unnecessary to study the other possibility of divergence of tf , δχmin →δχi> 0. Let us explore the second possibility, i.e., the vanishing of the denominator of the logarithm in Equation (29), which occurs when: tJ→1 λ2 lnδχi−δχmin δχtot . (A14) Note that the numerator of the logarithm is positive for this value of tJ because λ1>λ2 . When Equation (A14) holds, the left hand side (lhs) of Equation (30) vanishes. As a consequence, it is δχmax → 0 + , since the factor accompanying it on the rhs of Equation (A14) is basically the numerator of the logarithm in Equation (29) . Therefore, tf→∞ when δχmax →0+. Wrapping things up, our analysis above implies that it is always possible to connect two non-equilibrium steady states as long as δχmin < 0 and δχmax > 0, that is, the lower (upper) bound of the thermostat intensity, χmin ( χmax ) is below (above) the one corresponding to the final state—i.e., unity, with our choice of variables. Therefore, there are no additional regions in parameter space that do not allow for connecting the two NESS. As already said above, the HC case (heating, Ti< 1) is treated in a completely analogous way, with the same conclusion: χmin < 1 and χmax > 1, with the roles of these limitations exchanged with respect to the case CH. Appendix D. Bang-Bang for the Non-Linear Case In the linear case, Equation (10) , the number of switchings of the bang-bang control is given by the theorem in Appendix A. In the first Sonine approximation employed to describe the granular gas, we have two variables and thus only one switching. In this way, either the CH protocol or the HC protocol is that minimising the connection time between the initial and final NESS. In the non-linear case, Equation (8) , the optimal connection is also of bang-bang type [ 22 ]. The evolution equations—despite being non-linear in the temperature—are linear in the intensity of the driving, and Pontryagin’s maximum principle [ 53 , 54 ] ensures that the optimal control minimising the connection time is bang-bang. However, the number of switchings from one extreme value of χ to the other is not known. The simplest two-step bang-bang protocols were investigated in Ref. [ 22 ], but it was not proved that the two-step bang-bangs led to the minimum time. Here we present such a proof. In this Appendix D, we consider more complex bang-bang processes, with more than two steps. We show that the connection time for this more complex protocols is always longer than that for the two-step bang-bangs. For the non-linear case, we look into the movement of the system in the phase plane (A2 , T) (instead of (δA2 , δT) ). Let us first focus on the case Ti< 1, illustrated by Figure A4. The system starts from a point I= ( 1, Ti) and ends up at the target point F= ( 1, Tf= 1 ) . The heating curve—with χ=χmax =∞—passing through I (red solid line) and the cooling curve—with χmin =0—passing through F (blue solid line) intersect at the point J . The two-step bang-bang is formed by joining the portion of the heating curve joining I and J and the portion of the cooling curve joining J and F , i.e., the arcs IJ and JF. A four-step bang-bang is shaped as follows. Let us consider a point K belonging to the portion of the heating curve from I to F , interrupt the heating at this point and switch the driving to χmin = 0. Then the system starts to sweep the cooling curve passing through K (blue dashed line). At the point L , the cooling is interrupted and the driving switched to χmax . Then the system starts to sweep the heating curve passing through L (red dashed curve). The heating is interrupted when the latter heating curve reaches the point M , which belongs to the cooling curve passing through F . The arcs IK , KL , LM , and MF build a Entropy 2022,24, 131 24 of 29 four-step bang-bang. The point K must verify TK<TJ , otherwise, it is easy to show that the four-step bang-bang cannot reach the target NESS. The above picture entails that in order to show the optimality of the two-step bangbang, we have to establish that the time needed for going from K to L , tKL , is longer than the time needed for going from J to M , tJM . Making use of Equation (8) , the equation of motion in the phase plane with χmin =0 is given by: ˙ T=−T3/2(1+3 16as 2A2),˙ A2=2T1/2 AHCS 2−A2 AHCS 2−1, (A15) where AHCS 2=aHCS 2/as 2>1. Thus, in the phase plane, one has: dT dA2 =−(1+3 16 as 2A2)T 2AHCS 2−A2 AHCS 2−1 , (A16) which can be integrated to obtain an explicit expression for T(A2) over the cooling curve [22]. From Equation (A15), we have: dt =AHCS 2−1 2 dA2 (T(A2))1/2(AHCS 2−A2), (A17) and integrating it, we obtain the expressions for tJM and tKL: tJM =AHCS 2−1 2ZA2M A2J dA2 [T(A2)]1/2(AHCS 2−A2)(A18a) tKL =AHCS 2−1 2ZA2L A2K dA0 2 [T(A0 2)]1/2(AHCS 2−A0 2)(A18b) In order to compare tJM and tKL , it is convenient to rewrite them in terms of a common variable. This variable can be naturally defined by employing the bijection that the heating curves establish between the points belonging to the JM and LK arcs. In the limit as χ→∞ , Equation (8) leads to: dT dA2 =−βT 2A2 =⇒T2Aβ 2=const. (A19) Therefore, we define a variable ξin the following way: ξ=T2(A2)Aβ 2,ξi≡T2 i≤ξ≤ξL≡T2 LAβ 2L=ξM≡T2 MAβ 2M, (A20) and we have: tJM =AHCS 2−1 2ZξL ξi dξ ξ1/4 dA2 dξ Aβ/4 2 AHCS 2−A2(ξ)(A21a) tKL =AHCS 2−1 2ZξL ξi dξ ξ1/4 dA0 2 dξ A0 2 β/4 AHCS 2−A0 2(ξ)(A21b) In view of Figure A4, it is straightforward to see that: Aβ/4 2 AHCS 2−A2(ξ)<A0 2 β/4 AHCS 2−A0 2(ξ), (A22) since A2(ξ)<A0 2(ξ),∀ξ. Now, we show that the next inequality also holds: dA2(ξ) dξ<dA0 2(ξ) dξ. (A23) Entropy 2022,24, 131 25 of 29 Taking logarithms in Equation (A20), one obtains: 2dln T(A2) dA2 +β A2dA2=dξ ξ. (A24) Plugging Equation (A16) into the last equation, we obtain: Ξ(A2)dA2=dξ ξ, (A25) where Ξ(A2)≡ −(AHCS 2−1)1+3 16 as 2A2 AHCS 2−A2 +β A2 (A26) is a decreasing function of A2 (all of the terms in the sum of its derivative are negative), positive for A2<1 and negative for A2>1. The resulting integrals can be written as: tJM =AHCS 2−1 2ZξL ξi dξ ξ5/4 1 Ξ(A2(ξ)) Aβ/4 2 AHCS 2−A2(ξ), (A27a) tKL =AHCS 2−1 2ZξL ξi dξ ξ5/4 1 Ξ(A0 2(ξ)) A0 2 β/4 AHCS 2−A0 2(ξ). (A27b) These expressions let us conclude that tJM <tKL : A2(ξ)<A0 2(ξ) and thus the integrand in tJM is smaller than the integrand in tKL . Therefore, a two-step protocol is better than a four-step one. This line of reasoning can be easily extended two six-step, eightstep, etc., protocols, showing each of them to be worse than the previous one. Protocols with an odd number of switches can be described as a limiting case of protocols with an even number of switching: For example, the three-step protocol can be seen as a four-step protocol for which Mcoincides with F. Let us inspect the case Ti> 1. Again, we start by considering a two-step protocol, comprising the arcs IJ and JF . We also build a four-step protocol, comprising the arcs IK , KL , LM , and MF . The expression for the times are formally equal to those for the case Ti>1, given by Equation (A27): tKJ =AHCS 2−1 2ZξK 1 dξ ξ5/4 A2(ξ)β/4 |Ξ(A2(ξ))|(AHCS 2−A2(ξ)), tLM =AHCS 2−1 2ZξK 1 dξ ξ5/4 A0 2(ξ)β/4 |Ξ(A0 2(ξ))|(AHCS 2−A0 2(ξ)), where 1 =ξM=ξJ and ξL=ξK . However, the inequality in Equation (A22) is reversed, so it is not direct that tKJ <tLM . Still, it is possible to inspect the behaviour of F(A2) , defined as: F(A2) = Aβ/4 2 |Ξ(A2)|(AHCS 2−A2)(A28) in the interval 1 <A2<AHCS 2 : Therein, F(A2) monotonically decreases with A2 , and thus tKJ <tLM . As a consequence, the four-step bang-bang lasts longer than the two-step bang-bang. Six-step, eight-step, etc., protocols are even worse, so the two-step bang-bang is the optimal one also for Ti>1.