Optimising the relaxation route with optimal control

We look into the minimisation of the connection time between non-equilibrium steady states. As a prototypical example of an intrinsically non-equilibrium system, a driven granular gas is considered. For time-independent driving, its natural time scale for relaxation is characterised from an empirical -- the relaxation function -- and a theoretical -- the recently derived classical speed limits -- point of view. Using control theory, we find that bang-bang protocols -- comprising two steps, heating with the largest possible value of the driving and cooling with zero driving -- minimise 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 analysed 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 non-equilibrium 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.


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][2][3][4][5][6]. In the quantum realm, speed limits have been known for a long time: the socalled Mandelstam-Tamm [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 trade-off 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][2][3][4][5][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 trade-off 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 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 Margolus-Levitin bounds is tight [10].
Here, not only do we show how to speed up the connection between non-equilibrium steady states (NESS) but also how to optimise 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 non-linear-nor the velocity distribution function is Gaussian-even in the long-time limit, when the granular gas reaches a non-equilibrium steady state (NESS).
It must be stressed that there is no "thermodynamic" description for granular fluids. Extending thermodynamics concepts to them is far from trivial: inelastic collisions break time reversal invariance and make the system intrinsically out-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 [29][30][31].
Moreover, results derived under the assumption of Markovian dynamics [1][2][3][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,32,33].
The concept of a thermodynamic length was first in-troduced in the context of finite-time thermodynamics about forty years ago [34][35][36], employing an approach similar to that used for defining a statistical distance in Hilbert space for quantum mechanical systems [37]. More recently, the relation between the thermodynamic length and Fisher information was unveiled [38]. Later works showed how to employ this formalism to find optimal protocols, in the sense that the relevant physical quantities attain a minimal value [39,40]. Over the last few years, further work has linked the Fisher information I(t) with the so-called thermodynamic uncertainty relations [32,33,41], also showing that I(t) is related to the entropic acceleration, i.e. to the second time derivative of Shannon's entropy [32]. The natural time scale for connecting two NESS corresponding to different values of the driving can be characterised 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 ≥ L 2 /(2C) [5], where L and C are 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 minimise 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 [42]. Therefore, a non-vanishing minimum connection time emerges associated to a suitable timedependent χ(t) protocol for the driving. To work out the optimal connection, we leverage Pontryagin's maximum principle, a key result in control theory [43,44].
Our work shows the feasibility of beating theconstant 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 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 T f being larger or smaller than the initial one.
In addition, we argue that similar bang-bang processes also minimise the connecting time for a quite general class of systems. Despite the non-linear 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,45]), active lattice gases (diffusion coefficient [46], noise strength and (or) density [47]), 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 manuscript is organised as follows. In Sec. II, we introduce our model system and write down the evolution equations for the granular temperature and the excess kurtosis. The characteristic relaxation times for relaxation at constant driving are analysed 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 minimisation 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 analysing 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 m and 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 m 2 ξ 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 [48].
Our analysis is mainly done in the so-called first Sonine approximation for the kinetic equation. This approach characterises the gas in terms of the granular temperature T and the excess kurtosis a 2 . The latter incorporates non-Gaussianities in the velocity distribution function in the simplest possible way, it is the first non-trivial cumulant. The Sonine approximation accurately describes the granular gas in many different situations [48][49][50][51][52][53], 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 analysed does not vanish. Non-Gaussianities, i.e. the excess kurtosis in the Sonine approximation, only introduce corrections to the predicted behaviour in such a case.
We start by defining the granular temperature T and the excess kurtosis a 2 , Higher-order cumulants are neglected, which makes it possible to get a closed set of equations for T and a 2 .
In addition, non-linearities in a 2 are 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 value of the temperature and the excess kurtosis are given by [48,49] The temperature and the excess kurtosis obey the evolution equationṡ which are non-linear in the temperature but linear in a 2 , as a consequence of the Sonine approximation. The parameter B is only a function of α and d, namely where is the value of the excess kurtosis in the homogeneous cooling state-the long-time time-dependent state that the system tends to approach when cools freely, i.e. with no driving. Before proceeding further, we introduce dimensionless variables by taking adequate units for T , χ, and t: where T i is the initial value of the temperature. Consistently, velocities are made dimensionless with T i /m, v * = m/T i v. The excess kurtosis is already dimensionless, but for our purposes it is convenient to define the scaled variable We see in what follows that A 2 is basically non-negative, whereas a 2 changes sign with the inelasticity (specifically, a s 2 = 0 for α = 1/ √ 2). In the remainder of the paper, we always work with dimensionless variables-therefore, we drop the asterisks not to clutter our formulae. The corresponding evolution equations can be written aṡ Apart from a factor d, T is basically the (dimensionless) energy per particle. Therefore, the first term in f 1 , χ(1+ 3 16 a s 2 ) is the rate of energy input from the stochastic thermostat, while the second term, −T 3/2 (1+ 3 16 a s 2 A 2 ), 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, 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 χ f at t = 0. Then, the system relaxes to a new NESS with granular temperature T f corresponding to the noise intensity χ f ≡ T 3/2 f . Note that the stationary value of the excess kurtosis a s 2 is independent of the noise intensity and so is A s 2 , namely A s 2 = 1. Relaxation in this process has a certain characteristic time t R , at which the temperature has almost completely reachedcomplete relaxation only happens for infinite time-its steady state value.
To characterise the relaxation time from an empirical point of view, we define the relaxation function of the temperature as φ(t) = (T (t) − T f )/(1 − T f ), such that φ(t = 0) = 1 and φ(t → ∞) = 0. The granular temperature has almost relaxed to T f when φ(t R ) = 1, i.e. for a temperature T R (T f , ) = T f + (1 − T f ). We consider = 10 −4 for the sake of concreteness. This relaxation time t R can be estimated by numerically solving the system of equations (8). Figure 1 shows t R as a function of the final temperature T f for 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 T f is basically independent of α and d in our dimensionless time scale defined in Eq. (6) [54]. It is also observed that t R is a decreasing function of the final temperature and vanishes in the limit as T f → ∞.
The weak dependence of t R on (α, 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 where Ω is given by [55] Ω(x) = ln Figure 1 also shows t G R as a function of the final temperature T f (solid line). The agreement with the numerical estimates for t R is excellent over the whole range of temperatures considered, which covers four orders of magnitude, 0.01 ≤ T f ≤ 100. Eq. (10) entails that t G R vanishes algebraically in the high temperature limit T f 1, specifically So far, we have characterised the relaxation time from an empirical standpoint. Henceforth, we consider the classical speed limits that have been recently proposed in the literature [1][2][3][4][5][6]. Specifically, we analyse 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 and plays a central role in information geometry [56]. Therefrom, the statistical length is introduced as [37, 38] which represents the distance swept by the probability distribution in the time interval (0, t f ). Since the probability distribution is normalised for all times, P (v, t) moves on the unit sphere. As a result, the statistical   length L is bounded from below by the arc-length between P i (v) ≡ P (v, 0) and P f (v) ≡ P (v, t f ), i.e. the so-called Bhattacharyya angle [5,37] 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 where t f and are the operation time and the cost of the process, respectively [5]. Eq. (16) expresses a trade-off between time and cost operation, 2t f C ≥ L 2 ≥ Λ 2 . The bound provided by L is 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) = I G (t) = d/2(Ṫ /T ) 2 , T is a monotonic function of time, and both bounds are completely determined by T f -as detailed in Appendix A. With the definitions we have that [57] Making use of these expressions, the connecting time t rel f in a relaxation process verifies the inequality where t (1) It is the geodesic in probability space that the smallest bound t (2) R corresponds 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) R and t (2) R are shown in Fig. 1 for the two-dimensional case.
Consistently, we have that the Gaussian estimate t G R for the relaxation time lies above both of them, specifically t G R /t R changes fromapproximately-4 to 30 across the range 0.01 ≤ T f ≤ 100 [58]. Non-Gaussianities in the velocity distribution function will affect the speed limits t  R . However, we expect the smallness of a 2 to 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 = T 3/2 f -in a given time t f , as short as possible. A relevant question thus arises: whether or not it is possible to beat the characteristic relaxation time of the systemnot only t G R but also the classical speed limits t R and t (2) R for 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 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 , the system is really stationary at both the initial and final time, i.e.Ṫ (t = 0) =Ṫ (t = t f ) = 0 andȦ 2 (t = 0) =Ȧ 2 (t = t f ) = 0. First, we show that it is indeed possible to connect the two NESS in a finite time, by a reverse-engineering procedure [59]. We start from a certain function (protocol) T p (t) that connects the initial and final values of the temperature and, in addition, is stationary at both t = 0 and t = t f , i.e.
(23) We aim at finding a driving χ p (t) and a time evolution for the scaled kurtosis A 2p (t), such that (i) (T p (t), A 2p (t)) is a solution to Eq. (8) for the driving χ p (t), (ii) the boundary conditions for A 2 (t) are verified, A 2p (0) = A 2p (t f ) = 1, and (iii) the driving verifies the boundary conditions Now, we employ Eqs. (8a) and (8b) to write the driving in terms of (T p (t), A 2p (t)), Since we do not know A 2p (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 A 2p (t), which we can solve with the initial condition A 2p (0) = 1. Therefore, we need one free parameterto be included in our choice for T p (t)-to "tune" A 2p (t) to verify A 2p (t f ) = 1. Eqs. (23) and (24) ensure that in 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. The control problem
Let us consider the ESR connection problem from the following point of view. For a given-in general timedependent-choice of the driving intensity χ(t), the system of ODEs (8) predicts the corresponding time evolutions for the granular temperature T and the excess kurtosis A 2 . 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 t f , and ensure stationarity at the initial and final times, i.e.
The control function χ(t) is assumed to be piecewise continuous in the time interval [0, t f ]. 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 = T 3/2 f at t = 0, and T and A 2 are 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 A 2 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 optimisation problems, such as minimising the connection time, have a non-holonomic constraint. Therefore, we leverage Pontryagin's maximum principle [43,44] to solve the optimisation problem and find the optimal control χ(t) for the corresponding physical situation. For the sake of mathematical rigour, we also consider that the noise intensity is bounded from above, χ(t) ≤ χ max ; afterwards we will take the limit χ max → ∞.

B. Optimising the connection
Let us consider the following optimisation problem: we want to obtain the minimum time for making the connection between the two NESS, i.e. we want to minimise t f = t f 0 dt. In order to apply Pontryagin's procedure, we define a variable y 0 (t) such that the optimisation problem is equivalent to the minimisation of y 0 (t f ), i.e.
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 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 χ * = 0 or χ * = χ max , depending on the sign of ∂Π/∂χ. The optimal control jumps from χ * = 0 to χ * = χ 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,44]. The simplest situation is thus a two-step bang-bang process, with two possibilities: (i) high driving window From our study of the polynomial connection, we may guess that (i) is the optimal protocol for T f > 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 bang-bang controls we have just described above. For the sake of simplicity, we explicitly build such protocols for the case χ max 1 [60]. In general, we focus on the motion of point describing the state of the system in the phase space plane (A 2 , 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 A. Heating-cooling bang-bang Here, we analyse the bang-bang process in which the granular fluid is first heated, χ(t) = χ max 1, 0 ≤ t ≤ t J , and afterwards freely cools, χ(t) = 0, t J ≤ t ≤ t f . Taking the limit χ max 1 in Eq. (32) and solving the resulting separable ODE with initial condition (A 2,i , T i ) = (1, 1) in the (A 2 , T ) plane, we get Now we investigate the behaviour of the system in the second time window, t J ≤ t ≤ t f . Putting χ = 0 in Eq. (32) and taking into account Eq. (4), we arrive again at a separable first order ODE, the solution of which is given by For the final time, t = t f , we have that T = T f and A 2 = A 2f = 1. We obtain a relation between T J and T f by particularising Eq. (34) for the joining time t = t J , specifically In turn, T J and A 2J are related by as implied by Eq. (33). As a consequence, Eq. (35) gives a one to one relation between T f and T J -or T f and A 2J [61]. A qualitative plot of the motion of the system in the (A 2 , T ) plane is shown in Fig. 2. In the first part of the protocol, 0 ≤ t ≤ t J , the system is heated with χ(t) = χ max and follows Eq. (33). The second part of the bangbang process starts at a given point (A 2J , T J ) over this line. Therefrom, for t J < t < t f , 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 A 2 = 1 at the corresponding target point (A 2f = 1, T f ). In order to keep the system stationary for t ≤ 0 and t ≥ t f , the control function has sudden jumps at these points, Note that with this order of the bangs, the bang-bang protocol always leads the system to a final NESS with T f > 1. The impossibility of reaching T f < 1 can be physically understood in the following way: in the first part of the bang-bang process, the system always heats, T J > 1, and the corresponding excess kurtosis decreases in absolute value-the velocity distribution function becomes closer to a Gaussian, A 2J < 1. Therefore, the initial slope, i.e. at the point (A 2J , T J ), 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 /dA 2 and taking into account that A 2J < 1. Since evolution curves corresponding to differential initial points cannot intersect in the (A 2 , T ) plane, it must be concluded that T f > 1. To reach NESS with T f < 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.

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 ≤ t J , and afterwards is strongly heated, The same separable first-order ODEs in the (A 2 , T ) plane have to be solved, but with different initial conditions. In the cooling stage, the resulting evolution is For t J ≤ t ≤ t f , the system evolves with χ(t) = χ max 1, and we have that This equation is similar to (33), but here the second part of the protocol ends at the point (A 2f = 1, T f ). Since it starts from at t = t J from the point (A 2J , T J ), we get the relation In turn, T J and A 2J are related by the particularisation of Eq. (37) for t = t J , Figure 3 shows the motion of the system in the (A 2 , T ) plane for this bang-bang process. Therefore, it is analogous to 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), for 0 ≤ t ≤ t J . In its second part, starting from a given point (A 2J , T J ) over this line, the system evolves according to Eq. (38). This bang-bang process connects the initial NESS (A 2i = 1, T i = 1) with the final NESS (A 2f = 1, T f ), but now we have that T f ≤ 1. In order to keep the system stationary for t = 0 and t = t f , 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

