Optimizing the relaxation route with optimal control

PHYSICAL REVIEW RESEARCH 3, 023128 (2021) Optimizing the relaxation route with optimal control A. Prados * Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, 41080 Sevilla, Spain (Received 26 July 2020; revised 17 January 2021; accepted 15 April 2021; published 19 May 2021) We look into the minimization of the connection time between nonequilibrium steady states. As a prototypical example of an intrinsically nonequilibrium system, a driven granular gas is considered. For time-independent driving, its natural time scale for relaxation is characterized from an empirical (the relaxation function) and a theoretical (the recently derived classical speed limits) point of view. Using control theory, we find that bangbang protocols (comprising two steps, heating with the largest possible value of the driving and cooling with zero driving) minimize the connecting time. The bang-bang time is shorter than both the empirical relaxation time and the classical speed limit: in this sense, the natural time scale for relaxation is beaten. Information theory quantities stemming from the Fisher information are also analyzed over these optimal protocols. The implementation of the bang-bang processes in numerical simulations of the dynamics of the granular gas show an excellent agreement with the theoretical predictions. Moreover, general implications of our results are discussed for a wide class of driven nonequilibrium systems. Specifically, we show that analogous bang-bang protocols, with a number of bangs equal to the number of relevant physical variables, give the minimum connecting time under quite general conditions. DOI: 10.1103/PhysRevResearch.3.023128 I. INTRODUCTION Very recent developments make it possible to define the natural time scale for the dynamical evolution or, in other words, a speed limit, in classical systems from a fundamental point of view [1–6]. In the quantum realm, speed limits have been known for a long time: the so-called MandelstamTamm [7] and Margolus-Levitin [8] bounds. A recent review on the matter is provided by Ref. [9]. Roughly speaking, the quantum speed limit entails a tradeoff between operation time and uncertainty in energy, i.e., the time-energy uncertainty relation. This idea has been extended to classical systems with Markovian dynamics: taking advantage of the similarities of the mathematical structure of the respective Hilbert spaces, the different versions of speed limits in Refs. [1–6] have been derived. Very recently, a speed limit that is the classical analog of the Mandelstam-Tamm bound has been derived [6]. It is valid for a completely general dynamics, not necessarily Markovian, and includes, as a particular case, the one derived in Ref. [5] starting from the Cramér-Rao inequality. These speed limits in Refs. [5,6] can be understood as a tradeoff between time and cost in the considered process. It must be noted, however, that their being the most restrictive bounds on operation time has not been yet proved. Currently, this is *[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. an open question for the classical speed limits, whereas for their quantum counterparts it has been rigorously established that the unification of the Mandelstam-Tamm and MargolusLevitin bounds is tight [10]. The possibility of accelerating the dynamical evolution of a given physical system has been recently analyzed in different contexts, both for classical [11–20] and quantum systems [21–26] (for a recent review, see Ref. [27]). In the classical case, the focus has been put on engineering the connection between equilibrium states for Markovian systems, the dynamics of which is described by a Fokker-Planck or a master equation. This has especially been done in the simple harmonic potential case [11–15,18–20], for which the fact that the probability distribution remains Gaussian for all times strongly simplifies the mathematical treatment.1 Here, not only do we show how to speed up the connection between nonequilibrium steady states (NESS) but also how to optimize this connection. This is done in a system that is a benchmark for out-of-equilibrium systems, a driven granular gas. In the kinetic description, neither the dynamics is Markovian (the Boltzmann-Fokker-Planck equation is nonlinear) nor the velocity distribution function is Gaussian [even in the long-time limit, when the granular gas reaches a nonequilibrium steady state (NESS)]. It must be stressed that there is no “thermodynamic” description for granular fluids. Extending thermodynamic concepts to them is far from trivial: inelastic collisions break time-reversal invariance and make the system intrinsically out 1Very recently, the connection between two nonequilibrium steady states of a Brownian gyrator has been analyzed [20], but still the probability distribution remains exactly Gaussian in that case. 2643-1564/2021/3(2)/023128(22) 023128-1 Published by the American Physical Society A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) of equilibrium, which has many, some of them unexpected, implications. For example, Shannon’s entropy no longer increases monotonically in the nondriven, freely cooling, case and there is no clear formulation of the second principle for granular fluids [28–30]. Moreover, results derived under the assumption of Markovian dynamics [1–4] are in principle not valid in the framework of kinetic theory. Nevertheless, the very recent results based on information geometry apply because the underlying dynamics is very general [5,6]. Central to the latter approach is the concept of Fisher information I(t), which is the curvature of the Kullback-Leibler divergence and is related to entropy production for Markovian dynamics [3,5,31,32]. The concept of a thermodynamic length was first introduced in the context of finite-time thermodynamics about 40 years ago [33–35], employing an approach similar to that used for defining a statistical distance in Hilbert space for quantum mechanical systems [36]. More recently, the relation between the thermodynamic length and Fisher information was unveiled [37]. Later works showed how to employ this formalism to find optimal protocols, in the sense that the relevant physical quantities attain a minimal value [38,39]. Over the last few years, further work has linked the Fisher information I(t) with the so-called thermodynamic uncertainty relations [31,32,40], also showing that I(t) is related to the entropic acceleration, i.e., to the second time derivative of Shannon’s entropy [31]. The natural time scale for connecting two NESS corresponding to different values of the driving can be characterized both empirically and theoretically. Let us consider relaxation at constant driving: at t=0, the driving is instantaneously changed from its initial to its final value. From an empirical standpoint, the relaxation time in such a process can be measured by looking for the point over the relaxation curve at which the granular temperature equals its steady value, up to a certain small precision. From a theoretical standpoint, the relaxation time is bounded from below by the classical speed limit t⩾L2/(2C)[5], where Land Care the integrals over time of √I(t) and I(t), respectively. One of the main objectives of this paper is to engineer a protocol to minimize the connection time between the two NESS. Note that the existence of nonholonomic constraints impinges on the connecting time: it is not possible to have an arbitrarily short connecting time since this leads, in general, to the violation of the constraints.2Therefore, a nonvanishing minimum connection time emerges associated to a suitable time-dependent χ(t) protocol for the driving. To work out the optimal connection, we leverage Pontryagin’s maximum principle, a key result in control theory [42,43]. Our work shows the feasibility of beating the constant driving relaxation times, both the empirical one and the theoretical speed limits, with optimal control. Specifically, the optimal process comprises two time windows: one with the largest possible value of the driving χ=χmax, and the other with no 2This is a practical shortcoming of the usual “shortcut to adiabaticity” or “engineered swift relaxation” processes. The emergence of negative values of the stiffness of the harmonic trap for too fast protocols is a well-known issue of such transformations [18,41,44]. driving at all χ=0. In the context of control theory, the kind of processes in which the control function changes abruptly between its limiting values are known as bang-bang. Here, we have two bangs because the description of our system involves two variables (see below). The order of the bangs depends on the value of the target granular temperature Tfbeing larger or smaller than the initial one. In addition, we argue that similar bang-bang processes also minimize the connecting time for a quite general class of systems. Despite the nonlinear dependence on the relevant physical variables of the evolution equations, the latter are often linear in the “control function(s).” A few illustrative examples are a colloidal particle trapped in a harmonic potential [controls: stiffness of the trap and (or) temperature of the bath [11,44]], active lattice gases [diffusion coefficient [45], noise strength and (or) density [46]], and a particle in an electric field (intensity of electric field [16]). Moreover, in most of these situations the controls (stiffness, temperature, diffusion coefficient, noise strength) are non-negative and a nonholonomic constraint arises. Bang-bang protocols thus emerge as the optimal ones because of the linearity of Pontryagin’s Hamiltonian in the control function, with the number of bangs depending on the number of independent variables. This paper is organized as follows. In Sec. II,weintroduce our model system and write the evolution equations for the granular temperature and the excess kurtosis. The characteristic relaxation times for relaxation at constant driving are analyzed in Sec. III, including the classical speed limits. Section IV is devoted to the possibility of accelerating the connection between two NESS corresponding to different values of the driving. Therein, we put forward the control problem for the minimization of the connection time and show that the optimal processes are of bang-bang type. The bang-bang processes are explicitly built in Sec. V, and the associated physical properties over them (minimum connecting time, length, and cost) are derived in Sec. VI. Numerical simulations of the dynamics are presented and compared with our analytical predictions in Sec. VII. The generality of the bang-bang protocols is investigated in Sec. VIII. We illustrate the general situation by briefly analyzing the optimal connection for a colloidal particle trapped in a three-dimensional harmonic well. Finally, Sec. IX discusses the main results of our work, their implications for a wide class of driven systems, and possible future developments. The Appendices deal with some technicalities that complement the main text. II. EVOLUTION EQUATIONS We consider a uniformly heated granular gas of ddimensional hard spheres of mass mand diameter σ, with number density n. In addition to inelastic collisions, with restitution coefficient α, the gas particles are submitted to a white-noise force of variance m2ξ2, the so-called stochastic thermostat. In the low-density limit, the dynamics of the system is accurately described by the Boltzmann-Fokker-Planck equation [47]. Our analysis is mainly done in the so-called first Sonine approximation for the kinetic equation. This approach characterizes the gas in terms of the granular temperature Tand the excess kurtosis a2. The latter incorporates non-Gaussianities 023128-2 OPTIMIZING THE RELAXATION ROUTE WITH … PHYSICAL REVIEW RESEARCH 3, 023128 (2021) in the velocity distribution function in the simplest possible way; it is the first nontrivial cumulant. The Sonine approximation accurately describes the granular gas in many different situations [47–52], and we employ it here to investigate the classical speed limits. Nevertheless, at some points of the paper we will make use of the harsher Gaussian approximation, which, as a rule of thumb, works when the property being analyzed does not vanish. Non-Gaussianities, i.e., the excess kurtosis in the Sonine approximation, only introduce corrections to the predicted behavior in such a case. We start by defining the granular temperature Tand the excess kurtosis a2: T≡mv2 d,a2≡d d+2v4 v22−1.(1) Higher-order cumulants are neglected, which makes it possible to get a closed set of equations for Tand a2. In addition, nonlinearities in a2are dropped because the typical values of the excess kurtosis are quite small. In the long-time limit, the granular gas reaches a NESS. Therein, energy loss from collisions is compensated, in average, by the energy input from the stochastic thermostat. The stationary values of the temperature and the excess kurtosis are given by [47,48] T3/2 s=mξ2 ζ01+3 16 as 2≡χ, ζ0=2nσd−1(1−α2)πd−1 2 √md(d/2) , (2a) as 2=16(1 −α)(1 −2α2) 73 +56d−24dα−105α+30(1 −α)α2.(2b) The temperature and the excess kurtosis obey the evolution equations ˙ T=ζ0χ1+3 16as 2−T3/21+3 16a2,(3a) ˙a2=2ζ0 T(T3/2−χ)a2+BT3/2as 2−a2,(3b) which are nonlinear in the temperature but linear in a2,asa consequence of the Sonine approximation. The parameter B is only a function of αand d, namely, B=aHCS 2 aHCS 2−as 2 ,(4) where aHCS 2=16(1 −α)(1 −2α2) 25 +2α2(α−1) +24d+α(8d−57),(5) is the value of the excess kurtosis in the homogeneous cooling state (HCS): the long-time time-dependent state that the system tends to approach when it cools freely, i.e., with no driving. Before proceeding further, we introduce dimensionless variables by taking adequate units for T,χ, and t: T∗=T/Ti,χ ∗=χ/T3/2 i,t∗=ζ0T1/2 it,(6) where Tiis the initial value of the temperature. Consistently, velocities are made dimensionless with √Ti/m,v∗= √m/Tiv. The excess kurtosis is already dimensionless, but for our purposes it is convenient to define the scaled variable A2≡a2/as 2.(7) We see in what follows that A2is basically non-negative, whereas a2changes sign with the inelasticity (specifically, as 2=0forα=1/√2). In the remainder of the paper, we always work with dimensionless variables, therefore, we drop the asterisks not to clutter our formulas. The corresponding evolution equations can be written as ˙ T=f1(T,A2;χ),(8a) f1(T,A2;χ)≡χ1+3 16as 2−T3/21+3 16as 2A2,(8b) ˙ A2=f2(T,A2;χ),(8c) f2(T,A2;χ)≡2 T[(T3/2−χ)A2+BT3/2(1 −A2)].(8d) Apart from a factor d,Tis basically the (dimensionless) energy per particle. Thus, the first term in f1,χ(1 +3 16 as 2), is the rate of energy input from the stochastic thermostat, while the second term, −T3/2(1 +3 16 as 2A2), is the rate of energy dissipation in collisions. Equations (8) must be supplemented with suitable initial conditions. With our choice of units, the initial temperature equals unity. Since we are interested in processes that start from the NESS corresponding to the initial temperature, T(t=0) =A2(t=0) =1.(9) III. CHARACTERISTIC RELAXATION TIME Initially, our granular fluid is in the NESS corresponding to χi=1. A typical relaxation process is constructed by suddenly changing the noise intensity from χi=1 to a different value χfat t=0. Then, the system relaxes to a new NESS with granular temperature Tfcorresponding to the noise intensity χf≡T3/2 f. Note that the stationary value of the excess kurtosis as 2is independent of the noise intensity and so is As 2, namely, As 2=1. Relaxation in this process has a certain characteristic time tR, at which the temperature has almost completely reached (complete relaxation only happens for infinite time) its steady-state value. To characterize the relaxation time from an empirical point of view, we define the relaxation function of the temperature as φ(t)=[T(t)−Tf]/(1 −Tf), such that φ(t=0) =1 and φ(t→∞)=0. The granular temperature has almost relaxed to Tfwhen φ(tR)=1, i.e., for a temperature TR(Tf,)=Tf+(1 −Tf). We consider =10−4for the sake of concreteness. This relaxation time tRcan be estimated by numerically solving the system of equations (8). Figure 1shows tRas a function of the final temperature Tffor a couple of values of (α, d), namely, (0.3,2) (circles) and (0.8,3) (open triangles). Other (α, d) pairs are not shown because all the curves would be basically superimposed. Therefore, the “natural” time scale for the relaxation of the granular temperature to its final value Tfis basically independent of αand din our dimensionless time scale defined in 023128-3 A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) FIG. 1. Characteristic relaxation time as a function of the target temperature. The numerical value of tR(symbols) is obtained by integrating the system of equations (8) numerically for the considered pair of parameters (α,d). Note that tRdepends very weakly on (α, d) and is very well predicted by the Gaussian approximation tG R(black solid line), as given by Eq. (10). Also plotted are the speed limits t(1) R (blue broken) and t(2) R(red solid) for the relaxation process, as given by Eq. (21), for d=2. Eq. (6).3It is also observed that tRis a decreasing function of the final temperature and vanishes in the limit as Tf→∞. The weak dependence of tRon (α,d) suggests that it can be quite accurately predicted by the Gaussian approximation, in which the excess kurtosis is set to zero in Eq. (8). This yields tG R(Tf,)=TR 1 dT T 3 2 f−T3 2= TR Tf−1 Tf 3√Tf ,(10) where is given by [53] (x)=ln 1+x+x2 |1−x|2−2√3arctan1+2x √3.(11) Figure 1also shows tG Ras a function of the final temperature Tf (solid line). The agreement with the numerical estimates for tR is excellent over the whole range of temperatures considered, which covers four orders of magnitude, 0.01 ⩽Tf⩽100. Equation (10) entails that tG Rvanishes algebraically in the high-temperature limit Tf1, specifically tG R∼2|ln | 3T−1/2 f,Tf1.(12) So far, we have characterized the relaxation time from an empirical standpoint. Henceforth, we consider the classical speed limits that have been recently proposed in the literature [1–6]. Specifically, we analyze those in Ref. [5] within the framework of information geometry, which are valid for a general dynamics, not necessarily Markovian. We denote the one-particle PDF for the velocity by P(v,t). The Fisher information is defined as I(t)≡dv(∂tP(v,t))2 P(v,t)=(∂tln P(v,t))2⩾0 (13) 3We have made time dimensionless with ζ0, which depends on d and is proportional to (1 −α2). and plays a central role in information geometry [54]. Therefore, the statistical length is introduced as [36,37] L=tf 0 dt I(t),(14) which represents the distance swept by the probability distribution in the time interval (0,tf). Since the probability distribution is normalized for all times, P(v,t)movesonthe unit sphere. As a result, the statistical length Lis bounded from below by the arc length between Pi(v)≡P(v,0) and Pf(v)≡P(v,tf), i.e., the so-called Bhattacharyya angle [5,36] =2 arccos dvPi(v)Pf(v),L⩾. (15) The equality L=is attained only over the geodesic in probability space, along which the Fisher information I(t) remains constant, e.g., see Appendix E of Ref. [5] for details. It has recently been proved that the Cauchy-Schwartz inequality leads to the classical speed limits tf⩾ L2 2C⩾2 2C,(16) where tfand C≡1 2tf 0 dt I(t) (17) are the operation time and the cost of the process, respectively [5]. Equation (16) expresses a tradeoff between time and cost operation, 2tfC⩾L2⩾2. The bound provided by Lis tighter but, in general, depends on the whole dynamical evolution, whereas only depends on the initial and final distributions. For our system, the speed limits above can be exactly calculated within the Gaussian approximation. Therein, I(t)= IG(t)=d/2( ˙ T/T)2,Tis a monotonic function of time, and both bounds are completely determined by Tf, as detailed in Appendix A. With the definitions γ(T)≡2√T 1+Td/2 ,ϕ(T)≡T3/2−3T1/2+2, (18) we have that4 G=2 arccos γ(Tf),Lrel G=d 2|ln Tf|,Crel G=d 4ϕ(Tf). (19) Making use of these expressions, the connecting time trel fin a relaxation process verifies the inequality trel f⩾t(1) R⩾t(2) R,(20) where t(1) R=|ln Tf|2 ϕ(Tf),t(2) R=8[arccos γ(Tf)]2 dϕ(Tf).(21) 4Not only does Lbut also Cdepend on the specific protocol followed to connect the initial and final states. This is explicitly taken into account in our notation by writing Lrel and Crel for the relaxation process. 023128-4 OPTIMIZING THE RELAXATION ROUTE WITH … PHYSICAL REVIEW RESEARCH 3, 023128 (2021) It is the geodesic in probability space that the smallest bound t(2) Rcorresponds to. A relevant question is thus the attainability of the geodesic for the granular fluid: we discuss this issue in Appendix A. Both t(1) Rand t(2) Rare shown in Fig. 1for the twodimensional case. Consistently, we have that the Gaussian estimate tG Rfor the relaxation time lies above both of them, specifically tG R/t(1) Rchanges from, approximately, 4 to 30 across the range 0.01 ⩽Tf⩽100.5Non-Gaussianities in the velocity distribution function will affect the speed limits t(1) R and t(2) R. However, we expect the smallness of a2to introduce only slight changes to the results above, as was the case of the empirical relaxation time. IV. ENGINEERED SWIFT RELAXATION Our idea is engineering a protocol, by controlling the noise intensity χ(t), that connects the initial and final NESS, the ones corresponding to χi=1 and χf=T3/2 f, in a given time tf, as short as possible. A relevant question thus arises: whether or not it is possible to beat the characteristic relaxation time of the system, not only tG Rbut also the classical speed limits t(1) Rand t(2) Rfor the relaxation process. Note that the latter is possible only for time-dependent driving. In order to connect the two NESS, the solution to Eq. (8) must verify the initial conditions (9) and also T(t=tf)=Tf,A2(t=tf)=1.(22) Therefore, Eqs. (9) and (22) constitute the boundary conditions for our engineered swift relaxation (ESR) protocol. If a solution to Eq. (8) satisfies these boundary conditions and the control function χ(t) is such that χ(t=0) =1, χ(t=tf)= T3/2 f, the system is really stationary at both the initial and final time, i.e., ˙ T(t=0) =˙ T(t=tf)=0 and ˙ A2(t=0) =˙ A2(t= tf)=0. First, we show that it is indeed possible to connect the two NESS in a finite time, by a reverse-engineering procedure.6 We start from a certain function (protocol) Tp(t) that connects the initial and final values of the temperature and, in addition, is stationary at both t=0 and tf, i.e., Tp(0) =1,Tp(tf)=Tf,˙ Tp(t=0) =˙ Tp(t=tf)=0. (23) We aim at finding a driving χp(t) and a time evolution for the scaled kurtosis A2p(t), such that (i) (Tp(t),A2p(t))is a solution to Eq. (8) for the driving χp(t), (ii) the boundary conditions for A2(t) are verified, A2p(0) =A2p(tf)=1, and (iii) the driving verifies the boundary conditions χ(0) =1, χ(tf)=T3/2 f, which ensure stationarity at both t=0 and tf. 5It should be remarked that the value of the empirical relaxation time depends on the specific value chosen for . For instance, taking =10−2instead of 10−4makes the empirical relaxation time roughly one-half of the one plotted in Fig. 1. 6The idea is similar to that employed in Refs. [11,19] for connecting two equilibrium states. Now, we employ Eqs. (8a) and (8b) to write the driving in terms of (Tp(t),A2p(t)), χp(t)= ˙ Tp(t)+[Tp(t)]3/21+3 16 as 2A2p(t) 1+3 16 as 2 .(24) Since we do not know A2p(t) yet, χp(t) is not completely determined at this point. However, insertion of Eq. (24)into(8c) and (8b) gives us a closed equation for A2p(t), which we can solve with the initial condition A2p(0) =1. Therefore, we need one free parameter, to be included in our choice for Tp(t), to “tune” A2p(t)toverifyA2p(tf)=1. Equations (23) and (24) ensure that χp(0) =1, χp(tf)=T3/2 fin such a case: the solution found in this way is indeed stationary at the initial and final times and an ESR protocol has been successfully constructed. We show how to build a simple polynomial connection in Appendix B. A. Control problem Let us consider the ESR connection problem from the following point of view. For a given (in general, time-dependent) choice of the driving intensity χ(t), the system of ordinary differential equations (ODEs) (8) predicts the corresponding time evolutions for the granular temperature Tand the excess kurtosis A2. Therefore, χ(t) plays the role of a control function. We restrict ourselves to a certain set of admissible control functions, specifically those that make it possible to connect the two NESS in a certain time tf, T(0) =1,T(tf)=Tf,A2(0) =A2(tf)=1,(25) and ensure stationarity at the initial and final times, i.e., χ(0) =1,χ(tf)=T3/2 f.(26) The control function χ(t) is assumed to be piecewise continuous in the time interval [0,tf]. The presence of finite jumps in χ(t) is not problematic from a physical point of view: already in the “basic” relaxation process χjumps from 1 to χf=T3/2 f at t=0, and Tand A2are always continuous functions of time. Above, we have shown that there exist control functions χ(t) that do the job, at least for not too short connecting times (see also Appendices A2 and B). Here, we would like to consider the problem in the light of optimal control theory: our control verifies the inequality χ(t)⩾0 and thus the possible optimization problems, such as minimizing the connection time, have a nonholonomic constraint. Therefore, we leverage Pontryagin’s maximum principle [42,43]tosolve the optimisation problem and find the optimal control χ(t)for the corresponding physical situation. For the sake of mathematical rigor, we also consider that the noise intensity is bounded from above, χ(t)⩽χmax; afterwards, we will take the limit χmax →∞. B. Optimizing the connection Let us consider the following optimization problem: we want to obtain the minimum time for making the connection between the two NESS, i.e., we want to minimize tf=tf 0dt. In order to apply Pontryagin’s procedure, we define a variable 023128-5 A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) y0(t) such that the optimization problem is equivalent to the minimization of y0(tf), i.e., ˙y0=f0(T,A2,χ),f0(T,A2,χ)≡1,y0(tf)=tf.(27) To proceed, we introduce Pontryagin’s Hamiltonian (y,ψ,χ)≡ψ0f0(y,χ)+ψ1f1(y,χ)+ψ2f2(y,χ).(28) In this context, we employ the notation y1≡T,y2≡A2, y≡(y0,y1,y2), ψ≡(ψ0,ψ 1,ψ 2) to simplify some formulas. The variables yiand their conjugate momenta ψi,i=0,1,2, evolve following the Hamiltonian system ˙yi=∂ ∂ψi ,˙ ψi=−∂ ∂yi=− 2  j=0 ψj ∂fj ∂yi .(29) From the construction above, the functions fjdo not depend on y0and thus ˙ ψ0=0, ψ0=const. Pontryagin’s maximum principle states necessary conditions for optimal connection: in order that (χ∗(t),T∗(t),A∗ 2(t))be optimal, it is necessary that there exists a nonzero continuous vector function ψ∗(t)= (ψ∗ 0(t),ψ∗ 1(t),ψ∗ 2(t))corresponding to (χ∗(t),T∗(t),A∗ 2(t)) such that for all t,0⩽t⩽tf, (i) the canonical system (29) holds, (ii) if we define the supremum of as function of the control, H(y,ψ)=supχ(y,ψ,χ), we have that H(y∗(t),ψ∗(t))=(y∗(t),ψ∗(t),χ∗(t)),(30) and (iii) the two constants of motion ψ∗ 0and H∗≡ H(y∗(t),ψ∗(t))satisfy ψ∗ 0⩽0 and H∗=0. To find the supremum of with respect to χ, we calculate ∂/∂χ: either χ∗follows from the condition ∂/∂χ|χ∗=0 or lies at the boundaries of the interval [0,χ max]. Making use of Eqs. (8b), (8d), (27), and (28), we obtain ∂ ∂χ =ψ11+3 16as 2−2ψ2 A2 T,(31) which does not depend on χand thus does not allow for finding χ∗. This is a consequence of being a linear function of χand therefore either χ∗=0orχmax, depending on the sign of ∂/∂χ. The optimal control jumps from χ∗=0to χmax at those times for which ∂/∂χ changes from negative to positive, and vice versa. This kind of discontinuous optimal controls are commonly known as bang-bang [16,21,26,43]. The simplest situation is thus a two-step bang-bang process, with two possibilities: (i) high driving window χ∗(t)= χmax,0<t⩽tJ, followed by free cooling χ∗(t)=0, tJ< t<tf, and (ii) first free cooling χ∗(t)=0, 0 <t⩽tJ,followed by high driving χ∗(t)=χmax,tJ<t<tf. From our study of the polynomial connection, we may guess that (i) is the optimal protocol for Tf>1, but this ansatz has to be checked. V. BANG-BANG OPTIMAL CONTROLS In this section we carry out an in-depth study of the bangbang controls we have just described above. For the sake of simplicity, we explicitly build such protocols for the case χmax 1.7 In general, we focus on the motion of point describing the state of the system in the phase-space plane (A2,T): Eq. (8) is a system of first-order ODEs and trajectories in the phase-space plane cannot intersect. Making use of them, we arrive at 2 T dT dA2=χ1+3 16 as 2−T3/21+3 16 as 2A2 (T3/2−χ)A2+BT3/2(1 −A2)(32) A. Heating-cooling bang-bang Here, we analyze the bang-bang process in which the granular fluid is first heated, χ(t)=χmax 1, 0 <t⩽tJ, and afterwards freely cools, χ(t)=0, tJ<t<tf. Taking the limit χmax 1inEq.(32) and solving the resulting separable ODE with initial condition (A2i,Ti)=(1,1) in the (A2,T) plane, we get T2A1+3 16 as 2 2=T2 iA1+3 16 as 2 2,i=1,0⩽t⩽tJ.(33) Now we investigate the behavior of the system in the second time window tJ⩽t⩽tf. Setting χ=0inEq.(32) and taking into account Eq. (4), we arrive again at a separable first-order ODE, the solution of which is given by 2ln T Tf=3 16as 2AHCS 2−1(A2−1) +1+3 16as 2AHCS 2 ×AHCS 2−1ln AHCS 2−A2 AHCS 2−1, tJ⩽t⩽tf.(34) For the final time, t=tf, we have that T=Tfand A2=A2f = 1. We obtain a relation between TJand Tfby particularizing Eq. (34) for the joining time t=tJ, specifically 2ln Tf TJ=3 16as 2AHCS 2−1(1 −A2J) +1+3 16as 2AHCS 2AHCS 2−1 ×ln AHCS 2−1 AHCS 2−A2J.(35) In turn, TJand A2Jare related by T2 JA1+3 16 as 2 2J=1,(36) as implied by Eq. (33). As a consequence, Eq. (35)givesa one-to-one relation between Tfand TJor Tfand A2J.8 7The existence of the restriction χ⩽χmax is of practical nature, we assume that the intensity of the heat bath cannot be arbitrarily large, whereas the constraint χ⩾0 is of fundamental nature: in average, the stochastic forcing always increases the kinetic energy of the particles. 8In the first part of the bang-bang process, the system heats with χmax and thus ˙ T⩾0andTJ⩾1, which entails that A2J⩽1. In the 023128-6 OPTIMIZING THE RELAXATION ROUTE WITH … PHYSICAL REVIEW RESEARCH 3, 023128 (2021) 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 10 FIG. 2. Bang-bang protocol for Tf>1. A representative example of the motion of the granular gas in the (A2,T) plane is shown: specifically, we have considered α=0.8andd=2. Other values of (α,d) lead to a completely analogous picture. The bang-bang process connects the initial NESS with (A2i =1,Ti=1) with the final state (A2f =1,Tf>1) and comprises two parts: first heating (red dashed line) followed by cooling (solid blue lines). Different target points (A2f =1,Tf>1) over the vertical line A2=1 (dotted) are reached by starting the cooling part from different points (A2J,TJ) over the heating curve. A qualitative plot of the motion of the system in the (A2,T) plane is shown in Fig. 2. In the first part of the protocol, 0<t⩽tJ, the system is heated with χ(t)=χmax and follows Eq. (33). The second part of the bang-bang process starts at a given point (A2J,TJ) over this line. Therefore, for tJ<t<tf, the system freely cools with χ(t)=0 and thus follows Eq. (34). This part of the bang-bang finishes when the system hits the vertical line A2=1 at the corresponding target point (A2f =1,Tf). In order to keep the system stationary for t⩽0 and t⩾tf, the control function has sudden jumps at these points: χ(t⩽0) =1, χ(t)=χmax 1for0<t⩽tJ, χ(t)=0fortJ<t<tf,χ(t⩾tf)=T3/2 f. Note that with this order of the bangs, the bang-bang protocol always leads the system to a final NESS with Tf>1. The impossibility of reaching Tf<1 can be physically understood in the following way: in the first part of the bang-bang process, the system always heats, TJ>1, and the corresponding excess kurtosis decreases in absolute value: the velocity distribution function becomes closer to a Gaussian, A2J<1. Therefore, the initial slope, i.e., at the point (A2J,TJ), of the curve for the second part of the bang-bang process (blue solid in Fig. 2)is always larger than the slope of the curve for the first part (red dashed) at the same point. This can be shown by inspecting the corresponding expressions for dT/dA2and taking into account that A2J<1. Since evolution curves corresponding to different initial points cannot intersect in the (A2,T) plane, it must be concluded that Tf>1. To reach NESS with Tf<1, one intuitively thinks that inverting the bangs, i.e., first cooling and afterwards heating should be necessary. We prove this is indeed the case in the next section. limit as TJ→∞, the velocity distribution becomes Gaussian at the joining time, A2J→0. 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 0.0 0.2 0.4 0.6 0.8 1.0 FIG. 3. Bang-bang protocol for Tf<1. As a representative example we show the case (α=0.3,d=2), the qualitative picture of the motion of the possible in the (A2,T) plane is the same for other values of (α,d). The bang-bang process connects the initial NESS (A2i =1,Ti=1) with the target NESS (A2f =1,Tf<1). Again, it comprises two parts, but the order of the bangs is reversed, as compared with Fig. 2: first the system is cooled (blue solid line) and afterwards is heated (red dashed lines). Different target NESS over the vertical line A2=1 (dotted) are reached by starting the heating part from different points (A2J,TJ) over the cooling curve. B. Cooling-heating bang-bang Next, we look into the bang-bang protocol in which the granular fluid freely cools first, χ(t)=0, 0 <t⩽tJ, and afterwards is strongly heated, χ(t)=χmax,tJ<t<tf.The same separable first-order ODEs in the (A2,T) plane have to be solved, but with different initial conditions. In the cooling stage, the resulting evolution is 2lnT=3 16as 2AHCS 2−1(A2−1) +1+3 16as 2AHCS 2AHCS 2−1ln AHCS 2−A2 AHCS 2−1, 0⩽t⩽tJ.(37) For tJ<t<tf, the system evolves with χ(t)=χmax 1, and we have that T2A1+3 16 as 2 2=T2 f,tJ⩽t⩽tf.(38) This equation is similar to (33), but here the second part of the protocol ends at the point (A2f =1,Tf). Since it starts from at t=tJfrom the point (A2J,TJ), we get the relation T2 f=T2 JA1+3 16 as 2 2J.(39) In turn, TJand A2Jare related by the particularization of Eq. (37)fort=tJ: 2lnTJ=3 16as 2AHCS 2−1(A2J−1) +1+3 16as 2AHCS 2AHCS 2−1ln AHCS 2−A2J AHCS 2−1. (40) Figure 3shows the motion of the system in the (A2,T) plane for this bang-bang process. Therefore, it is analogous to 023128-7 A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) Fig. 2, but with the order of the bangs reversed. In the first part of the bang-bang process, the system follows the curve given by Eq. (37). In its second part, starting from a given point (A2J,TJ) over this line, the system evolves according to Eq. (38). This bang-bang process connects the initial NESS (A2i =1,Ti=1) with the final NESS (A2f =1,Tf), but now we have that Tf⩽1. In order to keep the system stationary for t=0 and tf, the control function has again sudden jumps at the initial and final times: at t=0+, it changes from 1 to 0; at t=t− f, it changes from χmax to χf=T3/2 f. VI. PHYSICAL PROPERTIES FOR THE BANG-BANG OPTIMAL CONTROLS Let us analyze in more detail the just described bangbang protocols, which drive the system from the initial NESS (A2i =1,Ti=1) to the final NESS (A2f =1,Tf= 1). The two-step bang-bang processes provide us with the minimum connecting time, and we obtain it as a function of Tfboth for Tf>1 and for Tf<1.9In addition, we calculate the statistical length and the cost for them. In the following, we investigate the cases Tf>1 and Tf<1 separately. A. Heating-cooling bang-bang: Tf>1 We start by analyzing the heating-cooling bang-bang process described in Sec. VA, which makes it possible to reach temperatures that are larger than the initial one Tf>1. It comprises two steps: (i) χ=χmax for 0 <t⩽tJ, and (ii) χ(t)=0fortJ<t<tf. 1. Minimum connecting time Along the first part of the heating-cooling bang-bang, i.e., in the time window 0 <t⩽tJ,wehave ˙ T∼χmax(1 +3 16 as 2). Therefore, we get tJ=TJ−1 χmax1+3 16 as 2→0,χ max →∞.(41) Note that tJ→0, but χmaxtJremains finite. In the second part of the process tJ<t<tf, the system freely cools with χ=0. Therefore, making use of Eq. (8) and taking into account that tJ→0, tf=TJ Tf dT T3/21+3 16 as 2A2(T),(42) in which A2(T) is implicitly given by Eq. (34): it is thus impossible to carry out this integral analytically, at least in an exact manner. 9In linear response, when |Tf−1|1, the suboptimality of bangbang protocols with more than two steps is a consequence of a theorem in the number of switchings (see for instance Theorem 10 in Sec. III.17 of Pontryagin’s book [42]). The formal proof for this specific nonlinear case is quite lengthy and will be published elsewhere. A physical argument for the number of steps of the optimal bang-bang for a general nonlinear case with nvariables is provided in Sec. VIII. 1 2 3 4 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 1234 0 0.5 1 1.5 FIG. 4. Minimum connection time as a function of the target temperature, for Tf>1. All lines correspond to d=2, and different values of the restitution coefficient are considered. From top to bottom, α=0.3 (black dashed), α=0.8 (magenta dotted), α=0.9 (green dotted-dashed), α=0.98 (orange dashed), and α=0.998 (brown dotted). Note that t(0) fvanish in the limit as Tf→1inall cases, whereas its high-temperature behavior depends on the inelasticity. For reference, the speed limits for relaxation t(1) Rand t(2) Rare plotted in the inset, which lie well above t(0) f. We can obtain an approximate analytical expression for the connecting time if we bring to bear that 3as 2/16 is quite small over the whole range of restitution coefficient 0 ⩽α⩽1, and A2is expected to be of the order of unity. Accordingly, denoting by t(0) fthe connecting time obtained by putting as 2=0in Eq. (42),10 we get t(0) f=T(0) J Tf dT T3/2=2T−1/2 f−T(0) J−1/2.(43) Above, T(0) Jmeans that TJmust be consistently put in terms of Tfby considering Eqs. (35) and (36)foras 2=0butA2=O(1), which yields Tf=T(0) JAHCS 2−1 AHCS 2−A(0) 2J AHCS 2−1 2 ,T(0) J=A(0) 2J−1/2.(44) Equations (43) and (44) provide us with the connecting time t(0) fas a function of the final temperature Tf, both of them are given in terms of A(0) 2J,0<A(0) 2J<1. Figure 4shows t(0) f, as given by Eqs. (43) and (44), as a function of the target temperature Tf. Over the scale of the figure, t(0) fis indistinguishable from the numerical integration of Eq. (42). Specifically, we have plotted the curves for the two-dimensional case and several values of the inelasticity. The minimum connection time given by control theory clearly beats the speed limits for relaxation t(1,2) R, given by Eq. (21), which are shown in the inset. It is observed that t(0) fdecreases 10Note that we are keeping A2and thus this is not equivalent to the Gaussian approximation, in which a2is completely disregarded from the very beginning. Therein, the optimal connecting time vanishes because only one bang with χmax suffices to reach the final temperature and tG f=(Tf−1)/χmax →0. 023128-8 OPTIMIZING THE RELAXATION ROUTE WITH … PHYSICAL REVIEW RESEARCH 3, 023128 (2021) 1.0 1.5 2.0 2.5 3.0 10 50 100 500 1000 FIG. 5. Acceleration factor as a function of the final temperature, for Tf>1. Different lines correspond to different values of the restitution coefficient α,ford=2, with the same code as in Fig. 4.Note the logarithmic scale in the vertical axis. as the restitution coefficient αincreases, vanishing in the elastic limit as α→1. Physically, this can be understood as follows: the system does not cool in the second part of the process for α→1. Thus, T(0) J→Tfand t(0) f→0. Mathematically, AHCS 2→1 in the elastic limit, which ensures that T(0) J→Tf. Asymptotic expressions for t(0) fcan be derived in some limits. First, in the high-temperature limit, T(0) Jbecomes large and A(0) 2Jsmall; therefore, we have that t(0) f∼2T−1/2 f⎡ ⎣1−AHCS 2−1 AHCS 2 AHCS 2−1 4⎤ ⎦,Tf1.(45) Note that the right-hand side vanishes in the elastic limit, in which AHCS 2→1. Second, we consider the linear response limit Tf−11. Therein, Eq. (44) implies that T(0) J−1∼ (Tf−1)1/2and then t(0) fvanishes as t(0) f∼AHCS 2−1 AHCS 21/2 (Tf−1)1/2,Tf−11.(46) Again, the factor AHCS 2−1 makes the right-hand side vanish in the elastic limit. As already commented above, the minimum value of the connecting time t(0) fbeats the speed limits in Eq. (21). Therefore, it entails a really large acceleration of the relaxation, as compared with the characteristic relaxation time tG Rgiven by Eq. (10). We can measure the acceleration factor in the bang-bang process by the ratio tG R/t(0) f.InFig.5,weplotthis ratio as a function of the target temperature for d=2 and the same values of the restitution coefficient as in Fig. 4. Specifically, relaxation is speeded up by more than one of order of magnitude for high temperatures and by a diverging amount as the final temperature approaches unity, i.e., in the linear response limit. For high target temperatures, tG R/t(0) fgoes to a constant value that depends on the inelasticity: both times vanish as T−1/2 f[see Eqs. (12) and (45)]. 2. Associated length and cost It is worth investigating the length Ltraversedbythe system in probability space and the cost Cof the bang-bang process. There is a tradeoff between operation and cost, as expressed by the “thermodynamic uncertainty relation” (16). As a result, it is expected that minimizing the operation time, as we have done, should entail a neat separation from the geodesic in probability space, for which L=, and an increase in the cost C. In Appendix C, we show that the Fisher information is given by I(t)=I(0)(t)1+Oas 2,I(0)(t)=d 2˙ T(t) T(t)2 .(47) In the following, we calculate the lowest-order contribution to the length and cost, i.e., L(0) =tf 0 dt I(0)(t)=d 2tf 0 dt  ˙ T(t) T(t) (48) and C(0) =1 2tf 0 dt I(0)(t)=d 4tf 0 dt˙ T(t) T(t)2 .(49) First, we look into the length swept by the probability distribution. Taking into account that TJ>Tf>1, we have to split the integral in Eq. (48) into two summands: the first one for the time interval [0,tJ], in which the temperature monotonically increases from Ti=1toTJ, and the second one for [tJ,tf], in which the temperature monotonically decreases from TJto Tf. Therefore, we have that L(0) =d 2tJ 0 dt ˙ T T−tf tJ dt ˙ T T =d 2TJ 1 dT T+TJ Tf dT T=d 2ln T2 J Tf.(50) We plot our estimate L(0) for the length over the optimal bang-bang connection as a function of the target temperature in Fig. 6. Also plotted are the length for the relaxation process Lrel G(blue broken) and the length over the geodesic G(solid red), given in Eq. (19). Despite the heating-cooling bang-bang minimizes the connection time, which is much shorter than the relaxation time tG R, the corresponding length is always larger than Lrel Gsince L(0) −Lrel G=√2dln TJ Tf⩾0,(51) with the equality holding when TJ=Tf. It is inelasticity that makes TJdifferent from Tfand thus increases the length swept by the system in probability space. Only in the elastic limit α→1wehavethatL(0) →Lrel G, as neatly observed in the plot. Second, we consider the cost of this heating-cooling bangbang process. Since the bang-bang process minimizes the connection time, a large value of the cost associated with the speed limit is expected. Similarly to what we have just done for the length, the smallness of non-Gaussianities allows us to estimate the cost with C(0). Again, by splitting the time 023128-9 A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) contrast with the vanishing of the connecting time for the granular case. The latter behavior stems from the collision rate being proportional to √T, which also accelerates the cooling step of the bang-bang in the granular case, whereas the relaxation rates of the harmonic modes κβare independent of the temperature. In case III, the incorporation of the third mode increases the connection time, t(III) f>t(II) f, but keeps the qualitative picture unchanged. IX. DISCUSSION In this work, we have applied information geometry and control theory ideas to a system described at the kinetic level. The resulting framework is physically appealing. On the one hand, information geometry put lower bounds on the operating time: the classical speed limits in Refs. [5,6] apply although the dynamics is not Markovian. On the other hand, control theory makes it possible to build protocols that entail large accelerations of the system dynamics, by minimizing the connection time. It is known that reverse engineering techniques, such as engineered swift equilibration [11–14,20], can connect states in times that are shorter than the empirical relaxation time. Here, we show that optimal control protocols not only are able to beat the empirical relaxation time for relaxation (by more than an order of magnitude) but also the recently derived classical speed limits, which are considerably shorter than the empirical time. The latter play in classical systems a role similar to that of the quantum speed limits, associated with the time-energy uncertainty relation, in quantum systems. This beating of the classical speed limit for relaxation does not represent a contradiction since the optimal control protocols involve a time-dependent driving. There appears a clear asymmetry between the cases Tf>1 and Tf<1: recall that in our dimensionless units the initial temperature equals unity. For the case Tf>1, the optimal connecting times are rather small, vanishing in the limits Tf→1 and Tf→∞. The smallness of the minimum connecting times can be understood in a physical way: in the Gaussian approximation, the minimum connecting time vanishes because the optimal protocol is clearly a pulse of very high noise intensity such that the granular temperature instantaneously changes from 1 to Tf. Therefore, it is non-Gaussianities, specifically, the excess kurtosis a2that is small, that make impossible to instantaneously connect the two NESS for Tf>1. The excess kurtosis decreases and therefore the state after the instantaneous heating pulse is not stationary. For the case Tf<1, the minimum connecting times are longer than those for heating. Again, this can be understood from the Gaussian approximation: therein, the optimal protocol is letting the system freely cool, i.e., with driving intensity χ=0. At difference with the case Tf>1, the minimum connection time for Tf<1 does not vanish because free cooling involves a finite time. Interestingly, both for Tf>1 and Tf<1 non-Gaussianities make the connecting times longer: this is physically understood by taking into account that nonGaussianities stem from the inelasticity of collisions. One of the main results of our paper is the emergence of optimal bang-bang protocols for minimizing the connecting time. In our case, the bang-bang processes comprise two steps (i.e., one switching): (i) instantaneous heating with a very high driving intensity χmax →∞and (ii) free cooling, i.e., no driving, χ=0. The order of the bangs is different for Tf>1 and Tf<1: heating-cooling for Tf>1, but cooling-heating for Tf<1. Qualitatively, this can be understood as follows: in both cases, the first step corresponds to what would be done in the Gaussian approximation. However, the existence of non-Gaussianities entail that the excess kurtosis does not have the stationary value at the end of the first step. This imbalance is somehow mended by the second step of the bang-bang. Indeed, bang-bang processes are expected to emerge as the optimal protocols, in the sense of minimizing the connection time, in a wide variety of physical situations, not only for the specific case of granular fluids. The evolution equations of the relevant physical properties typically include “control functions”: other quantities, the time dependence of which can be externally controlled. Often, the control function λ (stiffness of a harmonic trap, temperature of the bath, diffusion coefficient, noise strength, density, etc., depending on the physical context) verify that (i) the evolution equations of the physical properties are linear in them, and (ii) a nonholonomic constraint limits their physically acceptable values, λmin ⩽λ⩽λmax. Examples abound, from the trapped Brownian particle [11,44] or active systems [45,46] to a particle moving in an electric field [16]. Therein, the mathematical structure of Pontryagin’s principle ensures that the optimal controls minimizing the connecting time are of bang-bang type. Let us consider thus a physical system such that, at the macroscopic (or hydrodynamic, or thermodynamic, ...) level of description is described by nphysical quantities.14 The number nis thus small, in the examples above n= 1[11,16,44]orn=2[45,46], as is the case of the granular fluid in the Sonine approximation. A relevant question is as follows: How many bangs, i.e., how many time windows inside which either λ(t)=λmin or λ(t)=λmax, are necessary to make the optimal connection? For n=1, a simple bang, either λ(t)=λmin or λ(t)=λmax, of duration tfmakes it possible to adjust the final value of the single variable. For n=2, there will be a mismatch between the target value of the second variable and the one obtained with a simple bang that tunes the final value of the first variable. This makes it necessary to introduce a second step with the control being switched to the opposite limit: two bangs, either λ(t)=λmax followed by λ(t)=λmin or vice versa, of durations t1and t2, with tf= t1+t2, allow for matching the final values of two variables, and thus suffices for n=2. In general, we have to introduce n−1 jumps at intermediate times to allow for matching all n variables: the number of bangs equals the number of variables. Specifically, we have shown that the above picture is indeed the correct one by solving a simple but relevant physical situation: the compression and decompression of a Brownian particle trapped in a d-dimensional harmonic 14Also at the mesoscopic level of description (fluctuating hydrodynamics, stochastic thermodynamics, ...), which incorporates fluctuations of these quantities to the picture. 023128-16 OPTIMIZING THE RELAXATION ROUTE WITH … PHYSICAL REVIEW RESEARCH 3, 023128 (2021) potential, by controlling the temperature of the bath.15 For n=2 (axial symmetry for d=3 or 2), a two-step bang-bang process, qualitatively similar to that found in the granular fluid, carries out the optimal connection: heating-cooling or cooling-heating, depending on the final value of the temperature being larger or smaller than the initial one. For d=3, an optimal three-step bang-bang process arises: heating-coolingheating or cooling-heating-cooling, again depending on the final value of the temperature. Our work has been focused on the minimization of the connecting time tfbetween two NESS of the granular fluid. Not only is the optimization of the connection time between two given states a relevant problem from a fundamental point of view, but also has potential applications in different contexts. For example, minimizing the connection time in the adiabatic, in the sense of zero heat, branches is essential for building a finite-time version of Carnot’s heat engine that maximizes the delivered power [56]. Also, the optimization of the relaxation route to equilibrium or to a NESS is of interest in connection with behaviors such as the Mpemba effect, which is currently a very active field of research [52,57–61]. In addition, we have considered the associated statistical length Land cost Cover the optimal processes in the granular fluid. The length of the optimal bang-bang protocols is always longer that that of the relaxation process, which contributes to increase the bound for the connecting time; recall that tf⩾ L2/(2C). However, this is compensated by the cost, which is dominated by a term proportional to the maximum value of the noise intensity χmax →∞, i.e., the cost diverges. Therefore, the bound goes to zero and the connecting time may (and we have shown that this is indeed the case here) beat the speed limit for relaxation. Had we minimized the cost, we would have obtained an infinite operation time. Therein, the system would be for all times in the NESS corresponding to the instantaneous value of the noise intensity and thus the cost would vanish. The divergence of the operation time, when minimizing the cost, is the counterpart of the divergence of the cost, when minimizing the connection time. But the analogy ends there. We have shown in the granular gas that the minimum connecting time does not vanish despite the diverging cost: the cooling part of the bang-bang protocol involves a finite time. Our general arguments above, and the specific example of the d-dimensional oscillator, show that this will be also the case in many physical situations where a nonholonomic constraint is present, e.g., when controls are non-negative. Our approach opens interesting perspectives for further research. In the context of granular systems, it is far from trivial to rigorously prove the global stability of the long-time NESS. Indeed, there are strong signs, but not a formal proof, that it is the relative Kullback-Leibler divergence with respect to the stationary distribution, and not Shannon’s entropy, that acts as a Lyapunov functional [29,30,62,63]. In this sense, the 15Optical confinement makes it possible to control the time dependence of the effective temperature seen by the Brownian particle by randomly shaking the confining trap or by using Brownian particles with an inherent charge submitted to a random electric field [55]. role of the Fisher information for rigorously establishing the H theorem for granular gases is worth investigating. For Markovian dynamics, the cost Chas been related in general to entropy production [3,5,31,32]. In the realm of kinetic theory and, more specifically, for the granular case, the situation is far more complicated. Even admitting Shannon’s as the good definition of entropy for the granular case, there is not a clear-cut way of splitting entropy production into “irreversible” and “flux” contributions, as discussed in Ref. [28]. Therefore, elucidating the physical meaning of information geometry’s cost (beyond stating that it is the physical quantity appearing in the thermodynamic uncertainty relation for the connecting time) in granular fluids is a relevant problem that remains to be solved. Kinetic theory tools are not restricted to low density (or moderate density if using Enskog’s equation instead of Boltzmann’s) gases, either molecular or granular. They have also been successfully applied to other intrinsically nonequilibrium systems such as active matter [64–69]. Also, the classical kinetic approach holds for dilute ultracold gases: despite the very low temperatures, they are still far from the threshold at which quantum effects cease to be negligible [70–72]. Therefore, it is worth looking into the application of information-geometry concepts and, in general, the extension of our results to these physical contexts. ACKNOWLEDGMENTS I acknowledge the financial support from project PGC2018-093998-B-I00, funded by: FEDER/Ministerio de Ciencia e Innovación–Agencia Estatal de Investigación (Spain). Also, I thank C. A. Plata, A. Santos, and E. Trizac for discussions and their critical reading of the manuscript. APPENDIX A: CLASSICAL SPEED LIMITS AND GEODESIC PROPERTIES IN THE GAUSSIAN APPROXIMATION In this Appendix, we derive explicit expressions for the speed limits for the relaxation process and analyze some properties of the geodesic in probability space. The analysis is carried out in the Gaussian approximation, where the granular gas is completely described by the granular temperature T(t). The velocity distribution function is assumed to be the d-dimensional Maxwellian PG(v;T(t))=(2πT)−d/2exp − v2 2T.(A1) The temperature obeys the evolution equation ˙ T=χ−T3/2,(A2) which stems from Eqs. (8a) and (8b) with as 2=0. Throughout, we employ the subindex Gto denote those quantities calculated within the Gaussian approximation. 1. Classical speed limits for the relaxation process First, we calculate the Fisher information. Its definition (13) entails that IG(t)=[∂tln PG(v,T(t))]2G, where ...Gmeans average with the Gaussian distribution in 023128-17 A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) Eq. (A1). Making use of ∂tln PG(v;T(t))= ˙ T(t) 2T(t)−d+ v2 T,(A3) and the fact that v4=(d+2)v22/d=d(d+2)T2for a Gaussian distribution, it is readily shown that IG(t)=d 2˙ T(t) T(t)2 G .(A4) The subindex Gon the right-hand side means that we have to consistently evaluate T(t) within the Gaussian approximation, i.e., over the solution of Eq. (A2). Second, the Bhattacharyya angle is also obtained from its definition, Eq. (15). Specifically, we calculate the angle between the Gaussian distributions corresponding to the initial temperature (recall that Ti=1 with our choice of units) and the final one Tf. Taking into account that (i) the integrand is Gaussian and (ii) the d-dimensional integral factorizes into the product of didentical integrals, we have that dvPG(v,Ti=1)PG(v,Tf)=2√Tf 1+Tfd/2 (A5) and G=2 arccos 2√Tf 1+Tfd/2.(A6) We can also derive analytical expressions for the statistical length and the cost in the Gaussian case. In particular, we are interested here in the relaxation process between the initial and final NESS, with time-independent driving χ= χf=T3/2 f. In the Gaussian approximation, the evolution of the temperature is monotonic and therefore we can integrate over the temperature instead of over time. For the statistical length, we get Lrel G=d 2tf 0 dt ˙ T(t) T(t)=d 2tf 0 dt ˙ T(t) T(t) =d 2Tf 1 dT T=d 2|ln Tf|,(A7) whereas for the cost we have that Crel G=d 4tf 0 dt˙ T(t) T(t)2 =d 4Tf 1 dT ˙ T T2 =d 4Tf 1 dT T3/2 f−T3/2 T2=d 4T3/2 f−3T1/2 f+2. (A8) Let trel fbe the time for connecting the initial and final NESS in the relaxation process. The classical speed limits derived in Ref. [5] ensure that, within the Gaussian approximation we are employing, trel f⩾ L2 G 2CG ⩾2 G 2CG ,(A9) i.e., tf⩾|ln Tf|2 2T3/2 f−3T1/2 f+2 ⩾ 8 arccos22√Tf 1+Tfd/2 dT3/2 f−3T1/2 f+2.(A10) The above inequalities are equivalent to those in Eqs. (18)– (21)ofthemaintext. 2. Geodesic in the Gaussian approximation The geodesic in probability space can be further characterized. To do so, it is useful to introduce a rescaled time τ=t/tffor a process of duration tf, with 0 ⩽τ⩽1. The probability distribution over the geodesic P∗(v,τ) is obtained by minimizing the statistical length Lwith the constraint dvP(v,τ)=1, ∀τ. A straightforward but rather lengthy calculation leads to the result [5] P∗(v,τ)=√Pi(v)sin 2(1 −τ)+√Pf(v)sin 2τ sin  2. (A11) Over the geodesic, the cost is also minimized but it depends on the connecting time tf, i.e., on the parametrization of the geodesic [37]. Specifically, it is obtained that C∗=2/(2tf). It is worth stressing that Eq. (A11) is exact and general, valid for any dynamics, as shown in Ref. [5]. Let us analyze the geodesic in more detail, for the specific case of the granular fluid. The granular temperature over the geodesic is directly obtained by taking the second moment of the probability distribution in Eq. (A11): T∗(τ)=1 sin2G 2sin2G 2(1 −τ)+Tfsin2G 2τ +4Tf 1+Tf cos G 2sin G 2τ ×sin G 2(1 −τ).(A12) The first and second terms on the right-hand side come from Pi(v) and Pf(v), respectively, and are exact. The third term stems from the product √Pi(v)Pf(v) and has been written in the Gaussian approximation. Consistently, we have substituted with G, which is given as a function of Tfby Eq. (A6), in the three terms. Of course, T∗(τ)>0 because all terms on the right-hand side of Eq. (A12) are non-negative. A relevant issue is whether it is possible for the granular gas to move over the geodesic or not. We can answer this question in the Gaussian approximation we are employing. The temperature program over the geodesic follows from the time-dependent protocol for the driving χ∗ G(t)=˙ T∗(t)+[T∗(t)]3/2,(A13) making use of Eq. (A2). In the scaled time τ, the driving is thus χ∗ G(τ)=1 tf dT∗(τ) dτ+[T∗(τ)]3/2.(A14) The second term corresponds to the “quasistatic” driving: in the NESS, χ=T3/2. The first term is the finite-time contribution, which vanishes in the limit as tf→∞. 023128-18 OPTIMIZING THE RELAXATION ROUTE WITH … PHYSICAL REVIEW RESEARCH 3, 023128 (2021) Taking the derivative of Eq. (A12), after a little bit of algebra one gets dT∗ dτ=G 2sin 2G 2 Tf−1 Tf+1{Tfsin(Gτ)+sin[G(1 −τ)]}. (A15) The temperature evolution is monotonic over the geodesic: the sign of the derivative is completely encoded in Tf−1 because all the other terms are strictly positive. This introduces an asymmetry between the cases Tf>1 and Tf<1. For Tf>1, the finite-time contribution t−1 fdT∗(τ)/dτis always positive and there is a well-defined driving that makes the temperature sweep the geodesic curve, even in the limit tf→0+.For Tf<1, the finite-time contribution t−1 fdT∗(τ)/dτis negative: for short enough connecting time tf, it will become larger (in absolute value) than the quasistatic driving and make χG<0. The above discussion implies that the geodesic cannot be swept for too short connecting times for Tf<1. It is illustrative to consider the particularization of Eq. (A15)forτ=1 2to give an estimate for the connecting time such that χGbecomes negative: dT∗ dττ=1/2=G 2sin G 2(Tf−1).(A16) The condition for having χG(τ=1/2) <0is tf<G 2sin G 2 1−Tf [T∗(τ=1/2)]3/2,(A17) which can be ensured if tf<G 2sin G 2(1 −Tf).(A18) Note the time for which χGfirst becomes negative is longer than the right-hand side of Eq. (A18). APPENDIX B: SIMPLE ESR POLYNOMIAL CONNECTION Here we discuss how the ESR protocol is built from a simple polynomial. We need at least a fourth-order polynomial with five coefficients: four to adjust the boundary conditions (23) for the temperature, and one extra parameter to impose that A2p(tf)=1. To start with, it is adequate to employ the scaled time τ= t/tfintroduced in Appendix A2and to work with the thermal velocity vth ≡√T. Consistently, vth,p(τ)=√T p(τ), and we rewrite Eq. (24)as χp(τ)= v2 th,p(τ)2 tf dvth,p(τ) dτ+vth,p(τ)1+3 16 as 2A2p(τ) 1+3 16 as 2 . (B1) Insertion of this expression for the noise intensity into the evolution equation of the excess kurtosis gives, after some algebra, dA2p(τ) dτ=− 4 1+3 16 as 2 dln vth,p(τ) dτA2p(τ) −2tfB+3as 2A2p(τ) 16 +3as 2vth,p(τ)(A2p(τ)−1). (B2) (a) 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.5 2.0 2.5 (b) 0.0 0.2 0.4 0.6 0.8 1.0 −10 0 10 20 30 40 (c) 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.5 2.0 2.5 (d) 0.0 0.2 0.4 0.6 0.8 1.0 −10 0 10 20 30 40 50 FIG. 14. Thermal velocity vth,pand noise intensity χpfor the fourth-order polynomial connection, as a function of the normalized time τ=t/tf. All panels are for the two-dimensional case: (a) and (b) correspond to α=0.8 and (c) and (d) correspond to α=0.3. In each panel, three curves are plotted for different connection times: from bottom to top, tf=1 (solid black), tf=0.5 (dashed purple), and tf=0.25 (dotted green). For the shortest connection time, χp(t) becomes negative inside a certain time window. We solve this equation, with the initial condition A2p(0) =1, with the following fourth-order polynomial for the thermal velocity, vth,p(τ)=1+cτ2+(4vth −2c)τ3+(c−3vth)τ4, (B3) where vth ≡√T=√Tf−1. The parameter cis tuned to meet the boundary condition A2p(tf)=1: there is only one fourth-order polynomial making the connection. We have carried out the above procedure by numerically solving Eq. (B2) for the two-dimensional case, i.e., hard disks. We show the numerical results for Tf>1inFig.14, specifically for Tf=4(vth =1). The following qualitative behavior is observed: as the connecting time tfin decreased, the driving χp(t) goes to very high values before decreasing to lower, even negative, values. Evidently, the noise intensity χp(t) cannot become negative, so this means that the ESR connection cannot be done with a fourth-order polynomial for too short times.16 The observed behavior hints at the emergence of a minimum, nonvanishing, value of the connecting time for the ESR protocol, both for Tf>1 and Tf<1. This feeling is reinforced by employing higher-order polynomials. For example, in the fifth-order case, there is a monoparametric family of polynomials connecting the initial and final NESS. Nevertheless, χp(t) becomes negative for tfbelow a certain value, over the whole family of polynomials making the connection. 16A similar behavior is found for Tf<1, but the driving first decreases, taking also negative values for short enough tf,andafterwards increases to overshoot its final value. 023128-19 A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) APPENDIX C: FISHER INFORMATION IN THE SONINE APPROXIMATION Here, we look into the Fisher information I(t) within the Sonine approximation we have considered throughout. The velocity distribution function is expanded as P(v,t)=PG(v,T(t))1+∞  k=2 ak(t)Skv2 2T,(C1) where PG(v,T(t))is the Maxwellian distribution of Eq. (A1), ak(t) are the coefficients of the expansion, which are related to the cumulants, and finally Sk(x)≡L(d−2 2) k(x) are the associated Laguerre (or Sonine) polynomials [73]. The explicit expression for the first Sonine polynomials are [47] S0(x)=1,S1(x)=−x+d 2,(C2a) S2(x)=1 2x2−d+2 2x+d(d+2) 8.(C2b) The first Sonine approximation consists in keeping only the first term in the expansion, with coefficient a2that equals the excess kurtosis, and neglecting nonlinear contributions in a2in all the expressions derived from Eq. (C1).17 For our purposes, it is convenient to rewrite Eq. (C1) as follows: we introduce a dimensionless velocity c(v,T(t))=v/2T(t),(C3) and the order of unity quantity A2(t) defined in Eq. (7), so that P(v,t)=e−c2 [2πT(t)]d/21+as 2A2(t)S2(c2).(C4) Now we proceed to calculate the Fisher information. For that, we take into account ∂tf(c2)=df (c2) d(c2)∂tc2=− ˙ T(t) T(t)c2df(c2) d(c2)(C5) to write ∂tln P(v,t)=− ˙ T TS1(c2)+as 2˙ A2S2(c2)−A2 ˙ T Tc2dS2(c2) d(c2). (C6) In order to obtain I(t), Eq. (C6) is squared and averaged with the probability distribution (C4), neglecting nonlinear terms in as 2. After a little algebra, one gets I(t)=IG(t)+as 2A2˙ T T2 ×S2 1(c2)S2(c2)+2c2S1(c2)dS2(c2) d(c2),(C7) 17The first polynomial S1does not appear in the expansion because the Gaussian distribution gives the correct value for the temperature T(t), i.e., the corresponding coefficient a1(t) vanishes identically. where we have omitted the time dependence of T(t) and A2(t) to simplify the notation, and defined f(c)≡dcf(c)φ(c),φ(c)=π−d/2e−c2(C8) as the average of f(c) with the dimensionless Gaussian distribution φ(c). The averages in Eq. (C7) are thus ddimensionless integrals of polynomials with the Gaussian distribution, which result in I(t)=I(0)(t)1−d+2 2as 2A2(t),I(0)(t)=d 2˙ T(t) T(t)2 . (C9) There is no contribution coming from the term proportional to ˙ A2in Eq. (C6) because of the orthogonality of Sonine polynomials Sj(c2)Sk(c2)=0for j= k. Also, note that I(0)(t)= IG(t) because we no longer set the excess kurtosis to zero in the first Sonine approximation. Notwithstanding, the smallness of as 2implies that the main contribution to the Fisher information comes from I(0)(t). APPENDIX D: NORMAL MODES FOR THE d-DIMENSIONAL HARMONIC POTENTIAL Our starting point is the Fokker-Planck equation (65), for the harmonic potential case. The transformation to normal modes is orthogonal, i.e., there exists an orthogonal matrix of elements Cjk such that xj= d  β=1 Cjβξβ, d  β=1 CjβCkβ=δjk,(D1) which diagonalizes the symmetric matrix, with elements λjk, of the harmonic well Uh(x), d  β=1 d  β=1 CjβλjkCkβ=κβδβ,β,(D2) Uh(ξ)=1 2 d  β=1 κβξ2 β.(D3) Therefore, the Fokker-Planck equation can be rewritten in terms of the probability P(ξ,t)=P(x,t)as γ∂ tP(ξ,t)=∇ξ·[∇ξUh(ξ)P(ξ,t)]+kBT(t)∇2 ξP(ξ,t). (D4) At equilibrium, the initial distribution P(ξ,t=0) factorizes into the product of dGaussian distributions with zero mean and variances Ti/κβ, one for each normal mode. Since ∂ξβUh(ξ)=κβξβ, the joint distribution P(ξ,t) still factorizes into dGaussian distributions with zero mean for all times. Therefore, it is completely characterized by the variances of the modes ξ2 β, which obey the uncoupled equations γd dt ξ2 β=−2κβξ2 β+2kBT(t).(D5) It is convenient to go to dimensionless variables, by introducing suitable units for time, length, and temperature. We 023128-20 OPTIMIZING THE RELAXATION ROUTE WITH … PHYSICAL REVIEW RESEARCH 3, 023128 (2021) label the modes in such a way that κ1⩽···⩽κd. We define t∗=κ1 γt,ξ ∗ β=ξβ √kBTi/κ1 ,T∗(t)=T(t) Ti .(D6) With our choice of units, (ξ∗ 1)2i=1 and T∗ i=1. We can rewrite (D5)as d dt∗(ξ∗ β)2=−2κ∗ β(ξ∗ β)2+2T∗(t),(D7) where κ∗ β=κβ/κ1, i.e., κ∗ 1=1⩽···⩽κ∗ d.(D8) Equation (D7) is equivalent to (66) of the main text. Therein, we have dropped the asterisks in order not to clutter our formulas. [1] B. Shanahan, A. Chenu, N. Margolus, and A. del Campo, Quantum Speed Limits across the Quantum-to-Classical Transition, Phys. Rev. Lett. 120, 070401 (2018). [2] M. Okuyama and M. Ohzeki, Quantum Speed Limit is Not Quantum, Phys.Rev.Lett.120, 070402 (2018). [3] S. Ito, Stochastic Thermodynamic Interpretation of Information Geometry, Phys.Rev.Lett.121, 030605 (2018). [4] N. Shiraishi, K. Funo, and K. Saito, Speed Limit for Classical Stochastic Processes, Phys.Rev.Lett.121, 070601 (2018). [5] S. Ito and A. Dechant, Stochastic time-evolution, information geometry and the Cramer-Rao Bound, Phys. Rev. X 10, 021056 (2020). [6] S. B. Nicholson, L. P. García-Pintos, A. del Campo, and J. R. Green, Time-information uncertainty relations in thermodynamics, Nat. Phys. 16, 1211 (2020). [7] L. Mandelstam and I. Tamm, The Uncertainty Relation Between Energy and Time in Non-relativistic Quantum Mechanics, in Selected Papers,editedbyB.M.Bolotovskii,V.Y. Frenkel, and R. Peierls (Springer, Berlin, 1991), pp. 115–123. [8] N. Margolus and L. B. Levitin, The maximum speed of dynamical evolution, Phys. D (Amsterdam) 120, 188 (1998). [9] S. Deffner and S. Campbell, Quantum speed limits: From Heisenberg’s uncertainty principle to optimal quantum control, J. Phys. A: Math. Theor. 50, 453001 (2017). [10] L. B. Levitin and T. Toffoli, Fundamental Limit on the Rate of Quantum Dynamics: The Unified Bound Is Tight, Phys. Rev. Lett. 103, 160502 (2009). [11] I. A. Martínez, A. Petrosyan, D. Guéry-Odelin, E. Trizac, and S. Ciliberto, Engineered swift equilibration of a Brownian particle, Nat. Phys. 12, 843 (2016). [12] P. Muratore-Ginanneschi and K. Schwieger, An application of pontryagin’s principle to brownian particle engineered equilibration, Entropy 19, 379 (2017). [13] G. Li, H. T. Quan, and Z. C. Tu, Shortcuts to isothermality and nonequilibrium work relations, Phys. Rev. E 96, 012144 (2017). [14] M. Chupeau, S. Ciliberto, D. Guéry-Odelin, and E. Trizac, Engineered swift equilibration for Brownian objects: From underdamped to overdamped dynamics, New J. Phys. 20, 075003 (2018). [15] J. A. C. Albay, S. R. Wulaningrum, C. Kwon, P.-Y. Lai, and Y. Jun, Thermodynamic cost of a shortcuts-to-isothermal transport of a Brownian particle, Phys. Rev. Res. 1, 033122 (2019). [16] V. Martikyan, D. Guéry-Odelin, and D. Sugny, Comparison between optimal control and shortcut to adiabaticity protocols in a linear control system, Phys.Rev.A101, 013423 (2020). [17] K. Funo, N. Lambert, F. Nori, and C. Flindt, Shortcuts to Adiabatic Pumping in Classical Stochastic Systems, Phys. Rev. Lett. 124, 150603 (2020). [18] J. A. C. Albay, P.-Y. Lai, and Y. Jun, Realization of finite-rate isothermal compression and expansion using optical feedback trap, Appl. Phys. Lett. 116, 103706 (2020). [19] C. A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, Finitetime adiabatic processes: Derivation and speed limit, Phys. Rev. E101, 032129 (2020). [20] A. Baldassarri, A. Puglisi, and L. Sesta, Engineered swift equilibration of a Brownian gyrator, Phys. Rev. E 102, 030105(R) (2020). [21] X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. Guéry-Odelin, and J. G. Muga, Fast Optimal Frictionless Atom Cooling in Harmonic Traps: Shortcut to Adiabaticity, Phys. Rev. Lett. 104, 063002 (2010). [22] X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Shortcut to Adiabatic Passage in Twoand Three-Level Atoms, Phys. Rev. Lett. 105, 123003 (2010). [23] S. Deffner and E. Lutz, Quantum Speed Limit for NonMarkovian Dynamics, Phys.Rev.Lett.111, 010402 (2013). [24] S. Campbell and S. Deffner, Trade-Off Between Speed and Cost in Shortcuts to Adiabaticity, Phys. Rev. Lett. 118, 100601 (2017). [25] T.-N. Xu, J. Li, T. Busch, X. Chen, and T. Fogarty, Effects of coherence on quantum speed limits and shortcuts to adiabaticity in many-particle systems, Phys. Rev. Res. 2, 023125 (2020). [26] Y. Ding, T.-Y. Huang, K. Paul, M. Hao, and X. Chen, Smooth bang-bang shortcuts to adiabaticity for atomic transport in a moving harmonic trap, Phys. Rev. A 101, 063410 (2020). [27] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications, Rev. Mod. Phys. 91, 045001 (2019). [28] I. Bena, F. Coppex, M. Droz, P. Visco, E. Trizac, and F. van Wijland, Stationary state of a heated granular gas: Fate of the usual H-functional, Phys. A (Amsterdam) 370, 179 (2006). [29] U. M. B. Marconi, A. Puglisi, and A. Vulpiani, About an Htheorem for systems with non-conservative interactions, J. Stat. Mech. (2013) P08003. [30] M. I. García de Soria, P. Maynar, S. Mischler, C. Mouhot, T. Rey, and E. Trizac, Towards an H-theorem for granular gases, J. Stat. Mech. (2015) P11009. [31] S. B. Nicholson, A. del Campo, and J. R. Green, Nonequilibrium uncertainty principle from information geometry, Phys. Rev. E 98, 032106 (2018). [32] Y. Hasegawa and T. Van Vu, Uncertainty relations in stochastic processes: An information inequality approach, Phys.Rev.E 99, 062126 (2019). 023128-21 A. PRADOS PHYSICAL REVIEW RESEARCH 3, 023128 (2021) [33] P. Salamon and R. S. Berry, Thermodynamic Length and Dissipated Availability, Phys. Rev. Lett. 51, 1127 (1983). [34] P. Salamon, J. D. Nulton, and R. S. Berry, Length in statistical thermodynamics, J. Chem. Phys. 82, 2433 (1985). [35] T. Feldmann, B. Andresen, A. Qi, and P. Salamon, Thermodynamic lengths and intrinsic time scales in molecular relaxation, J. Chem. Phys. 83, 5849 (1985). [36] W. K. Wootters, Statistical distance and Hilbert space, Phys. Rev. D 23, 357 (1981). [37] G. E. Crooks, Measuring Thermodynamic Length, Phys. Rev. Lett. 99, 100602 (2007). [38] D. A. Sivak and G. E. Crooks, Thermodynamic Metrics and Optimal Paths, Phys. Rev. Lett. 108, 190602 (2012). [39] E. J. Kim, U. Lee, J. Heseltine, and R. Hollerbach, Geometric structure and geodesic in a solvable model of nonequilibrium process, Phys. Rev. E 93, 062127 (2016). [40] A. Dechant, Multidimensional thermodynamic uncertainty relations, J. Phys. A: Math. Theor. 52, 035001 (2019). [41] C. A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, Optimal work in a harmonic trap with bounded stiffness, Phys. Rev. E 99, 012140 (2019). [42] L. S. Pontryagin, Mathematical Theory of Optimal Processes (CRC Press, Boca Raton, FL, 1987). [43] D. Liberzon, Calculus of Variations and Optimal Control Theory: A Concise Introduction (Princeton University Press, Princeton, NJ, 2012). [44] M. Chupeau, B. Besga, D. Guéry-Odelin, E. Trizac, A. Petrosyan, and S. Ciliberto, Thermal bath engineering for swift equilibration, Phys. Rev. E 98, 010104(R) (2018). [45] M. Kourbane-Houssene, C. Erignoux, T. Bodineau, and J. Tailleur, Exact Hydrodynamic Description of Active Lattice Gases, Phys.Rev.Lett.120, 268003 (2018). [46] A. Manacorda and A. Puglisi, Lattice Model to Derive the Fluctuating Hydrodynamics of Active Particles with Inertia, Phys. Rev. Lett. 119, 208003 (2017). [47] T. P. C. Van Noije and M. H. Ernst, Velocity distributions in homogeneous granular fluids: The free and the heated case, Granul. Matter 1, 57 (1998). [48] J. M. Montanero and A. Santos, Computer simulation of uniformly heated granular fluids, Granular Matter 2,53 (2000). [49] M. I. García de Soria, P. Maynar, and E. Trizac, Universal reference state in a driven homogeneous granular gas, Phys. Rev. E 85, 051301 (2012). [50] E. Trizac and A. Prados, Memory effect in uniformly heated granular gases, Phys. Rev. E 90, 012204 (2014). [51] A. Prados and E. Trizac, Kovacs-Like Memory Effect in Driven Granular Gases, Phys. Rev. Lett. 112, 198001 (2014). [52] A. Lasanta, F. Vega Reyes, A. Prados, and A. Santos, When the Hotter Cools More Quickly: Mpemba Effect in Granular Fluids, Phys. Rev. Lett. 119, 148001 (2017). [53] T. P. C. van Noije, M. H. Ernst, E. Trizac, and I. Pagonabarraga, Randomly driven granular fluids: Large-scale structure, Phys. Rev. E 59, 4326 (1999). [54] S.-i. Amari, Information Geometry and Its Applications,Applied Mathematical Sciences Vol. 194 (Springer, Tokyo, 2016). [55] I. A. Martínez, E. Roldán, J. M. R. Parrondo, and D. Petrov, Effective heating to several thousand kelvins of an optically trapped sphere in a liquid, Phys. Rev. E 87, 032159 (2013). [56] C. A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, Building an irreversible Carnot-like heat engine with an overdamped harmonic oscillator, J. Stat. Mech. (2020) 093207. [57] Z. Lu and O. Raz, Nonequilibrium thermodynamics of the Markovian Mpemba effect and its inverse, Proc. Natl. Acad. Sci. USA 114, 5083 (2017). [58] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvión, A. Gordillo-Guerrero, D. Iñiguez, A. Lasanta, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane et al., The Mpemba effect in spin glasses is a persistent memory effect, Proc. Natl. Acad. Sci. USA 116, 15350 (2019). [59] A. Gal and O. Raz, Precooling Strategy Allows Exponentially Faster Heating, Phys. Rev. Lett. 124, 060602 (2020). [60] A. Kumar and J. Bechhoefer, Exponentially faster cooling in a colloidal system, Nature (London) 584, 64 (2020). [61] A. Lapolla and A. Godec, Faster Uphill Relaxation in Thermodynamically Equidistant Temperature Quenches, Phys. Rev. Lett. 125, 110602 (2020). [62] C. A. Plata and A. Prados, Global stability and H-theorem in lattice models with nonconservative interactions, Phys. Rev. E 95, 052121 (2017). [63] A. Megías and A. Santos, Kullback-Leibler Divergence of a Freely Cooling Granular Gas, Entropy 22, 1308 (2020). [64] A. Baskaran and M. C. Marchetti, Enhanced Diffusion and Ordering of Self-Propelled Rods, Phys. Rev. Lett. 101, 268101 (2008). [65] A. Baskaran and M. Cristina Marchetti, Nonequilibrium statistical mechanics of self-propelled hard rods, J. Stat. Mech. (2010) P04019. [66] T. Ihle, Kinetic theory of flocking: Derivation of hydrodynamic equations, Phys.Rev.E83, 030901(R) (2011). [67] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143 (2013). [68] T. Ihle, Chapman-Enskog expansion for the Vicsek model of self-propelled particles, J. Stat. Mech. (2016) 083205. [69] L. L. Bonilla and C. Trenado, Contrarian compulsions produce exotic time-dependent flocking of active particles, Phys. Rev. E 99, 012612 (2019). [70] D. S. Lobser, A. E. S. Barentine, E. A. Cornell, and H. J. Lewandowski, Observation of a persistent non-equilibrium state in cold atoms, Nat. Phys. 11, 1009 (2015). [71] D. Guéry-Odelin and E. Trizac, Ultracold atoms: Boltzmann avenged, Nat. Phys. 11, 988 (2015). [72] M. Hohmann, F. Kindermann, T. Lausch, D. Mayer, F. Schmidt, E. Lutz, and A. Widera, Individual Tracer Atoms in an Ultracold Dilute Gas, Phys. Rev. Lett. 118, 263401 (2017). [73] M. Abramowitz, I. A. Stegun, and R. H. Romer, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Am.J.Phys.56, 958 (1988). 023128-22