Optimal Control of Uniformly Heated Granular Fluids in Linear Response

Abstract

We present a detailed analytical investigation of the optimal control of uniformly heated granular gases in the linear regime. The intensity of the stochastic driving is therefore assumed to be bounded between two values that are close, which limits the possible values of the granular temperature to a correspondingly small interval. Specifically, we are interested in minimising the connection time between the non-equilibrium steady states (NESSs) for two different values of the granular temperature by controlling the time dependence of the driving intensity. The closeness of the initial and target NESSs make it possible to linearise the evolution equations and rigorously—from a mathematical point of view—prove that the optimal controls are of bang-bang type, with only one switching in the first Sonine approximation. We also look into the dependence of the optimal connection time on the bounds of the driving intensity. Moreover, the limits of validity of the linear regime are investigated.

1. Introduction

The study of granular media, beyond its own theoretical interest, is particularly important for industrial applications such as improving their transport or storage. Granular materials are discrete clusters of macroscopic particles that exhibit two fundamental features. First, collisions between particles are inelastic, so that energy is not conserved: It monotonically decreases with time if there is no external mechanism that injects energy into the system. Second, thermal energy is many orders of magnitude lower than the characteristic potential energy, making thermal fluctuations largely irrelevant for the behaviour of granular systems [1].

In the simplest model for granular fluids, particles are d-dimensional smooth, hard spheres of mass m that undergo inelastic binary collisions. In each collision, the tangential component of the relative velocity is unchanged, whereas the normal component is reversed and shrunk by a factor $\alpha$, $0\le \alpha \le 1$, which is termed the restitution coefficient. Energy is only kinetic, and the energy dissipated in each collision is thus proportional to $1-{\alpha }^2$—the elastic limit corresponds to $\alpha =1$. In the undriven system, after a few collisions per particle, the so-called homogeneous cooling state (HCS) is reached [2,3,4,5,6,7], in which the system remains homogeneous and the granular temperature T—basically the average kinetic temperature—monotonically decreases following an algebraic decay, the Haff law [2].

In order to allow the system to reach a stationary state, an energy injection mechanism is needed. A simple but also relevant situation is the uniformly heated granular fluid [8,9] that we consider throughout this work. Therein, independent white noise forces act on the particles of the granular fluid, the intensity of which is characterised by a parameter $\chi \ge 0$ related to the variance of the stochastic force. The granular fluid reaches a non-equilibrium steady state (NESS) in the long-time limit, in which the system remains homogeneous. Therein, the energy injected by the stochastic thermostat balances—in average—the energy loss in collisions and the value of the granular temperature depends on the intensity of the driving, whereas higher-order cumulants of the velocity are independent thereof. The steady state is always stable, at variance with the HCS, which is known to be unstable for perturbations of large enough wavelength [4,10,11] (On another note, velocity correlations stemming from the dissipative character of collisions are relevant in certain physical situations, e.g., for the understanding of the total energy fluctuations in the HCS [12,13]). The uniformly heated granular gas has been extensively studied, both its properties at the NESS [8,9,14,15] and its dynamical evolution [16,17,18].

Granular systems are intrinsically out-of-equilibrium systems. Their dissipative dynamics entail that their velocity distribution function (VDF) is non-Gaussian, even in the long-time limit in which a hydrodynamic state, independent of the initial condition, is reached. This is true for both the HCS in the undriven case and the NESS in the uniformly heated situation. The non-Gaussianities of the VDF are essential to understand the behaviour of granular fluids and are incorporated to the picture by implementing a Sonine expansion [3] of the Enskog–Fokker–Planck equation. This leads to an infinite hierarchy of equations for the cumulants, which is typically closed by introducing the so-called first Sonine approximation: Only the fourth cumulant or excess kurtosis $a_2$ is retained—higher order cumulants are neglected. Therein, the granular temperature and the excess kurtosis obey a system of two coupled ordinary differential equations, the accuracy of which for describing the dynamical evolution of the granular fluid has been validated in many works, e.g., [8,9,15,17,18,19,20,21]. In this context, especially relevant are those analysing memory effects such as the Kovacs hump or the Mpemba crossing, in which non-Gaussianities are key to facilitate their emergence [18,19,20,21].

Only very recently has the possibility of controlling the dynamical evolution of granular systems been analysed [22]. This might be surprising at first sight, since the control of physical systems has been considered for some time in different physical contexts, such as quantum mechanics [23,24,25,26] and statistical mechanics [27,28,29,30,31,32]. A paradigmatic case of control of a mesoscopic system is that of an optically trapped colloidal particle [27,28,31,33,34,35,36,37,38,39]. When the confining potential is harmonic, the time dependence of the stiffness of the trap $\kappa \left(t\right)$ can be externally controlled, and one aims at optimising the connection between two given equilibrium states, corresponding to different values of the stiffness of the trap—i.e., the colloidal particle is being confined or deconfined. Here, optimising means that some relevant physical observable (irreversible work, entropy production, connection time, …) is minimised. The time-dependent stiffness $\kappa \left(t\right)$ plays the role of the control function—sometimes together with the temperature of the bath, which can be changed in an effective way by adding a random force [40,41]. The control problem is greatly simplified by the following three features. First, the initial and target states are equilibrium states, so that their corresponding probability distribution functions (PDFs) are perfectly known. Second, the PDF is Gaussian for all times, so that it is completely characterised by its average and variance. Third, the evolution equations for the average and the variance are exactly solvable in closed form.

The delay in posing the problem of controlling granular systems probably stems from the challenging character of the control problem in this case, both at the conceptual and mathematical level. None of the three simplifying features above, holding for the harmonically trapped Brownian particle, is present in granular fluids. First, the initial and target states are NESS, and their PDFs are only approximately known. Second, the PDF is non-Gaussian for all times. Third, the evolution equations are non-linear and thus not exactly solvable. It is interesting to compare the situation in the granular case described above with the one appearing in other paradigmatic system, the Brownian gyrator [42,43,44]. Although the initial and final states are also NESSs in that case, the PDF is Gaussian for all times, and the evolution equations for the relevant moments can be exactly solved. Thus, the control problem of this system is simpler, although only non-optimal connections have been worked out, to the best of our knowledge [45].

One may thus pose the problem of connecting two NESSs of the granular fluid corresponding to different values of the driving intensity $\chi$, ${\chi }_{\mathrm{i}}$ and ${\chi }_{\mathrm{f}}$, i.e., to different values of the granular temperature $T_{\mathrm{i}}$ and $T_{\mathrm{f}}$. The control function here is the intensity of the driving $\chi \left(t\right)$. We are interested in the time optimisation problem, i.e., to find the protocol $\chi \left(t\right)$—starting from (and ending at) the desired initial (and target) NESS—that minimises the connection time between the initial and final states. This kind of time optimisation problem is important from a fundamental point of view and also has relevance for applications. For the connection between equilibrium states, related problems emerge in the optimisation of irreversible heat engines [46], the analysis of the Mpemba effect [21,47,48,49], and the optimisation of the relaxation route to equilibrium [50,51,52].

The limiting situation in which all the power of the stochastic thermostat is available, i.e., $0\le \chi <\infty$, was investigated in Ref. [22] within the first Sonine approximation. Despite the challenges mentioned above, the unboundedness of the control makes it possible to give analytical predictions for the connecting time—the evolution equations are heavily simplified in the limiting cases $\chi =0$ and $\chi =\infty$. In addition, the analytical predictions were compared with numerical simulations, and an excellent agreement was found. Indeed, this was expected because the accurateness of the first Sonine approximation is a well-established fact for the smooth hard-sphere granular gas, even for the study of subtle behaviours such as the Kovacs or Mpemba memory effects [18,19,20,21]. This is the reason why we follow a purely analytical approach in this paper.