VI. PHYSICAL PROPERTIES FOR THE BANG-BANG OPTIMAL CONTROLS
Let us analyse in more detail the just described bangbang protocols, which drive the system from the initial NESS (A 2i = 1, T i = 1) to the final NESS (A 2f = 1, T f = 1). The two-step bang-bang processes provide us with the minimum connecting time, and we obtain it as a function of T f both for T f > 1 and for T f < 1 [62]. In addition, we calculate the statistical length and the cost for them.
In the following, we investigate the cases T f > 1 and T f < 1 separately.
A. Heating-cooling bang-bang: We start by analysing the heating-cooling bang-bang process described in Sec. V A, which makes it possible to reach temperatures that are larger than the initial one, T f > 1. It comprises two steps:

Minimum connecting time
Along the first part of the heating-cooling bang-bang, i.e. in the time window 0 ≤ t ≤ t J , we haveṪ ∼ χ max 1 + 3 16 a s 2 . Therefore, we get Note that t J → 0, but χ max t J remains finite.
In the second part of the process, t J ≤ t ≤ t f , the system freely cools with χ = 0. Therefore, making use of Eq. (8) and taking into account that t J → 0, in which A 2 (T ) is implicitly given by Eq. (34): it is thus impossible to carry out this integral analytically, at least in an exact manner.
We can obtain an approximate analytical expression for the connecting time if we bring to bear that 3a s 2 /16 is quite small over the whole range of restitution coefficient, 0 ≤ α ≤ 1, and A 2 is expected to be of the order of unity. Accordingly, denoting by t Above, T   large and A (0) 2J small; therefore we have that (45) Note that the rhs vanishes in the elastic limit, in which A HCS 2 → 1. Second, we consider the linear response limit, Again, the factor A HCS 2 − 1 makes the rhs vanish in the elastic limit.
As already commented above, the minimum value of the connecting time t (0) f beats the speed limits in Eq. (21). Therefore, it entails a really large acceleration of the relaxation, as compared with the characteristic relaxation time t G R given by Eq. (10). We can measure the acceleration factor in the bang-bang process by the ratio t G R /t (0) f . In Fig. 5, we plot this 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, t G R /t

Associated length and cost
It is worth investigating the length L traversed by the system in probability space and the cost C of the bangbang process. There is a trade-off between operation and cost, as expressed by the "thermodynamic uncertainty relation" (16). As a result, it is expected that minimising 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 In the following we calculate the lowest order contribution to the length and cost, i.e. and First, we look into the length swept by the probability distribution. Taking into account that T J > T f > 1, we have to split the integral in Eq. (48) into two summands: the first one for the time interval [0, t J ], in which the temperature monotonically increases from T i = 1 to T J , and the second one for [t J , t f ], in which the temperature monotonically decreases from T J to T f . Therefore, we have that 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 L rel G (blue broken) and the length over the geodesic Λ G (solid red), given in Eq. (19). Despite the heating-cooling bang-bang minimises the connection time, which is much shorter than the relaxation time t G R , the corresponding length is always larger than L rel G , since with the equality holding when T J = T f . It is inelasticity that makes T J different from T f and thus increases the length swept by the system in probability space. Only in the elastic limit α → 1 we have that L (0) → L rel G , as neatly observed in the plot.
Second, we consider the cost of this heating-cooling bang-bang process. Since the bang-bang process minimises 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 integral into the subintervals [0, t J ] and [t J , t f ] and integrating over temperature in each subinterval, one gets We have taken into account that . Also, and consistently, we have put a s 2 = 0 in the evolution equation for the temperature of Eq. (8), since we are evaluating C (0) to the lowest order.
Our main conclusion is thus that C (0) diverges for the optimal bang-bang connection, If χ max 1 but not infinite, the above equation gives the leading behaviour of the cost. In that case, t J = O(χ −1 max ) is small and the cooling part still rules the operation time, t f − t J t J , whereas the heating pulse still prevails for the cost.
Despite the divergence of C (0) , the energy input from the stochastic thermostat remains finite. Since the energy of the granular fluid is proportional to the temperature, we identify energy with temperature in the following discussion. The stochastic forcing is switched off in the cooling step of the bang-bang process, so the energy input comes from the heating step and equals In fact, for a fixed connection time t f , the bang-bang process minimises the energy input by the stochastic thermostat-and thus also the energy dissipated in collisions, since the total energy increment is given by T f − 1 for the initial and final NESS. The reason is simple: the only change is Pontryagin's scheme is that of the function f 0 , in which f 0 = 1 for time minimisation-see (27)-is substituted with f 0 = χ 1 + 3 16 a s 2 for energy input minimisation. Still, Pontryagin's Hamiltonian Π is linear in the control χ and therefore the bang-bang protocol also emerges as the optimal solution in this case.
B. Cooling-heating bang-bang: T f < 1 Now we turn our attention to the case in which the target temperature is smaller than the initial one, T f < 1. Similarly to what we have done in the previous section, we consider the two-step bang-bang process but with the order of the bangs reversed: χ = 0 for 0 ≤ t ≤ t J , and χ(t) = χ max for t J ≤ t ≤ t f , as described in Sec. V B.