In this work, we analytically investigate the more realistic case in which the driving intensity is bounded between two values, ${\chi }_{min}\le \chi \le {\chi }_{max}$. In order to make analytical progress, we consider the linear response regime, in which ${\chi }_{min}$ and ${\chi }_{max}$ are close: This allows us to linearise the evolution equations and make exact—in the linear response limit—predictions for the optimal connecting time as a function of the bounds $({\chi }_{min},{\chi }_{max})$. The linearisation of the equations also allows us to employ rigorous mathematical results of optimal control theory (OCT) and to check that the underlying hypotheses are fulfilled, a program that was unattainable in the non-linear case [22]. Moreover, we also explore the limits of validity of the linear response regime, by taking the double limit $({\chi }_{min}\ll 1,$ ${\chi }_{max}\gg 1)$ and comparing the obtained behaviour with those for the non-linear case with unbounded driving [22].

The structure of this paper is as follows. In Section 2, we put forward the model, write the evolution equations for the temperature and the excess kurtosis, and linearise them around the final NESS. Section 3 is devoted to the derivation of the optimal controls, in the sense of minimising the connection time. The trajectories of the temperature and the excess kurtosis—both as functions of time and in the phase plane—for the optimal controls are analysed in Section 4. The dependence of the minimum connection time on the bounds of the driving is the subject of study of Section 5. We investigate the limits of the validity of the linear response approximation as the bounds in the driving are loosened in Section 6. Finally, a discussion of the obtained results is presented in Section 7. The Appendices deal with some technicalities that are omitted in the main text.

2. The Model

Our system is a granular fluid with number density n, comprising N d-dimensional ($d=2,3$) hard-spheres of mass m and diameter $\sigma$ (hard discs in $d=2$). Specifically, we consider smooth inelastic hard spheres. Collisions between them are binary, and the post-collisional velocities $(v_1^{\prime },v_2^{\prime })$ are given in terms of the pre-collisional ones $(v_1,v_2)$ by:

$$v_1^{\prime }=v_1-\frac{1+\alpha }{2}(v_{12}·\widehat{\mathit{\sigma }})\widehat{\mathit{\sigma }},v_2^{\prime }=v_2+\frac{1+\alpha }{2}(v_{12}·\widehat{\mathit{\sigma }})\widehat{\mathit{\sigma }},$$

(1)

where $\widehat{\mathit{\sigma }}$ is the unit vector along the direction joining the centre of the particles and $\alpha$ is the restitution coefficient, $0\le \alpha \le 1$. In addition, the system is heated by a stochastic thermostat, i.e., a white-noise force ${\mathit{F}}_i$ independently acts on every particle verifying $〈{\mathit{F}}_i\left(t\right)〉=0$, $〈{\mathit{F}}_i\left(t\right){\mathit{F}}_j\left(t\right)〉=m^2{\xi }^2{\delta }_{ij}\delta (t-t^{\prime })$, $\forall i,j=1,\dots ,N$, and $\forall (t,t^{\prime })$.

In the first Sonine approximation that we employ throughout, the system is described by two variables, the granular temperature T and the excess kurtosis $a_2$. Their definitions in terms of moments of the velocity are:

$$〈v^2〉=\frac{dT}{m},a_2=\frac{d}{d+2}\frac{〈v^4〉}{{〈v^2〉}^2}-1.$$

(2)

As stated in the introduction, the system reaches an NESS in the long-time limit due to the balance—on average–of the energy input and dissipation. The stationary values of T and $a_2$ are given by:

$$T_{\mathrm{s}}^{3/2}=\frac{m{\xi }^2}{{\zeta }_0(1+\frac{3}{16}{a}_2^{\mathrm{s}})}\equiv \chi ,{\zeta }_0=\frac{2n{\sigma }^{d-1}(1-{\alpha }^2){\pi }^{\frac{d-1}{2}}}{\sqrt{m}d\mathsf{\Gamma }(d/2)},$$

(3)

$$a_2^{\mathrm{s}}=\frac{16(1-\alpha )(1-2{\alpha }^2)}{73+56d-24d\alpha -105\alpha +30(1-\alpha ){\alpha }^2},$$

(4)

Note that $a_2^{\mathrm{s}}$ is independent of the thermostat intensity, as measured by $\chi$, it only depends on $(d,\alpha )$. From the kinetic equation, the following coupled system of ordinary differential equations (ODEs) are obtained—see, e.g., [17,19]:

$$\dot{T}={\zeta }_0\left[\chi (1+\frac{3}{16}{a}_2^{\mathrm{s}})-T^{3/2}(1+\frac{3}{16}{a}_2)\right],$$

(5a)

$$\dot{a_2}=\frac{2{\zeta }_0}{T}\left[(T^{3/2}-\chi )a_2+B{T}^{3/2}(a_2^{\mathrm{s}}-a_2)\right],$$

(5b)

where the parameter B is given by [17,19]

$$B=\frac{a_2^{\mathrm{HCS}}}{a_2^{\mathrm{HCS}}-a_2^{\mathrm{s}}},a_2^{\mathrm{HCS}}=\frac{16(1-\alpha )(1-2{\alpha }^2)}{25+2{\alpha }^2(\alpha -1)+24d+\alpha (8d-57)};$$

(6)

$a_2^{\mathrm{HCS}}$ is the value of the excess kurtosis in the HCS [8,9].

3. Optimal Control in Linear Response

Above, we have considered that the driving intensity of the thermostat $\chi$ is constant. In general, it may be time-dependent, a certain given function of time $\chi \left(t\right)$ that determines the externally enforced driving program. Looking at the evolution Equation (5a,b) for $(T,a_2)$ in the light of OCT, this means that $\chi \left(t\right)$ is the control function. In this paper, we consider the following control problem: the connection of two NESS, i.e., bringing the system from an initial state $(T_{\mathrm{i}},a_{2\mathrm{i}}=a_2^{\mathrm{s}})$, to a target, final one $(T_{\mathrm{f}},a_{2\mathrm{f}}=a_2^{\mathrm{s}})$, by engineering a suitable driving program $\chi \left(t\right)$. Moreover, we would like to perform this connection in the shortest possible time. The case in which all the power of the thermostat is available has been considered in Ref. [22]. Here, we analyse the more realistic case in which the driving intensity is bounded between two limiting values, ${\chi }_{min}\le \chi \left(t\right)\le {\chi }_{max}$, with ${\chi }_{min}\ge 0$ and ${\chi }_{max}<\infty$.

In order to solve the control problem analytically, we restrict ourselves to the linear response regime, i.e., ${\chi }_{min}$ and ${\chi }_{max}$ are close—and so are $T_{\mathrm{i}}$ and $T_{\mathrm{f}}$. To look into the dynamics of the system, it is preferable to introduce scaled variables as follows:

$$t^{*}={\zeta }_0{T}_{\mathrm{f}}^{1/2}t,T^{*}=\frac{T}{T_{\mathrm{f}}},{\chi }^{*}=\frac{\chi }{T_{\mathrm{f}}^{3/2}},A_2=\frac{a_2}{a_2^{\mathrm{s}}}.$$

(7)

In this way, we have defined dimensionless time $t^{*}$, granular temperature $T^{*}$, and driving ${\chi }^{*}$; moreover, scaling the excess kurtosis with its steady value simplifies our analysis (Both $a_2^{\mathrm{HCS}}$ and $a_2^{\mathrm{s}}$ change sign for $\alpha =1/\sqrt{2}$, so that $a_2$ typically changes sign with the inelasticity. On the other hand, the scaled variable $A_2$ always remains positive). For these scaled variables, we have the following evolution equations:

$$\dot{T}=\chi (1+\frac{3}{16}{a}_2^{\mathrm{s}})-T^{3/2}(1+\frac{3}{16}{a}_2^{\mathrm{s}}{A}_2),$$

(8a)

$$\dot{A_2}=\frac{2}{T}\left[(T^{3/2}-\chi )A_2+B{T}^{3/2}(1-A_2)\right].$$

(8b)

We have omitted the superscript * in the dimensionless variables in order to simplify the notation since; from now on, these are the variables used. The term $\chi (1+\frac{3}{16}{a}_2^{\mathrm{s}})$ on the right hand side (rhs) of (8a) represents the energy injection due to the action of the thermostat, while the term $-T^{3/2}(1+\frac{3}{16}{a}_2^{\mathrm{s}}{A}_2)$ collects the energy losses due to the inelastic collisions. Of course, if $\chi$ is kept constant and equal to its target value, i.e., $\chi \left(t\right)=1$ $\forall t\ge 0$, the system reaches the NESS $(T_{\mathrm{s}}=1,A_2^{\mathrm{s}}=1)$ in the long-time limit, consistently with our discussion in the previous section.

In linear response, we thus write:

$$T=1+\delta T,A_2=1+\delta A_2,\chi =1+\delta \chi ,$$

(9)

with $\delta T\ll 1$, $\delta A_2\ll 1$, and $\delta \chi \ll 1$. Note that, to be consistent, we must assume that $\delta {\chi }_{min},\delta {\chi }_{max}\ll 1$. This allows us to linearise the evolution equations of T and $A_2$ as follows:

$$\frac{d}{dt}\left(\begin{array}{c}\delta T\\ \delta A_2\end{array}\right)=\left(\begin{array}{c}\beta \\ -2\end{array}\right)\delta \chi +\left(\begin{array}{cc}-\frac{3}{2}\beta & 1-\beta \\ 3& -2B\end{array}\right)\left(\begin{array}{c}\delta T\\ \delta A_2\end{array}\right),$$

(10)

where we have defined

$$\beta \equiv 1+\frac{3}{16}{a}_2^{\mathrm{s}}.$$

(11)

Now it is $\delta \chi$ that plays the role of the control function, $\delta {\chi }_{min}\le \delta \chi \le \delta {\chi }_{max}$.

Once the evolution equations are linearised, the control problem is stated as follows: We would like to bring the system from the initial NESS corresponding to

$$\delta T(t=0)=\delta T_{\mathrm{i}},\delta A_2(t=0)=0,$$

(12)

to the target NESS

$$\delta T\left(t_{\mathrm{f}}\right)=0,\delta A_2\left(t_{\mathrm{f}}\right)=0,$$

(13)

in the minimum possible time $t_{\mathrm{f}}$. Moreover, the system remains stationary for $t<0$ and $t>t_{\mathrm{f}}$: This means that for $t<0$, we have prepared the system in the NESS with the initial value of the temperature by driving it with the corresponding intensity and that for $t\ge t_{\mathrm{f}}$, the driving intensity for the target temperature is applied, i.e.,

$$\delta \chi \left(t\right)=\delta {\chi }_{\mathrm{i}}=\frac{3}{2}\delta T_{\mathrm{i}},t<0,\delta \chi \left(t\right)=0,t\ge t_{\mathrm{f}}.$$

(14)

Equation (10) is linear in both the variables $(\delta T,\delta A_2)$ and the control function $\delta \chi$, and therefore the rigorous theorems for linear control systems are applicable—see, for example, chapter III of Ref. [53]. For our specific situation in which $\delta {\chi }_{min}\le \delta \chi \le \delta {\chi }_{max}$, these theorems ensure that the optimal protocol that minimises the connection time $t_{\mathrm{f}}$ is of bang-bang type with at most one change. That is, $\delta \chi \left(t\right)$ is piece-wise continuous, taking either the value $\delta {\chi }_{max}$ or $\delta {\chi }_{min}$ and presenting, at most, one jump between these two values in the time window $(0,t_{\mathrm{f}})$ (More specifically, this result stems from Theorem 10 in Section 17 of Ref. [53], and we check that the hypotheses of this theorem are fulfilled in Appendix A). These kind of bang-bang optimal protocols arise in different physical situations [22,23,54,55,56]. In general, bang-bang protocols emerge as the optimal ones when Pontryagin’s Hamiltonian is linear in the controls—i.e., when the evolution equations are linear in the controls, although they may be non-linear in the relevant physical variables [30,57,58,59].

To be able to determine the optimal protocol, we must distinguish two cases according to the initial temperature, a global cooling process ($T_{\mathrm{i}}>T_{\mathrm{f}}=1$, $\delta T_{\mathrm{i}}>0$) and a global heating process ($T_{\mathrm{i}}<T_{\mathrm{f}}=1$, $\delta T_{\mathrm{i}}<0$):