Minimum connecting time
In the second part of the process, a line of reasoning similar to the one leading to Eq. (41) gives us that This means that t f → t J , the second part of the process is instantaneous in the limit as χ max → ∞. On the other hand, the system freely cools in the first part of the process, and then where A 2 (T ) in now given by Eq. (37). Again, the integral cannot be carried out analytically but it is possible to derive an approximate expression for t f by recalling that a s 2 is small and A 2 = O(1). In this way, we obtain where Once more, the last two equations give the connecting time t We show the behaviour of t (0) f as a function of the target temperature in Fig. 7, for T f < 1. All curves correspond to d = 2 but different values of the restitution coefficient α. In this case, t (0) f beats the speed limit t (1) R for relaxation for high enough T f -in the limit as R remains finite-but lies above it for T f 0.15, approximately. This contrasts with the situation for T f > 1, shown in Fig. 4. Physically, this asymmetry between the cases T f > 1 and T f < 1 can be understood as stemming from the non-holonomic constraint χ ≥ 0, which limits the rate at which the system can be cooled down-whereas no such limit exists for T f > 1 because we have considered that χ max → ∞.
At difference with the case T f > 1, t (0) f depends very weakly on α [64]. Since the excess kurtosis is small, we can obtain a rough estimate of the behaviour of the system by completely neglecting it-the so-called Gaussian approximation, which we have already employed in Sec. IV. Therein, it is clear that only one bang with χ = 0 suffices: the fastest way of reaching a temperature T f below the initial one is to turn off the stochastic thermostat. Putting χ and a s 2 to zero in Eq. (8), we obtain the Gaussian estimate for the connecting time t G f = 2 T −1/2 f − 1 , which is also plotted in Fig. 7.
We can obtain asymptotic expressions for t , which leads to In the elastic limit, A HCS 2 → 1 + , and thus t but it is of the same order of magnitude-for example, t f /t G f ∼ (A HCS 2 ) 1/4 for T f 1. This is in accordance with the behaviour observed in Fig. 7. In the linear response limit, 1 − T f 1, the behaviour is completely similar to that of T f > 1: Eq. (46) still holds replacing T f − 1 with its absolute value.