• For $\delta T_{\mathrm{i}}>0$, CH protocol: In the time window $[0,t_J)$, the driving $\delta \chi$ is set to its minimum value $\delta {\chi }_{min}$ (cooling), whereas in the time window $[t_J,t_{\mathrm{f}})$, it is set to its maximum $\delta {\chi }_{max}$ (heating):

  $$\delta \chi \left(t\right)=\left\{\begin{array}{cc}\delta {\chi }_{\mathrm{i}},\hfill & t<0,\hfill \\ \delta {\chi }_{min},\hfill & 0\le t<t_J,\hfill \\ \delta {\chi }_{max},\hfill & t_J\le t<t_{\mathrm{f}},\hfill \\ 0,\hfill & t\ge t_{\mathrm{f}}.\hfill \end{array}\right.$$

  (15)
• For $\delta T_{\mathrm{i}}<0$, HC protocol: In the time window $[0,t_J)$, the driving is set to its maximum value ${\chi }_{max}$, whereas in the time window $[t_J,t_{\mathrm{f}})$, it is set to its minimum ${\chi }_{min}$:

  $$\chi \left(t\right)=\left\{\begin{array}{cc}\delta {\chi }_{\mathrm{i}},\hfill & t<0,\hfill \\ \delta {\chi }_{max},\hfill & 0\le t<t_J,\hfill \\ \delta {\chi }_{min},\hfill & t_J\le t<t_{\mathrm{f}},\hfill \\ 0,\hfill & t\ge t_{\mathrm{f}}.\hfill \end{array}\right.$$

  (16)

The switching time $t_J$, $0\le t_J\le t_{\mathrm{f}}$, will be determined later as a function of the parameters of the problem, i.e., as a function of $(\delta {\chi }_{\mathrm{i}},\delta {\chi }_{min},\delta {\chi }_{max})$. As already stated above, the values of $\delta \chi$ for $t<0$ and for $t\ge t_{\mathrm{f}}$ ensure that the system starts from the NESS with $T=T_{\mathrm{i}}$ and, after the application of the bang-bang protocol, remains in the target NESS with $T=T_{\mathrm{f}}=1$, for both the CH and HC protocols.

At first, there is no clear reason to assign the CH protocol to the case $\delta T_{\mathrm{i}}>0$ and the HC protocol to the case $\delta T_{\mathrm{i}}<0$ (aside from the analogy with the full-thermostat-power case analysed in Ref. [22]). In order to show that this is indeed the case, one needs to study the behaviour of the trajectories swept by the point in the phase plane $(\delta A_2,\delta T)$. We defer this analysis until Section 4.

3.1. $T_i>T_f$ = 1 Cooling-Heating Bang-Bang

In this section, we integrate the solution of the system in two time windows: the first one, $[0,t_J)$, when $\delta {\chi }_{min}$ is applied, and a second one $[t_J,t_{\mathrm{f}})$, when $\delta {\chi }_{max}$ is applied. We also determine the time $t_J$ as well as the value of the variables $\delta T_J$ and $\delta A_{2J}$ at that time. The point $(\delta A_{2J},\delta T_J)$ constitutes the set of initial conditions for the control system in the second window.

Equation (10) is inhomogeneous, due to the term proportional to $\delta \chi$ on its rhs. (It is only homogeneous when $\delta \chi =0$, i.e., when the control is set to the constant value ${\chi }_{\mathrm{f}}=T_f^{3/2}$ corresponding to the final temperature.) Over each time window, $\delta \chi \left(t\right)=\delta {\chi }_{\mathrm{ext}}$, where the subscript “ext” includes both bangs, $\delta {\chi }_{\mathrm{ext}}=\delta {\chi }_{min}$ (first window) and $\delta {\chi }_{\mathrm{ext}}=\delta {\chi }_{max}$ (second window). The inhomogeneity can be thus understood as the system relaxing towards the NESS corresponding to $\delta {\chi }_{\mathrm{ext}}$. Let us denote by $T_{\mathrm{ext}}=1+\delta T_{\mathrm{ext}}$ the temperature corresponding to the NESS reached when the system is driven with constant intensity ${\chi }_{\mathrm{ext}}=1+\delta {\chi }_{\mathrm{ext}}$. Since the steady value of the excess kurtosis does not depend on the driving intensity, we have only to subtract:

$$\delta T_{\mathrm{ext}}={\left({\chi }_{\mathrm{ext}}\right)}^{2/3}-1=\frac{2}{3}\delta {\chi }_{\mathrm{ext}}+O{\left(\delta {\chi }_{\mathrm{ext}}\right)}^2$$

(17)

from $\delta T$ to make the system homogeneous. Thus, we define:

$$\delta \tilde{T}\equiv \delta T-\delta T_{\mathrm{ext}}.$$

(18)

The homogeneous system for $\delta \tilde{T}$ and $\delta A_2$ reads:

$$\frac{d}{dt}\left(\begin{array}{c}\delta \tilde{T}\\ \delta A_2\end{array}\right)=\left(\begin{array}{cc}-\frac{3}{2}\beta & 1-\beta \\ 3& -2B\end{array}\right)\left(\begin{array}{c}\delta \tilde{T}\\ \delta A_2\end{array}\right).$$

(19)

The eigenvalues $(-{\lambda }_1,-{\lambda }_2)$ and eigenvectors $(v_1,v_2)$ of this system are given by:

$$\begin{array}{ccc}\hfill {\lambda }_1& =\frac{1}{2}\left(\frac{3}{2}\beta +2B+k\right)>0,\hfill & \hfill v_1=\frac{1}{6}\left(\begin{array}{c}2{\lambda }_2-3\beta \\ 6\end{array}\right),\end{array}$$

(20a)

$$\begin{array}{ccc}\hfill {\lambda }_2& =\frac{1}{2}\left(\frac{3}{2}\beta +2B-k\right)>0,\hfill & \hfill v_2=\frac{1}{6}\left(\begin{array}{c}2{\lambda }_1-3\beta \\ 6\end{array}\right),\end{array}$$

(20b)

where we have introduced the parameter:

$$k\equiv {\lambda }_1-{\lambda }_2=\sqrt{{\left(\frac{3}{2}\beta -2B\right)}^2+12(1-\beta )}>0.$$

(21)

With the definitions above, both ${\lambda }_1$ and ${\lambda }_2$ are positive, and ${\lambda }_1>{\lambda }_2$.

Now we can write the solution in both time windows, separately, because $\delta \tilde{T}$ is different over each one. In the first step of the bang-bang, $t\in [0,t_J)$, where $\delta \chi \left(t\right)=\delta {\chi }_{min}$:

$$\left(\begin{array}{c}\delta \tilde{T}\\ \delta A_2\end{array}\right)=\left(\begin{array}{c}\delta T-\frac{2}{3}\delta {\chi }_{min}\\ \delta A_2\end{array}\right)=C_1{v}_1{e}^{-{\lambda }_{1}t}+C_2{v}_2{e}^{-{\lambda }_{2}t}.$$

(22)

In the second step of the bang-bang, $t\in [t_J,t_{\mathrm{f}})$, where $\delta \chi \left(t\right)=\delta {\chi }_{max}$:

$$\left(\begin{array}{c}\delta \tilde{T}\\ \delta A_2\end{array}\right)=\left(\begin{array}{c}\delta T-\frac{2}{3}\delta {\chi }_{max}\\ \delta A_2\end{array}\right)=C_3{v}_1{e}^{-{\lambda }_{1}t}+C_4{v}_2{e}^{-{\lambda }_{2}t}.$$

(23)

The constants $(C_1,C_2,C_3,C_4)$ are obtained by imposing the initial conditions in each time window. For $t=0$, we have the initial condition (12), which determines $C_1$ and $C_2$:

$$C_2=-C_1=3\frac{\delta T_{\mathrm{i}}-\delta T_{min}}{k}=2\frac{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{k}.$$

The point at the final time $t_J$ of the first time window is:

$$\begin{array}{cc}\hfill \delta T_J& =\delta {\tilde{T}}_J+\delta T_{min}=C_1\left[v_1\left(1\right)e^{-{\lambda }_1{t}_J}-v_2\left(1\right)e^{-{\lambda }_2{t}_J}\right]+\frac{2}{3}\delta {\chi }_{min},\hfill \end{array}$$

(24a)

$$\begin{array}{cc}\hfill \delta A_{2J}& =C_1\left(e^{-{\lambda }_1{t}_J}-e^{-{\lambda }_2{t}_J}\right).\hfill \end{array}$$

(24b)

The initial conditions for Equation (23) are supplied by $(\delta T_J,\delta A_{2J})$. Therefore, we can obtain $(C_3,C_4)$ as a function of the switching time $t_J$:

$$\left(\begin{array}{c}\delta T_J-\delta T_{max}\\ \delta A_{2J}\end{array}\right)=C_3{v}_1{e}^{-{\lambda }_1{t}_J}+C_4{v}_2{e}^{-{\lambda }_2{t}_J}.$$

(25)

Note that $(C_3,C_4)$ also depend on the bounds of the driving $(\delta {\chi }_{min},\delta {\chi }_{max})$ through $\delta T_{min}$ and $\delta T_{max}$. By imposing that we have to reach the target state, i.e., Equation (13), we write:

$$\left(\begin{array}{c}-\delta T_{max}\\ 0\end{array}\right)=C_3{v}_1{e}^{-{\lambda }_1{t}_{\mathrm{f}}}+C_4{v}_2{e}^{-{\lambda }_2{t}_{\mathrm{f}}}.$$

(26)

Equations (25) and (26) are four equations for the four unknowns $(C_3,C_4,t_J,t_{\mathrm{f}})$, which thus provide us with the solution to the control problem. Solving for $C_3$ and $C_4$, we obtain:

$$C_3=\frac{2\delta {\chi }_{max}}{k}{\left(\frac{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{e}^{{\lambda }_1{t}_J}-1}{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{e}^{{\lambda }_2{t}_J}-1}\right)}^{\frac{{\lambda }_1}{k}},C_4=-C_3{e}^{-k{t}_{\mathrm{f}}}.$$

(27)

We have introduced the total amplitude of the allowed interval for the driving:

$$\delta {\chi }_{\mathrm{tot}}\equiv \delta {\chi }_{max}-\delta {\chi }_{min}={\chi }_{max}-{\chi }_{min}>0.$$

(28)

In this way, the final time $t_{\mathrm{f}}$ is given as a function of the switching time $t_J$:

$$\begin{array}{cc}\hfill t_{\mathrm{f}}=& \frac{1}{k}ln\left(\frac{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{e}^{{\lambda }_1{t}_J}-1}{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{e}^{{\lambda }_2{t}_J}-1}\right)=t_J+\frac{1}{k}ln\left(\frac{1-\frac{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_1{t}_J}}{1-\frac{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_2{t}_J}}\right),\hfill \end{array}$$

(29)

which is in turn given by the solution of the implicit equation:

$$\begin{array}{cc}\hfill \delta {\chi }_{\mathrm{tot}}& {\left(1-\frac{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_2{t}_J}\right)}^{\frac{{\lambda }_1}{k}}=\delta {\chi }_{max}{\left(1-\frac{\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_1{t}_J}\right)}^{\frac{{\lambda }_2}{k}}.\hfill \end{array}$$

(30)

The set of Equations (29) and (30) provides an analytical solution for the minimum connection time $t_{\mathrm{f}}$ in the CH protocol, which is valid in the linear approximation we are considering in this paper.

3.2. $T_i<T_f$ = 1 Heating-Cooling Bang-Bang

Let us now start from an initial state with $\delta T_{\mathrm{i}}<0$. The analysis of this case is similar to that just carried out for $\delta T_{\mathrm{i}}>0$, but the order of the bangs is reversed. In the first time window, $[0,t_J)$, the maximum driving $\delta {\chi }_{max}$ is applied, whereas in the second time window, $[t_J,t_{\mathrm{f}})$, the minimum driving $\delta {\chi }_{min}$ is applied. Therefore, the homogenisation procedure for the temperature $\delta T$ is also reversed. In the first window $t\in [0,t_J)$ we have:

$$\left(\begin{array}{c}\delta \tilde{T}\\ \delta A_2\end{array}\right)=\left(\begin{array}{c}\delta T-\frac{2}{3}\delta {\chi }_{max}\\ \delta A_2\end{array}\right)=C_5{v}_1{e}^{-{\lambda }_{1}t}+C_6{v}_2{e}^{-{\lambda }_{2}t},$$

(31)

whereas in the second window $t\in [t_J,t_{\mathrm{f}})$ it is:

$$\left(\begin{array}{c}\delta \tilde{T}\\ \delta A_2\end{array}\right)=\left(\begin{array}{c}\delta T-\frac{2}{3}\delta {\chi }_{min}\\ \delta A_2\end{array}\right)=C_7{v}_1{e}^{-{\lambda }_{1}t}+C_8{v}_2{e}^{-{\lambda }_{2}t}.$$

(32)

The initial conditions are given by Equation (12). Inserting them into Equation (31), we obtain:

$$C_5=-C_6=3\frac{\delta T_{max}-\delta T_{\mathrm{i}}}{k}=2\frac{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{k}.$$

(33)

The evaluation of Equation (31) at time $t_J$ gives the initial condition for the second time window. Taking into account the difference in the $\delta \tilde{T}$ variables in Equations (31) and (32), due to our switching the value of the driving intensity at $t=t_J$, one obtains $(C_7,C_8)$ in terms of $t_J$ and also of the bounds $(\delta {\chi }_{min},\delta {\chi }_{max}$) in complete analogy with the CH protocol. In addition, at the final time $t_{\mathrm{f}}$, one must impose that the system reaches the target NESS, i.e., Equation (13), which provides the two extra equations needed to determine the switching time $t_J$ and the connection time $t_{\mathrm{f}}$ as functions of the system parameters. The result is:

$$C_7=\frac{2\delta {\chi }_{min}}{k}{\left(\frac{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{e}^{{\lambda }_1{t}_J}-1}{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{e}^{{\lambda }_2{t}_J}-1}\right)}^{\frac{{\lambda }_1}{k}},C_8=-C_7{e}^{-k{t}_{\mathrm{f}}},$$

(34)

for $(C_7,C_8)$ in terms of $(t_J,t_{\mathrm{f}})$:

$$\begin{array}{cc}\hfill & t_{\mathrm{f}}=\frac{1}{k}ln\left(\frac{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{e}^{{\lambda }_1{t}_J}-1}{\frac{\delta {\chi }_{\mathrm{tot}}}{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{e}^{{\lambda }_2{t}_J}-1}\right)=t_J+\frac{1}{k}ln\left(\frac{1-\frac{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_1{t}_J}}{1-\frac{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_2{t}_J}}\right),\hfill \end{array}$$

(35)

for the minimum connection time in terms of the switching time, and the following implicit equation:

$$\begin{array}{cc}\hfill \delta {\chi }_{\mathrm{tot}}& {\left(1-\frac{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_2{t}_J}\right)}^{\frac{{\lambda }_1}{k}}=-\delta {\chi }_{min}{\left(1-\frac{\delta {\chi }_{max}-\delta {\chi }_{\mathrm{i}}}{\delta {\chi }_{\mathrm{tot}}}{e}^{-{\lambda }_1{t}_J}\right)}^{\frac{{\lambda }_2}{k}},\hfill \end{array}$$

(36)

for $t_J$. In complete analogy with the CH case, the set of Equations (35) and (36) gives the minimum connection time for the HC protocol in the linear approximation. Note that exchanging $\delta {\chi }_{min}\leftrightarrow \delta {\chi }_{max}$ (which entails $\delta {\chi }_{\mathrm{tot}}\to -\delta {\chi }_{\mathrm{tot}}$) leads from Equations (35) and (36) to Equations (29) and (30)—and vice versa.

Above, we have derived analytical expressions for different physical variables of interest, for both the CH and HC bang-bang protocols, in the linear response approximation. More specifically, we have (i) the complete description of the trajectory followed by the system in the phase plane, i.e., the time evolution of the point $(\delta A_2\left(t\right),\delta T\left(t\right))$, and (ii) the switching time $t_J$ and the minimum connection time $t_{\mathrm{f}}$. The linear response approximation has allowed us to obtain analytical predictions as functions of all the relevant physical parameters: not only of the initial temperature $T_{\mathrm{i}}$, as measured by $\delta {\chi }_{\mathrm{i}}=3\delta T_{\mathrm{i}}/2$, but also of the bounds of the driving $({\chi }_{min},{\chi }_{max})$, as measured by $(\delta {\chi }_{min},\delta {\chi }_{max})$. Therefore, it is interesting to inspect the behaviour of the obtained expressions as a function of $(\delta {\chi }_{\mathrm{i}},\delta {\chi }_{min},\delta {\chi }_{max})$ in order to understand the response of the system to the optimal control designed.

4. Trajectories for the Temperature and the Excess Kurtosis

In this section, we look into the trajectories of the temperature and the excess kurtosis, to understand the need for a two-step bang-bang protocol on a physical basis. The time evolution of $\delta T$ and $\delta A_2$ is presented in Figure 1, for both the CH and the HC cases ($\delta T$ solid lines, $\delta A_2$ dashed lines). First, let us analyse the CH protocol (left panel), i.e., $\delta T_{\mathrm{i}}>0$ (also $\delta {\chi }_{\mathrm{i}}>0$). Therein, $\delta T$ relaxes to $\delta T_{min}$ under the action of $\delta {\chi }_{min}$ in the time window $[0,t_J)$ (without reaching it, since the relaxation at constant driving lasts for an infinite time). Simultaneously, $\delta A_2$ starts to increase from its steady value, equal to zero, because Equation (10) implies the following:

$${\left.\frac{d}{dt}\delta A_2\right|}_{t=0^{+}}=-2\delta {\chi }_{min}+3\delta T_{\mathrm{i}}=2\left(\delta {\chi }_{\mathrm{i}}-\delta {\chi }_{min}\right),$$

(37)

which is non-negative because $\delta {\chi }_{min}\le 0\le \delta {\chi }_{\mathrm{i}}$. (Otherwise the connection of the two NESS would be impossible, as rigorously proven in the next section—and in agreement with physical intuition.) A decrease (an increase) in the granular temperature makes the VDF separate from (closer to) the Gaussian shape, i.e., the scaled excess kurtosis $A_2$ correspondingly increases (decreases). Once the target temperature $T_{\mathrm{f}}=1$ is reached inside this first time window, i.e., the temperature curve crosses the horizontal axis $\delta T=0$, the action of $\delta {\chi }_{min}$ cannot be interrupted by setting the thermostat intensity to unity because $\delta A_2>0$ and the system is not in the target NESS. This is why we must let $\delta T$ continue to drop to a value $\delta T_J$ such that $\delta T_J<0$, associated with a kurtosis value $\delta A_{2J}>0$. This point $(\delta A_{2J},\delta T_J)$ is determined by the condition that, at the end of the subsequent relaxation with $\delta {\chi }_{max}$ in the time window $[t_J,t_{\mathrm{f}})$, $\delta A_2$ and $\delta T$ must simultaneously reach their target value (zero). Second, we analyse the HC case (right panel), the discussion is completely analogous and thus we summarise it in the following. In the first time window $[0,t_J)$ with $\delta {\chi }_{max}$, the horizontal axis $\delta T=0$ is crossed at some time smaller than $t_J$, but it is necessary to continue applying $\delta {\chi }_{max}$ to overshoot it, since $\delta A_2<1$ for that time and the system has not reached the target NESS. Once more, the point $(\delta A_{2J},\delta T_J)$ is determined by the condition that, at the end of the second time window with $\delta {\chi }_{min}$, both variables simultaneously vanish.

The need of a two-step protocol, and the order of the bangs, can also be understood—maybe more clearly—by looking at the trajectories in phase space. The trajectories of the phase space $(\delta A_2,\delta T)$ for the CH and the HC cases are shown in Figure 2. The target state is the origin $(0,0)$, so the optimal trajectories must end up thereat. Since the optimal protocols are of bang-bang type, with at most one switch, there are two possibilities: The system approaches the origin following either the heating curve with $\delta {\chi }_{max}$ (red solid line in the left panel) or the cooling curve with $\delta {\chi }_{min}$ (blue solid line in the right panel). These two curves are uniquely defined because the origin is not a fixed point of the evolution equations for $\delta {\chi }_{max}$, nor for $\delta {\chi }_{min}$.

The initial NESS does not lie on either of these two curves—they do not contain any NESS apart from the target state $(0,0)$, thus, the necessity of having a two-step bang-bang is clear. Recall that, for the linear case, there is a theorem ensuring that there is at most one switching.

In Appendix B, we rigorously show that the CH (HC) protocol is the one making it possible to connect the initial NESS with $\delta T_{\mathrm{i}}>0$ ($\delta T_{\mathrm{i}}<0)$.

5. Minimum Connection Time as a Function of the Bounds in the Driving Intensity

This section is devoted to studying the behaviour of the minimum connection time $t_{\mathrm{f}}$ (and also of the switching time $t_J$) as a function of the bounds in the driving intensity. The analysis is carried out for both the CH ($\delta T_{\mathrm{i}}>0$) and HC ($\delta T_{\mathrm{i}}<0$) protocols. We will use the variables without ‘$\delta$’ in order to keep the discussion clearer.

A first question that naturally arises is the range of values of ${\chi }_{min}$ and ${\chi }_{max}$ allowing to connect the initial and target states. In the non-linear regime and in the limit case $({\chi }_{min}=0,$ ${\chi }_{max}=\infty )$, it is always possible to connect two NESS corresponding to temperatures $T_{\mathrm{i}}$ and $T_{\mathrm{f}}$ [22]. However, it is not obvious at all that this is possible when not all the power of the thermostat is available, i.e., in our case with bounds in the driving: ${\chi }_{min}>0$ and ${\chi }_{max}<\infty$. For example, given $T_{\mathrm{i}}>1$, it is unclear whether there appears some change in the behaviour of the connecting time when the upper bound ${\chi }_{max}$ crosses the value ${\chi }_{\mathrm{i}}>1$. Accordingly with our approach throughout, we intend to study this problem within the linear response approximation.

5.1. CH Protocol

First, we consider the CH protocol, $T_{\mathrm{i}}>1$ or $\delta T_{\mathrm{i}}>0$. Figure 3 illustrates the dependence of $t_J$ and $t_{\mathrm{f}}$ on the bounds in the driving. Fixing the value of ${\chi }_{max}$, we can look into their behaviour as functions of ${\chi }_{min}$ (left panel). As the cooling capacity of the thermostat decreases, i.e., as ${\chi }_{min}$ increases, the minimum connection time $t_{\mathrm{f}}$ increases. This is logical, since the class of admissible control functions is being shrunk and the optimal connection thus lasts longer. In addition, the switching time $t_J$ increases: The cooling step of the bang-bang must be longer to compensate for the decrease of cooling power. Both times diverge in the limit as ${\chi }_{min}\to 1^{-}$, where the cooling power of the thermostat is vanishingly small, and thus the cooling step of the bang-bang process takes an infinite time. Now we fix the value of ${\chi }_{min}$ and study the behaviour as functions of ${\chi }_{max}$ (right panel). Analogously, as the heating capacity of the thermostat decreases, i.e., as ${\chi }_{max}$ decreases, $t_{\mathrm{f}}$ increases, because the class of admissible controls becomes smaller. On the other hand, the behaviour of $t_J$ is reversed, and $t_J$ increases with ${\chi }_{max}$. This is also logical, since the first step of the bang-bang is the cooling one, and as the heating capacity of the thermostat is increased, the cooling step must take a longer time. In this case, it is only $t_{\mathrm{f}}$ that diverges in the limit as ${\chi }_{max}\to 1^{+}$. The lack of heating capacity makes the duration of the second (heating) step diverge, since the time needed to relax towards $T_{\mathrm{f}}=1$ is infinity for a constant value of the driving $\chi ={\chi }_{\mathrm{f}}=1$. There is no change in behaviour in the connection time when ${\chi }_{max}$ crosses the value ${\chi }_{\mathrm{i}}$, the driving intensity corresponding to the initial value of the temperature. This is neatly observed in the inset, where a zoom of the graph for drivings ${\chi }_{max}\in [1,{\chi }_{\mathrm{i}}]$ is plotted.

An important point is the divergence of the connection time as ${\chi }_{min}\to 1^{-}$ (for fixed ${\chi }_{max}$) and as ${\chi }_{max}\to 1^{+}$ (for fixed ${\chi }_{min}$). Therefore, if ${\chi }_{\mathrm{f}}=1$ lies outside the interval $[{\chi }_{min},{\chi }_{max}]$, the target NESS is unreachable. In other words, the bounds in the driving must verify $({\chi }_{max}\ge 1,{\chi }_{min}\le 1)$, i.e., $(\delta {\chi }_{max}\ge 0,\delta {\chi }_{min}\le 0)$, to make it possible to connect the initial and target states. In other words, $T_{\mathrm{f}}$ must belong in the interval $[T_{min},T_{max}]$. In fact, it is possible to rigorously show that the connecting time only diverges when either ${\chi }_{min}\to 1^{-}$ or ${\chi }_{max}\to 1^{+}$, see Appendix C for details.

5.2. HC Protocol

Figure 4 illustrates the situation for the HC protocol ($T_{\mathrm{i}}<1$ or $\delta T_{\mathrm{i}}<0$). Note that the panels are basically the horizontal reflections of those in Figure 3, with the roles of ${\chi }_{min}$ and ${\chi }_{max}$ being exchanged. Therefore, the line of reasoning for physically understanding the observed behaviours is completely similar to the one in the previous section, and it will not be repeated here. We would only like to highlight the increase in the minimum connection time as the bounds become tighter, due to the shrinking of the set of admissible control functions, and its divergence for ${\chi }_{min}\to 1^{-}$ (fixed ${\chi }_{max}$) and ${\chi }_{max}\to 1^{+}$ (fixed ${\chi }_{min}$), which marks the impossibility of reaching a target state with temperature $T_{\mathrm{f}}=1$ lying outside the interval $[T_{min},T_{max}]$.

6. Validity of the Linear Response Approximation

The results obtained and analysed in the previous sections are quite general. On the one hand, we have derived expressions for the relevant physical quantities as functions of the bounds in the driving intensity $({\chi }_{min},{\chi }_{max})$ (or $(\delta {\chi }_{min},\delta {\chi }_{max})$). On the other hand, the linear response approximation limits the results, since we have assumed that the system remains always close to the target NESS. Therefore, it is relevant to investigate the possible validity of our results beyond the strictly linear framework.

In Ref. [22], it was shown that the minimum connection time in the non-linear case—for a full-strength thermostat $0\le \chi <\infty$, which we denote here by $t_{\mathrm{f}}^{n\ell }$—is given by:

$$t_{\mathrm{f}}^{n\ell }\sim {\left(\frac{2}{3B}\right)}^{1/2}|\delta {\chi }_{\mathrm{i}}{|}^{1/2},\left|\delta {\chi }_{\mathrm{i}}\right|\ll 1,$$

(38)

when the initial and final states are close—as expressed by the condition $|\delta {\chi }_{\mathrm{i}}|\ll 1$ (It must be remarked that the non-dimensionalisation of time in Ref. [22], $t^{*}={\zeta }_0{T}_{\mathrm{i}}^{1/2}t$, differs from ours in Equation (7) by a factor $\sqrt{T_{\mathrm{i}}/T_{\mathrm{f}}}$. This factor does not affect the lowest order asymptotic expression in Equation (38), since the introduced corrections are higher-order). Actually, Equation (38) does not have to hold for the linear case developed in this paper because we are considering that the driving intensity $\chi$ is restricted to a small interval, ${\chi }_{min}=1+\delta {\chi }_{min}\le \chi =1+\delta \chi \le {\chi }_{max}=1+\delta {\chi }_{max}$. (Recall that $\delta {\chi }_{max}\ge 0$, whereas $\delta {\chi }_{min}\le 0$.) Notwithstanding, we may progressively separate the bounds from unity and compare our linear response predictions with Equation (38). More specifically, it is interesting to take the limit ${\chi }_{min}\to 0$ and ${\chi }_{max}\to \infty$ and analyse the possible convergence of our minimum connection time (for both the CH and HC cases) to the time given by Equation (38).

In order to further explore this possible convergence, we have approximated $t_{\mathrm{f}}$ to first order in $\delta {\chi }_{\mathrm{i}}$. In Appendix E, it is shown that for short connecting times $t_{\mathrm{f}}\ll 1$, one has:

$$t_{\mathrm{f}}\sim {\left(\frac{2}{{\lambda }_1{\lambda }_2}\right)}^{1/2}{\left(\frac{\delta {\chi }_{\mathrm{tot}}\left|\delta {\chi }_{\mathrm{i}}\right|}{-\delta {\chi }_{min}\delta {\chi }_{max}}\right)}^{1/2},$$

(39)

where

$${\lambda }_1{\lambda }_2=3\beta B+3(\beta -1)=3B\left(1+\frac{3}{16}{a}_2^{\mathrm{s}}\right)+\frac{9}{16}{a}_2^{\mathrm{s}}.$$

(40)

Note that Equation (39) is valid to the lowest order in $t_{\mathrm{f}}$—terms of the order of $t_{\mathrm{f}}^2$ have been neglected—but no assumption has been made with regard to $\delta {\chi }_{min}$ and $\delta {\chi }_{max}$.

In order to make the comparison between the non-linear (with full-strength thermostat) and linear (with bounds in the driving) expressions above, we have represented the minimum connecting time for different values of the bounds ${\chi }_{min}$ and ${\chi }_{max}$ in Figure 5 and Figure 6—for the CH and HC cases, respectively. Therein, we show the linear response expressions Equation (29) (CH case) and (35) (HC case), together with the approximate linear expression (39), and the non-linear expression (38) for a full-strength thermostat. We observe how the times given by Equations (29) and (35), as well as their approximations (39), rapidly converge to Equation (38) as the bounds separate from unity. This convergence is qualitatively similar in the CH and HC cases. There are no significant differences between them up to this point.

Let us look at the convergence towards the non-linear expression (38) in more detail. First, we consider the CH case in Figure 7, which can be seen as a zoom of Figure 5—for values of the bounds such that the linear time is close to the non-linear one. It is clearly observed that as the bounds of the driving separate from unity, the linear response prediction approaches the non-linear expression (38) “from above”: The connection times of the linear theory are longer than those for the non-linear case. This is consistent, since Equation (38) was obtained for the largest possible set of control functions, i.e., $({\chi }_{min}=0,{\chi }_{max}=\infty )$: The loosest the restrictions on the control functions are, the shortest the minimum connection time is.

Now we have a closer look at the HC case in Figure 8. As the bound in the controls move away from unity, the minimum connection time is also very close to the non-linear expression (38). However, the convergence “from above” observed in the CH case is broken. In fact, the linear prediction is neatly below the non-linear one for the data corresponding to the most extreme values of the bounds. This marks a first physical limit for the range of controls that can be used in the linear approach: Beyond the values ${\chi }_{min}$ and ${\chi }_{max}$ such that the linear prediction for the minimum connection time become smaller than that provided by the non-linear prediction (38), the linear theory is clearly not valid. Recall that the latter was obtained for the full strength of the thermostat, $({\chi }_{min}=0,{\chi }_{max}=\infty )$, so for a smaller set of controls, the minimum connection time must be longer.

We have illustrated the breakage of the convergence “from above” in the HC protocol for the particular case $\alpha =0.9$ and $d=2$. This behaviour is robust: It occurs for all $\alpha$, and also for $d=3$. On the other hand, in the CH protocol, the inversion of the natural convergence “from above” never comes about. This asymmetry between the CH and HC protocols stems from the physical limit that ${\chi }_{min}$ has: While ${\chi }_{max}$ can be as large as desired, ${\chi }_{min}$ must always be non-negative. This entails that, when applying a CH protocol, the granular temperature T evolves between the values 0 and 1 for all times, which prevents the system from presenting important deviations from the linear response behaviour. However, for the HC protocol, the temperature can reach arbitrarily large values under the action of a high enough driving ${\chi }_{max}$, which makes the linear response approximation no longer valid. In fact, if we had studied the system from a purely mathematical point of view and removed the physical restriction ${\chi }_{min}\ge 0$ (letting it vary between $-\infty$ and $+\infty$), this asymmetry between the CH and HC cases would have disappeared.

It is interesting to remark that as $\alpha$ decreases (i.e., as the inelasticity increases), less extreme values of ${\chi }_{max}$ and ${\chi }_{min}$ are needed to provoke the inversion. In other words, the linear approximation breaks down for less extreme bounds. For example, let us consider $\alpha =0.1$, which can be regarded as a high-inelasticity case—as opposed to the low-inelasticity case $\alpha =0.9$. Fixing ${\chi }_{min}=0$, the connection times of the linear approximation become shorter than those given by Equation (38) for $\delta {\chi }_{max}\ge 9.2$ in the range of $\delta T_{\mathrm{i}}\in [-0.01,0]$ (or, equivalently, $\delta {\chi }_{\mathrm{i}}\in [-0.015,0]$), smaller than the value $\delta {\chi }_{max}\ge 9.8$ for $\alpha =0.9$. This trend with $\alpha$ of the bounds leading to the inversion of the convergence “from above” can be understood by recalling that $a_2^{\mathrm{s}}$ is a decreasing function of $\alpha$. Consequently, the importance of the heating term in the evolution equation of the temperature (8a), $\chi \left(1+\frac{3}{16}{a}_2^{\mathrm{s}}\right)$, increases as $\alpha$ is lowered: A smaller value of $\chi$ is needed to obtain the same value of the heating term.

7. Discussion

Our work improves the understanding of the optimal control of driven granular gases. The results obtained in this paper complement and enrich those obtained in Ref. [22] for a full-strength thermostat. The inclusion of bounds in the driving, ${\chi }_{min}\le \chi \le {\chi }_{max}$ raises non-trivial questions that have been answered by our study, such as the range of initial and target temperatures that can be connected. Our investigation has been carried out in the linear response regime, i.e., the initial and target states are close enough—and so are the bounds of the driving ${\chi }_{min}$ and ${\chi }_{max}$. This allows us to linearise the evolution equations around the final (target) NESS.

The linear response approximation leads to a set of evolution equations that are linear both in the control function and the dynamical variables—more precisely, in their deviations from their target values. Therefore, we obtain a linear control problem that can be completely solved. A rigorous mathematical theorem ensures that the optimal control, minimising the connection time, is of bang-bang type with at most one switching: i.e., the optimal control comprises two time windows $[0,t_J)$ and $[t_J,t_{\mathrm{f}})$, with the control being equal to one of its limiting values, either ${\chi }_{min}$ or ${\chi }_{max}$, in the first time window $[0,t_J)$ and changing to the other limiting value at the switching time $t_J$. Therefore, two types of bang-bang protocols arise, depending on the order of the bangs: ${\chi }_{min}$ followed by ${\chi }_{max}$, which we have termed CH, or ${\chi }_{max}$ followed by ${\chi }_{min}$, which we have termed HC. We have shown that the CH protocol is the optimal one when the initial temperature $T_{\mathrm{i}}$ is larger than the final one $T_{\mathrm{f}}$, whereas the HC protocol is the optimal one in the reverse situation, $T_{\mathrm{i}}<T_{\mathrm{f}}$.

We have investigated the behaviour of the connection time as a function of the bounds in the driving intensity. This study has allowed us to elucidate the range of initial and final temperatures that can be connected. We have shown that the final temperature has to lie between the temperatures $T_{min}$ and $T_{max}$, where $T_{min}$ ($T_{max}$) is the steady temperature corresponding to the constant driving ${\chi }_{min}$ (${\chi }_{max}$). On the other hand, the initial temperature may lie outside the interval $[T_{min},T_{max}]$, and the connection is still possible: The minimum connection time is still finite when the upper bound ${\chi }_{max}$ crosses the initial temperature $T_{\mathrm{i}}>T_{\mathrm{f}}$ (or the lower bound crosses the initial temperature $T_{\mathrm{i}}<T_{\mathrm{f}}$) (It could be argued that, still, the most relevant physical situation corresponds to the case $T_{\mathrm{i}}\in [T_{min},T_{max}]$ because one needs to prepare the system in the initial NESS).

In addition, we have explored the limits of the validity of the linear response approximation we have employed throughout. We have done this by loosening the restrictions on the bounds ${\chi }_{min}$ and ${\chi }_{max}$. Specifically, we have analysed the behaviour of our (linear response) prediction for the minimum connection time, $t_{\mathrm{f}}$, as ${\chi }_{min}$ is decreased to very small values, and ${\chi }_{max}$ is increased to very large values. This behaviour has been compared with the minimum connecting time for the non-linear case $t_{\mathrm{f}}^{n\ell }$, which was obtained when the thermostat has its full strength, ${\chi }_{min}=0$ and ${\chi }_{max}=\infty$. Specifically, we have compared the linear time with the asymptotic expression for $t_{\mathrm{f}}^{n\ell }$ for small temperature jumps—in which the corresponding connection times are also very small [22]. For the CH case, we have found that $t_{\mathrm{f}}$ tends to $t_{\mathrm{f}}^{n\ell }$ always “from above”, $t_{\mathrm{f}}>t_{\mathrm{f}}^{n\ell }$. This is logical, since the largest set of controls—such as that of the full-power thermostat—should lead to the shortest connection times. However, for the HC case, we have found that the tendency from above towards $t_{\mathrm{f}}$ is broken for large enough values of ${\chi }_{max}$. This marks a limit of validity for the linear response approximation in this case. The asymmetry between the CH and HC protocols can be understood on a physical basis: In the latter case, heating precedes cooling, and thus the temperature departs from the vicinity of $T_{\mathrm{f}}$ for high enough ${\chi }_{max}$—whereas in the former, cooling precedes heating, and the system remains closer to the target state even when ${\chi }_{max}$ becomes large.

In order to further look into the behaviour described in the previous paragraph, we have looked into the regime of short connecting times $t_{\mathrm{f}}\ll 1$ within the linear response framework. Note that our linear response predictions for $t_{\mathrm{f}}$, as given by Equations (29) and (35), contain all the powers of $\delta {\chi }_{\mathrm{i}}/\delta {\chi }_{\mathrm{tot}}$, $\delta {\chi }_{max}/\delta {\chi }_{\mathrm{tot}}$, $\delta {\chi }_{min}/\delta {\chi }_{\mathrm{tot}}$. Linear response assumes that both $\delta {\chi }_{\mathrm{i}}\ll 1$, $\delta {\chi }_{max}\ll 1$, and $\delta {\chi }_{min}\ll 1$, but the ratios between one another are in principle of the order of unity. It is only when $\delta {\chi }_{\mathrm{i}}$ is much smaller than $\delta {\chi }_{max}$ and $\delta {\chi }_{min}$ that the connection time becomes small. In this regime, we have obtained a simple approximate expression for $t_{\mathrm{f}}$ valid to the lowest order in $\delta {\chi }_{\mathrm{i}}$, which also depends on the bounds in the driving. This approximate expression always gives connection times that are longer than that for the full-power thermostat, both for the CH and HC cases. This means that the inversion of the tendency “from above” towards $t_{\mathrm{f}}^{n\ell }$ comes from higher-order terms in the ratios $\delta {\chi }_{\mathrm{i},max,min}/\delta {\chi }_{\mathrm{tot}}$.

Our work also opens the door to finding new optimal controls for other non-equilibrium systems. For example, let us look at a colloidal particle moving in the vicinity of a minimum of the trapping potential—which can be thus considered to be harmonic. In that case, the temperature of the thermal bath in which the particle is immersed plays the role of the driving intensity. Interestingly, the temperature of the bath can be effectively increased by adding a random forcing that can be modelled as a Gaussian white noise [40,41]. In this way, the effective temperature changes from $T_{min}$ (room temperature) to $T_{max}$ (thousands of kelvins). The similitude of the mathematical framework, linear evolution equations and bounded control, makes it appealing to analyse the optimal connection—also in the sense of minimising the connection time—in that case and compare the corresponding results with those derived here.

It is also relevant to extend the optimal controls derived here—and also of those for a full-power thermostat in Ref. [22]—to more complex situations, such as rough granular gases [60,61,62,63,64]. Translational and rotational velocities are correlated, even for small roughness [60,61]. In addition, the number of relevant variables increases: The values of the translational and rotational temperatures are, in general, different, and additional cumulants emerge (See, for example, [62,63,64] for the general kinetic framework). The impact of these features on the optimal controls is thus an interesting prospect for future work.