Associated length and cost
Let us evaluate the length and cost of the coolingheating bang-bang protocol. We start by calculating the length L (0) : here, T J < T f < 1 and once more the integral Eq. (48) has to be split into the two time subintervals inside which the temperature is monotonic, which yields Figure 8 is devoted to the comparison of the length L (0) over the optimal bang-bang connection with the lengths for the relaxation process L rel G and over the geodesic Λ G , which are given by Eq. (19). Similarly to the case T f > 1, L (0) is always larger than its value for relaxation L rel Again, the equality only holds when T J = T f , which corresponds to the elastic limit. Now we move on to the cost of the cooling-heating bang-bang. One more time, we split the time integral and change to integrate over temperature in each part of the bang-bang process. For the case T f < 1 we are analysing, the order of the bangs is reversed: Then we obtain that Taking into account that χ max 1, the second term on the rhs is negligible against the first one and thus In complete analogy to the T f > 1 case, Eq. (63) continues to give the leading behaviour of the cost when χ max 1 but not infinite, and the connection time and cost are still dominated by the cooling and heating bangs, respectively. Since the stochastic thermostat only heats the system in the time interval [t J , t f ], the energy input is now given by which is also finite. For the same reasons as in the case T f > 1, here the bang-bang protocol also minimises the energy input from the stochastic thermostat for a fixed connecting time t f .

VII. NUMERICAL SIMULATIONS
In order to check our theoretical predictions, we have carried out numerical simulations of the dynamics of the granular gas. Specifically, we have carried out Direct Simulation Monte Carlo (DSMC) for the two-dimensional case and two different values of the restitution coefficient, α = 0.3-for which a s 2 is positive-and α = 0.8-for which a s 2 is negative. In all cases, we start from a high temperature state with a Maxwellian velocity distribution function and switch on the stochastic thermostat with a certain intensity ξ i : the granular gas relaxes towards the corresponding NESS, in which the temperature T i and the noise intensity are related by Eq. (2a). Recall that we have employed T i to non-dimensionalise the temperature, so in our units T i = 1. From this initial NESS, we implement the bang-bang protocols developed in the previous sections.
In the case T f > 1, we proceed as follows in each trajectory of the simulation. First, the system is instantaneously heated from T i = 1 to T J : we make the velocities of all particles change as v i → v i + η i , where η i are independent Gaussian distributed random variables of a certain variance-the larger the variance, the larger the temperature increment T J − 1 and the smaller the excess kurtosis a 2J . Second, starting from the previously generated configuration, we let the system freely cool (ξ = 0) until a 2 in the trajectory equals the steady value a s 2 . This determines the connecting time t f , at which the temperature in the trajectory equals T f . At this time, we switch on the stochastic thermostat again but with an intensity ξ f such that the system remains stationary for t > t f : taking advantage of the theoretical prediction ξ ∝ T 3/4 s , as given by Eq. (2a), we set ξ f = ξ i T 3/4 f . The quantities T J , a 2J , t f , T f , and ξ f fluctuate from one realisation to another. A typical trajectory of the case T f > 1 is depicted in Fig. 9. Specifically, the realisation shown corresponds to d = 2 and α = 0.8 in a system with N = 10 6 particles. It is neatly observed how the system remains stationary after the stochastic forcing is switched on at t f . Note that fluctuations in the excess kurtosis are much larger than those of the temperature-which are basically not seen in the scale of the figure. Figure 10 shows the connecting time t f as a function of the target temperature T f . Once more, we consider the two-dimensional case and two different values of the restitution coefficient, α = 0.3 and α = 0.8. The simulation results are averaged over 100 trajectories and compared with the theoretical prediction (43), showing a very good agreement. The simulation curve is smoother for α = 0.3 than for α = 0.8, because |a s 2 | is larger for the former. In the case T f < 1, the cooling-heating bang-bang trajectory is generated in the following way. First, the system freely cools from the initial configuration, with T i = 1, until reaches a certain configuration with T J < 1 and a larger-in absolute value-excess kurtosis a 2J . Therefrom, we instantaneously heat the system by changing the velocities as v i → v i + M j=1 η ij , where η ij are independent Gaussian distributed random variables of a certain-small-variance. Note that, in contrast to the heating-cooling case described before, this is not done in one step but several. This recurrent procedure stops when the excess kurtosis-the absolute value of which is decreasing-equals its steady value a s 2 : this fixes the number of steps M over the considered trajectory. At this point, the temperature of the system is T f and, again, the stochastic forcing is switched on with intensity ξ f = ξ i T 3/4 f -this makes the system stationary for longer times. A typical trajectory of the case T f < 1 is depicted in Fig. 11. Specifically, the realisation shown corresponds to d = 2 and α = 0.3 in a system with N = 10 6 particles.
We compare the numerical results for the connecting time with the theoretical prediction, as given by Eq. (43), in Figure 12. Again, simulations correspond to d = 2, and α = (0.3, 0.8). The agreement between theory and simulations is excellent. Relative fluctuations seem to be smaller than in Fig. 10, but it has to be taken into account that the connection times here are longer.  This plot is similar to that in Fig. 9, but the order of the bangs is reversed: first, the granular gas freely cools (ξ = 0) in the time interval (0, t f ) and second, at t = t f , the system is instantaneously heated. Again, the thermostat is switched on with intensity χ f = T 3/2 f at t = t f and the system remains stationary for t > t f .

VIII. PREVALENCE OF OPTIMAL BANG-BANG PROTOCOLS
In this section, we discuss the emergence of similar bang-bang protocols as the optimal ones-in the sense of minimising the connection time-in a wide class of physical systems, not only for the specific case of granular fluids. The evolution equations of the relevant physical quantities are often linear in the "control functions" [11,16,[45][46][47]. In addition, the controls (stiff- ness, temperature, diffusions coefficient, noise strength) are non-negative in most of these situations, which gives rise to non-holonomic constraints.
Non-holonomic restrictions make it necessary to resort to Pontryagin's maximum principle to solve the problem of minimising the connection time. When the evolution equations involve the controls λ linearly, Pontryagin's Hamiltonian Π is also linear in them. Therefore, the the maximum of Π is reached at the limit values of the control, λ min and λ max : the bang-bang protocols arise in this way. Once the bang-bang protocols are established as the optimal ones, a relevant question is the number of "bangs" that are needed for the connection. The answer to this question in each specific physical system depends on the number n of variables that characterise its state. If the evolution equations are also linear in these physical variables-either exactly or as a consequence of linearising them around a steady state, a rigorous mathematical theorem ensures that the optimal bang-bang control has at most n − 1 switchings [65].
If the evolution equations are non-linear in the physical quantities, as is the case of granular fluids, the following argument supports the arising of a completely similar picture. A simple bang-with two possibilities, λ(t) = λ min or λ(t) = λ max -suffices for n = 1. Its duration t f makes it possible to adjust the final value of the only one variable. Two bangs-with two possibilities, λ max − λ min or λ min − λ max -suffice for n = 2, with only one switching at an intermediate time t J . The joining time t J together with the connecting time t f -or the duration of the bangs t 1 = t J and t 2 = t f −t J -make it possible to tune the final value of the two variables. This is the case we have found here when studying the granular fluid in the Sonine approximation. In general, n bangs-with two possibilities, starting with either λ max or λ min -suffice for a generic value of n, with n − 1 switchings at intermediate times t J1 , t J2 , . . . , t Jn−1 . These n − 1 joining times together with the connecting time t f -or the duration of the n bangs t 1 = t J1 , t 2 = t J2 − t J1 , . . . , t n = t f − t Jn−1 -make it possible to tune the final value of the n variables.
A. Brownian particle moving in a d-dimensional harmonic potential As a proof of concept, let us consider a Brownian particle trapped in a d-dimensional harmonic potential. Our analysis will be carried out in the overdamped limit, in which the probability distribution P (x, t) of the particle's position obeys the Fokker-Planck equation where γ is the friction coefficient and U h (x) = 1 2 d j,k=1 λ jk x j x k . We want to connect two equilibrium states in the minimum possible time, by controlling the temperature T (t) of the bath inside which the colloidal particle is immersed.
It is fitting to describe the particle position in terms of the three normal coordinates ξ i , such that the harmonic potential is diagonal in them. Therefore, we write where κ β > 0 are the eigenvalues of the matrix of elements k ij . The time evolution of the particle is uniquely characterised by the value of the variances of the normal modes, σ 2 β ≡ ξ 2 β . In dimensionless variables, they evolve according to Similarly to what we did in the analysis of the granular fluid, the units that non-dimensionalise variables are chosen to simplify our formulas: here, T i = 1 and κ 1 = 1. Without loss of generality, we choose to label the modes such that κ 1 = 1 ≤ . . . ≤ κ d . See Appendix D for details. In the initial and final equilibrium states, with respective initial temperatures T i = 1 and T f , we have The connection between the initial and final states is done by controlling the temperature of the bath T (t), which obeys the non-holonomic constraint T ≥ 0. Since the evolution equations (66) are linear in the control T (t), the optimal connection is of bang-bang type: it comprises several time windows with either T (t) = T max or T (t) = T min = 0. An exhaustive analysis of the optimal connection for this system, investigating in detail the behaviour of the connection time throughout the whole space of parameters (κ, T max , T f ), where κ = (κ 1 , . . . , κ d ), is outside the scope of this paper. Here, we focus on the general trends of the optimal connecting time as a function of the final temperature T f for some specific choices of κ, in order to show the generality of the arguments we have put forward above. For the sake of simplicity, we consider the limiting case T max → ∞-in analogy with our study of the granular fluid, in which the maximum noise strength χ max → ∞.
One-dimensional case.-The simplest situation is that of the one-dimensional case, where we have only one equation to control. Thus, one bang suffices of duration t (I) f . It is interesting to remark that the same situation appears in d = 2 or d = 3 for the isotropic or central symmetry situation, in which all the κ β are identical: both (66) and (67) do not depend on the mode and all σ β (t) are also identical. From an effective point of view, we have n = 1 and the control problem is equivalent to that of the one-dimensional case.
On the one hand, in the "heating" bang, (I) f = T f . In the limit as T max → ∞, this expression simplifies to Therefore, we have that Thus, any T f > 1 can be reached by an instantaneous, infinitely strong, jump to T max → ∞ that is stopped as soon as σ 2 1 attains its final value. On the other hand, in the "cooling" bang, T (t) = 0, σ 2 1 = e −2t and the connection time is determined by As a consequence of the non-holonomic constraint T (t) ≥ 0, we have that the optimal time is finite for T f < 1. Two-dimensional case.-Next we look into a the twodimensional case, with κ 2 = κ 1 . We have two equations to control and thus n = 2. It is clear that only one bang is not enough: in a "cooling" bang of duration t Therefore, we need two bangs: either "heating-cooling" or "cooling-heating", consistently with our general discussion for n = 2. Note that for d = 3 but with axial symmetry, for example σ 2 (t) = σ 3 (t), the control problem is equivalent to the two-dimensional situation considered here.
First, we analyse the "heating-cooling" bang-bang. In the first step, we have an instantaneous heating with T (t) = T max → ∞ that leads to The duration t 1 = t J of the heating bang vanishes in the limit as T max → ∞. In the second step, we put T (t) = 0 this means that T f > e −2t (III) if, from going from I to III, one incorporates a new value of κ while keeping the previous ones-i.e. if κ 1 and κ 2 in III are the same as in II. This is logical, the incorporation of additional variables reduces the size of the set of "admissible" controls-i.e. those connecting the initial and final states-and thus increases the minimum connection time. Figure 13 shows the connection time as a function of the final temperature, for three sets of the model parameters (κ 1 = 1, κ 2 , κ 3 ). In accordance with the physical picture of the previous paragraph, it is neatly observed that t = 0 for T f > 1 whereas it increases as the final temperature decreases for T f < 1, diverging in the limit as T f → 0. In case II, the observed behaviour is also qualitatively similar with that of the Sonine description of the granular system, the main difference being that t (II) f is no longer zero for T f > 1. The existence of two variables makes it impossible to connect the two states instantaneously, because the cooling part also involves a finite time due to the non-holonomic constraint T ≥ 0-analogous to χ ≥ 0 in the granular case. The main difference is the finite value of t (II) f for the twodimensional oscillator-or a three-dimensional oscillator with axial symmetry-in the limit as T f → ∞, in contrast with the vanishing of the connecting time for the granular case. The latter behaviour 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 minimising the connection time.
It is known that inverse engineering techniques-such as engineered swift equilibration [11][12][13][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 given by (69) and (70) for the one-dimensional oscillator, κ1 = 1 (or d-dimensional one with central simmetry, κ β = 1 ∀β); (b) t (II) f given by (73) and (75) for the two-dimensional oscillator (or three-dimensional with axial symmetry), with κ = (1, 2); and (c) t (III) f given by (78) and (80) for the three-dimensional case, with κ = (1, 2, 3).
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 T f > 1 and T f < 1-recall that in our dimensionless units the initial temperature equals unity. For the case T f > 1, the optimal connecting times are rather small, vanishing in the limits T f → 1 and T f → ∞. 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 T f . Therefore, it is non-Gaussianities-specifically, the excess kurtosis a 2 that is small-that make impossible to instantaneously connect the two NESS for T f > 1. The excess kurtosis decreases and therefore the state after the instantaneous heating pulse is not stationary.
For the case T f < 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 T f > 1, the minimum connection time for T f < 1 does not vanish because free cooling involves a finite time. Interestingly, both for T f > 1 and T f < 1 non-Gaussianities make the connecting times longer: this is physically understood by taking into account that non-Gaussianities stem from the inelasticity of collisions.
One of the main results of our paper is the emergence of optimal bang-bang protocols for minimising 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 T f > 1 and T f < 1: heating-cooling for T f > 1, but cooling-heating for T f < 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 minimising 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 non-holonomic constraint limits their physically acceptable values, λ min ≤ λ ≤ λ max . Examples abound, from the trapped Brownian particle [11,45] or active systems [46,47] to a particle moving in an electric field [16]). Therein, the mathematical structure of Pontryagin's principle ensures that the optimal controls minimising 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 n physical quantities [67]. The number n is thus small, in the examples above n = 1 [11,16,45] or n = 2 [46,47]-as is the case of the granular fluid in the Sonine approximation. A relevant question is: 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 t f makes 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 viceversa-of durations t 1 and t 2 , with t f = t 1 + t 2 , 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/decompression of a Brownian particle trapped in a d-dimensional harmonic potential, by controlling the temperature of the bath [68]. For n = 2 (axial symmetry for d = 3 or d = 2), a two-step bang-bang process, qualitatively similar to that found in the granular fluid, carries out the optimal connectionheating-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 bangbang process arises: heating-cooling-heating or coolingheating-cooling, again depending on the final value of the temperature.
Our work has been focused on the minimisation of the connecting time t f between two NESS of the granular fluid. Not only is the optimisation 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, minimising 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 maximises the delivered power [69]. Also, the optimisation of the relaxation route to equilibrium or to a NESS is of interest in connection with behaviours such as the Mpemba effect, which is currently a very active field of research [53,[70][71][72][73][74].
In addition, we have considered the associated statistical length L and cost C over 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 t f ≥ L 2 /(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 minimised 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 timewhen minimising the cost-is the counterpart of the divergence of the cost-when minimising 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 ddimensional oscillator, show that this will be also the case in many physical situations where a non-holonomic 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 [30,31,75,76]. In this sense, the role of the Fisher information for rigorously establishing the H-theorem for granular gases is worth investigating.
For Markovian dynamics, the cost C has been related in general to entropy production [3,5,32,33]. 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. [29]. Therefore, elucidating the physical meaning of information geometry's costbeyond 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 densityor 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 non-equilibrium systems such as active matter [77][78][79][80][81][82]. 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 [83][84][85]. Therefore, it is worth looking into the application of information-geometry concepts and, in general, the extension of our results to these physical contexts. nomial with five coefficients: four to adjust the boundary conditions (23) for the temperature, and one extra parameter to impose that A 2p (t f ) = 1.
To start with, it is adequate to employ the scaled time τ = t/t f introduced in Appendix A 2 and to work with the thermal velocity v th ≡ √ T . Consistently, v th,p (τ ) = √ T p (τ ), and we rewrite Eq. (24)  (B1) Insertion of this expression for the noise intensity into the evolution equation of the excess kurtosis gives, after some algebra, We solve this equation-with the initial condition A 2p (0) = 1-with the following 4-th order polynomial for the thermal velocity, v th,p (τ ) = 1+cτ 2 +(4∆v th −2c)τ 3 +(c−3∆v th )τ 4 . (B3) The parameter c is tuned to meet the boundary condition A 2p (t f ) = 1: there is only one 4-th 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 T f > 1 in Fig. 14, specifically for T f = 4 (∆v th = 1). The following qualitative behaviour is observed: as the connecting time t f in 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 4-th order polynomial for too short times [86].
The observed behaviour hints at the emergence of a minimum, non-vanishing, value of the connecting time for the ESR protocol, both for T f > 1 and T f < 1. This feeling is reinforced by employing higher order polynomials. For example, in the 5-th order case, there is a mono-parametric family of polynomials connecting the initial and final NESS. Nevertheless, χ p (t) becomes negative for t f below a certain value, over the whole family of polynomials making the connection.

Appendix C: Fisher information in the Sonine approximation
Here, we look into the Fisher information I(t) within the Sonine approximation we have considered through- out. The velocity distribution function is expanded as where P G (v, v th (t)) is the Maxwellian distribution of Eq. (A1), a k (t) are the coefficients of the expansionwhich are related to the cumulants, and finally S k (x) ≡ L ( d−2 2 ) k (x) are the associated Laguerre (or Sonine) polynomials [87]. The explicit expression for the first Sonine polynomials are [48] S 0 (x) = 1, The first Sonine approximation consists in keeping only the first term in the expansion, with coefficient a 2 that equals the excess kurtosis, and neglecting nonlinear contributions in a 2 in all the expressions derived from Eq. (C1) [88]. For our purposes, it is convenient to rewrite Eq. (C1) as follows: we introduce a dimensionless velocity c(v, T (t)) = v/ 2T (t), and the order of unity quantity A 2 (t) defined in Eq. (7), so that P (v, t) = e −c 2 [2πT (t)] d/2 1 + a s 2 A 2 (t)S 2 (c 2 ) .
Now we proceed to calculate the Fisher information. For that, we take into account that In order to obtain I(t), Eq. (C6) is squared and averaged with the probability distribution (C4)-neglecting nonlinear terms in a s 2 . After a little algebra, one gets where we have omitted the time dependence of T (t) and A 2 (t) to simplify the notation, and defined There is no contribution coming from the term proportional toȦ 2 in Eq. (C6) because of the orthogonality of Sonine polynomials, S j (c 2 )S k (c 2 ) = 0 for j = k. Also, note that I (0) (t) = I G (t) because we no longer set the excess kurtosis to zero in the first Sonine approximation. Notwithstanding, the smallness of a s 2 implies that the main contribution to the Fisher information comes from I (0) (t).