Kovacs-like memory effect in athermal systems: linear response analysis

We analyse the emergence of Kovacs-like memory effects in athermal systems within the linear response regime. This is done by starting from both the master equation for the probability distribution and the equations for the physically relevant moments. The general results are applied to a general class of models with conserved momentum and non-conserved energy. Our theoretical predictions, obtained within the first Sonine approximation, show an excellent agreement with the numerical results.


Introduction
The equilibrium state of physical systems is characterised by the value of a few macroscopic variables, for example pressure, volume and temperature in fluids. This characterisation of the equilibrium state is complete, in the sense that different samples sharing the same values of the macroscopic variables respond identically to an external perturbation. On the contrary, a system in a nonequilibrium state, even if it is stationary, is not completely characterised by the value of the macroscopic variables: the response to an external perturbation depends also on additional variables or, equivalently, on its entire thermal history. Therefore, it is often said that the response depends on the way the system has been aged.
A pioneering work in the field of memory effects in nonequilibrium systems was carried out by Kovacs [1]. The Kovacs experiment showed that pressure, volume and temperature did not completely characterise the state of a sample of polyvinyl acetate that had been aged for a long time at a certain temperature T 1 . The pressure was fixed during the whole experiment and the time evolution of the volume was recorded. At a certain time t w , the temperature was suddenly changed to T, for which the equilibrium value of the volume equalled its instantaneous value at t w . Counterintuitively (thinking in equilibrium terms), the volume did not remain constant. Instead, it displayed a hump, passing through a maximum before tending back to its equilibrium (and initial) value.
We look into the Kovacs experiment in a more detailed way in Figure 1. In recent studies of the effect in glassy systems, the relevant physical variable is the energy instead of the volume [2][3][4][5][6][7][8]. The system is equilibrated at a "high" temperature T 0 and at t = 0 the temperature is suddenly quenched to a lower temperature T, after which the relaxation function φ(t) of the energy E is recorded. Specifically, where E e is the average equilibrium energy at temperature T. Then, a similar procedure is followed, equilibrating the system again at T 0 , but at t = 0 the temperature is changed to an even lower value T 1 , T 1 < T < T 0 . The system relaxes isothermally at T 1 for a certain time t w , such that E (t = t w ) equals E e . At this time t w , the temperature is increased to T but the energy does not remain constant: it displays the behaviour that is qualitatively shown by K(t). At first, K(t) increases from zero until a maximum is attained for t = t k , and only afterwards it goes back to zero. Similarly to Figure 1. Scheme of the Kovacs experiment described in the text. The dashed curve on the right, labelled by φ(t), represents the direct relaxation from T 0 to T. The dashed curve on the left stands for the part of the relaxation from T 0 to T 1 that is interrupted at t = t w by the second temperature jump, changing abruptly the temperature from T 1 to T. After this second jump, the system follows the non-monotonic response K(t) (solid line), which reaches a maximum at t = t k and, afterwards, approaches φ(t) for very long times. the relaxation function, we have defined K(t) = E(t) − E e , for t ≥ t w . Note that K(t) ≤ φ(t) for all times, with the equality being only asymptotically approached in the long time limit.
For molecular (thermal) systems, the equilibrium distribution has the canonical form and it has been shown that in linear response theory [6] where the final temperature T and the waiting time t w are related by In linear response, the relaxation function φ(t) is proportional to the equilibrium time correlation function E(0)E(t) e − E 2 e (fluctuation-dissipation theorem), which decays monotonically in time [9]. The linear response results above make it possible to understand the crux of the observed Kovacs hump in experiments [6]: (i) the inequality 0 ≤ K(t) ≤ φ(t), which assures that the hump always has a positive sign (from now on, "normal" behaviour), (ii) the existence of only one maximum in the hump and (iii) the increase of the maximum height and the shift of its position to smaller times as t w is decreased. Nevertheless, it must be noted that the experiments, both real [1] and numerical [2][3][4][5]7] are mostly done out of the linear response regime: thus, it seems that the validity of these results extends beyond expectations. In fact, it has been checked in simple models that the linear approximation still gives a fair description of the hump for not-so-small temperature jumps [8].
Very recently, the investigation of Kovacs-like effects in athermal systems has been started. A granular fluid provides a prototypical example of an athermal system, which is intrinsically out-of-equilibrium [10]. A physical mechanism that inputs energy into the system and balances in average the energy loss in collisions, for instance the so-called stochastic thermostat [11], must be considered to reach a nonequilibrium steady state (NESS). Moreover, in general, fluctuations are far more important in granular systems than in molecular systems because of their smallness. The number of particles N, ranging from 10 2 to 10 4 , is large enough to make it possible to apply the methods of statistical mechanics but definitely much smaller than Avogadro's number.
The simplest case is that of uniformly driven granular gases considered in [12,13]. The value of the kinetic energy (granular temperature T g ) at the NESS is controlled by the driving intensity ξ of the stochastic thermostat. Therefore, a Kovacs-like protocol can be implemented in a completely analogous way to the one described above, with the changes T → ξ, E → T g . One of the main differences found is the emergence of "anomalous" Kovacs behaviour for large enough dissipation, when K(t) becomes negative and displays a minimum instead of a maximum. It must be stressed that these results have been obtained in the nonlinear regime, that is, for driving jumps ξ 0 − ξ, ξ 0 − ξ 1 that are not small.
More recently, Kovacs-like behaviour has been observed in other, more complex, athermal systems. This is the case of disordered mechanical systems [14] and also of active matter [15]. In the latter, a "giant" Kovacs hump has been reported, in the sense that the numerically observed maximum is much larger than the one predicted by the extrapolation of the linear approximation expression (1) to the considered protocol. Moreover, an alternative derivation of (1) has been provided in the supplemental material of [15]. This derivation holds for athermal systems, since it does not make use of either the explicit form of the probability distribution or the relationship between response functions and time correlations at the steady state, but is restricted to discrete-time dynamics at the macroscopic (average) level of description.
The objectives of our paper are twofold. Firstly, we put forward a rigorous and general derivation of the linear response expression for the Kovacs hump for systems with a realistic continuous time dynamics. This is done at both the mesoscopic and macroscopic levels of description, starting from the master equation for the probability distribution and from the hierarchy of equations for the moments, respectively. Our proof is also valid for athermal systems, since no hypothesis is needed with regard to the form of either the stationary probability distribution or the fluctuation-dissipation relation. Secondly, we apply our results to a simple class of dissipative models that mimic the shear component of a granular fluid, recently introduced [16][17][18]. Therein, we obtain explicit expressions for the Kovacs hump within the first Sonine approximation and compare the theoretical predictions with numerical results. We also discuss the compatibility of the non-monotonic behaviour displayed by the granular temperature in the Kovacs experiment and the monotonic approach to the NESS of the nonequilibrium entropy or H-functional [19][20][21].

General Markovian dynamics
Here, we consider a general system, whose state is completely characterised by a vector x with M components, x = {x 1 , x 2 , . . . , x M }. For example, in a one-dimensional Ising chain of N spins x i = σ i = ±1, and M = N; for a gas comprising N particles with positions r i and velocities v i , M = 6N and x = {r 1 , v 1 , . . . , r N , v N }. As is customary and for the sake of simplicity, from now on we use a notation suitable for systems in which the states can be labelled with a discrete index α, 1 ≤ α ≤ Ω. For example, this is the case of the Ising system above, where Ω = 2 N . The generalisation for a continuous index is straightforward, by changing sums into integrals and Kronecker delta by Dirac delta [9].
At the mesoscopic level of description, we assume that x is a Markov process and its dynamics is governed by a master equation for the probabilities P(x α , t), where W(x α |x α ; ξ) are the transition rates from state x α to state x α , and our notation explicitly marks their dependence on some control parameter ξ. Equation (3) can be formally written as where |P (t) is a vector (column matrix) whose components are the probabilities P(x α , t) and W(ξ) is the linear operator (square matrix) that generates the dynamical evolution of |P (t) , Let us assume that the Markovian dynamics is ergodic (or irreducible [9]), that is, all the states are dynamically connected through a chain of transitions with non-zero probability. Therefore, there is a unique steady solution of the master equation |P s (ξ) , which verifies and depends on the parameter ξ controlling the system dynamics. Note that ergodicity does not imply detailed balance, so we can have non-zero currents in the steady state. In general, this means that the system approaches a NESS in the long time limit, not an equilibrium state. Now, we consider the system evolving from certain initial state at time t 0 , characterised by the distribution |P (t 0 ) . The formal solution of the master equation can be written as This is the starting point for our derivation of the expression for the Kovacs effect in the linear response approximation, which is carried out in the next section. The time evolution of any physical property Y is obtained right away. Let us denote the value of Y for a given configuration x by Y(x): its expected or average value is given by where |Y is a ket whose components are Y(x α ), and Y | its corresponding bra (row matrix with the same components 1 ). By substituting (7) into (8), it is obtained that in which Y s (ξ) is the average value at the steady state of the system corresponding to ξ. Now we investigate the relaxation of the system from the steady state for ξ 0 = ξ + ∆ξ to the steady state for ξ. We do so in linear response, that is, ∆ξ is considered to be small and we neglect all terms beyond those linear in ∆ξ. Thus, at t = 0 we have that the probability distribution is Substitution of (10) into (9) yields the formal expression for the relaxation of Y in linear response, In order to have an order of unity function, one may define a normalised relaxation function 1 We are assuming that Y is a real quantity for all the configurations.
Sometimes, the relaxation function is further normalised by considering φ Y (t)/φ Y (t = t 0 ), see for instance [6,22], but such a multiplying factor is clearly not physically relevant and will not be introduced here.

The Kovacs protocol. Linear response analysis from the master equation.
Here, it is a Kovacs-like protocol that we introduce, by considering that the parameter ξ controlling the dynamics is changed in the following stepwise manner: Therefore, since ξ 0 is kept for an infinite time, at t = 0 the system is prepared in the corresponding steady state, P (t = 0) = P s (ξ 0 ). Our idea is to consider that the jumps ξ 1 − ξ 0 and ξ − ξ 1 are small, in the sense that all expressions can be linearised in the magnitude of these jumps. Note that this protocol is completely analogous to the Kovacs protocol described in the introduction, see Figure 1, but with ξ playing the role of the temperature.
We start by analysing the relaxation in the first time window, 0 < t < t w . Therein, we apply (7) with the substitutions t 0 → 0 and ξ → ξ 1 , that is, The final distribution function, at t = t w , is the initial condition for the next stage, t > t w , in which the system relaxes towards the steady state corresponding to ξ. Making use again of (7) with t 0 → t w , It must be stressed that the above expressions, Equations (14) and (15), are exact, no approximation has been made. The linear response approximation is introduced now: we assume that both jumps ξ 0 − ξ 1 and ξ − ξ 1 are small, so we can expand both |P s (ξ 0 ) − P s (ξ 1 ) and |P s (ξ 1 ) − P s (ξ) , similarly to what was done in Equation (10). Namely, Note that in both (16a) and (16b), the derivatives are evaluated at ξ; the difference introduced by evaluating them at either ξ 1 or ξ 0 are second order in the deviations. Then, This equation can be further simplified: since the two terms on its rhs are first order in the jumps, we can substitute W(ξ 1 ) with W(ξ), with the result This is the superposition of two responses: the first term on the rhs gives the relaxation from ξ 0 to ξ 1 , which starts at t = 0, whereas the second term stands for the relaxation from ξ 1 to ξ, which starts at t = t w . We have chosen to write −(ξ − ξ 1 ) in the second term because ξ > ξ 1 in the Kovacs protocol. The same structure in Equation (18) is transferred to the average value. Taking into account Equation (9), in which we recognise the relaxation function in linear response, defined in Equation (12). We can also normalise the response in this experiment, by defining a function K(t) as follows, It is understood that, as ξ 0 − ξ → 0, both prefactors ξ 0 −ξ 1 ξ 0 −ξ and ξ−ξ 1 ξ 0 −ξ are kept of the order of unity. A few comments on Equation (20) are in order. Hitherto, no restriction has been imposed on the state of the system at t = t w ; therefore, (20) is valid for arbitrary (ξ 0 , ξ 1 , ξ), provided that the jumps are small enough and the ratio of the jumps is of the order of unity. The function K(t) corresponds to a Kovacs-like experiment when ξ is chosen as a function of t w in such a way that Alternatively, one may consider that (21) defines t w as a function of ξ.
The complete analogy between Equations (20) and (21) and Equations (1) and (2) is apparent. Nevertheless, we have made use neither of the explicit form of the steady state distribution (in general non-canonical) nor of the relation between response functions and time correlation functions (fluctuation-dissipation relation), which were necessary in [6] to demonstrate Equation (1). Therefore, the proof developed here is more general, being valid for any steady state, equilibrium or nonequilibrium, and thus it specifically holds in athermal systems. Also, it must be noted that it can be easily extended to the Fokker-Planck, or the equivalent Langevin, equation.

Linear response from the equations for the moments
In this section, we do not start from the equation for the probability distribution, but from the equations for the relevant physical properties of the considered system. For example, one may think of the hydrodynamic equations for a fluid or the law of mass action equations for chemical reactions. Of course, these equations can be derived in a certain "macroscopic" approximation [9], which typically involves neglecting fluctuations, from the equation for the probability distribution by taking moments, but this is not our approach here. Anyhow, we borrow this term to call our starting point "equations for the moments".
We denote the relevant moments by z i , i = 1, . . . , J, where J is the number of relevant moments. The equations for the moments have the general form where f i are continuous, in general nonlinear, functions of the moments. This is a key difference between moment equations and the master (or Fokker-Planck) equation, since the latter is always linear in the probability distribution. Therefore, unlike the master equation, Equation (22) cannot be formally solved for arbitrary initial conditions. Notwithstanding, in the linear response approximation, we show here that a procedure similar to the one carried out in the previous section leads to the same expression for the Kovacs hump. We assume that there is only one steady solution of Equation (22) that is globally stable, at which the corresponding values of the moments are z s i (ξ). Linearisation of the dynamical system around the steady state gives We are using a notation completely similar to that in the previous section, |z is a vector, represented by a column matrix with components z i , and M(ξ) is a linear operator, represented by a square matrix with elements Note that the dimensions of these matrices are much smaller than those for the master equation, since J is of the order of unity and does not diverge in the thermodynamic limit. In general, M ij = M ji and the operator M is not Hermitian. However, we do not need M to be Hermitian to solve the linearised system in a formal way, as shown below. Analogously to what was done for the master equation, the formal solution of Equation (23) is In particular, if the initial condition is chosen to correspond to the steady state for ξ 0 = ξ + ∆ξ, one has The response for any of the relevant moments can be extracted by projecting the above result onto the "natural" basis |u i , whose components are u ij = δ ij . Then, the normalised linear response function for z i can be defined by Note the utter formal analogy of expression (27) with Equation (12), which was obtained from the master equation. The proof of the expression for the Kovacs hump follows exactly the same line of reasoning, and the result is exactly that in Equations (20) and (21), and thus it is not repeated here.

Definition of the model. Kinetic description
Here, we briefly put forward a class of models for the shear component of the velocity in a granular gas, focusing on the features that are needed for the discussion of memory effects. A more complete description of the model, including its physical motivation, can be found in [16][17][18]23].
The system is defined on a 1d lattice: there is a particle with velocity v l at each site l. The system configuration is thus given by v ≡ {v 1 , ..., v N }. The dynamics proceeds through inelastic nearest-neighbour binary collisions: each pair (l, l + 1) collides inelastically with a characteristic rate proportional to |v l − v l+1 | β , with β ≥ 0. For β = 0, nearest-neighbour particles collide independently of their relative velocity (the so-called Maxwell-molecule model [24]), whereas for β = 1 and β = 2 we have a collision rate analogous to that of hard spheres and very hard spheres, respectively [10,25]. The pre-collisional velocities are transformed into the post-collisional ones by the operatorb l , where 0 < α ≤ 1 is the normal restitution coefficient. Momentum is always conserved in collisions, with equality holding only in the elastic limit α = 1.
In order to have a steady state, a mechanism that injects energy into the system, and thus compensates in average the energy loss in collisions, must be introduced. For the sake of simplicity, we consider here that the system is "heated" by a white noise force, that is, the so-called stochastic thermostat [11,26]. The velocity change introduced by the stochastic forcing is and ξ i (τ) are Gaussian white noises of amplitude χ for i, j = 1, . . . , N. This version of the stochastic thermostat conserves total momentum, which is needed to assure the existence of a well-defined steady state [27,28]. We do not write here the master-Fokker-Planck equation for the N-particle P N (v, τ) probability density of finding the system in state v at time τ. Instead, we directly write the "kinetic" equation for the one-particle distribution function, namely . Therefrom, all the one-site velocity moments, which are the relevant physical quantities for our present purposes, can be derived, v n l (τ) ≡ +∞ −∞ dv v n P 1 (v; l, τ). As usual in kinetic theory, a closed equation for P 1 can be written only after introducing the Molecular Chaos assumption: spatial correlations at different sites are of the order of N −1 and then negligible.
We analyse here homogeneous states; if the initial state is homogeneous, the above dynamics preserves homogeneity. Moreover, we employ the usual notation in kinetic theory, f (v, τ) ≡ P 1 (v, £ l; τ) and consider the quasi-elastic limit, that is, ≡ 1 − α 2 1. This allows us to write a simpler expression for the inelastic collision term. Proceeding along the same lines as in [17,21], it is obtained that [29] We can define a dimensionless time scale, t = ω τ, over which where ξ is the rescaled strength of the noise, ξ = χ ω . Note that this kinetic equation is, like Boltzmann's or Enskog's, nonlinear in f (v, t).
The main physical magnitude is the granular temperature T, which we define as We used the notation T g for the granular temperature in the introduction, to differentiate it from the usual thermodynamic temperature T. Since the latter plays no role in our system, we employ the usual notation T for the granular temperature hereafter. It is also customary to define the thermal velocity Taking moments in (32) and making the change of variables above, one gets The first term on the rhs stems from collisions and cools the system, in the sense that it always make the granular temperature decrease. The second term stems from the stochastic thermostat and heats the system and, thus, in the long time limit a NESS is attained in which both terms counterbalance each other.

First Sonine approximation
The equation for the granular temperature is not closed in general and then an expansion in Sonine (or Laguerre) polynomials is typically introduced, where L (m) k (x) are the associated Laguerre polynomials [30]. In kinetic theory, m = d 2 − 1, with d being the spatial dimension, and often the notation S k (x) ≡ L ( d 2 −1) k is used. Here, we mainly use the so-called first Sonine approximation, in which (i) only the term with k = 2 is retained and (ii) nonlinear terms in a 2 are neglected. The coefficient a 2 is the excess kurtosis, c 4 = 3(1 + a 2 )/4.
Although the linearisation in a 2 is quite standard in kinetic theory, we derive firstly the evolution equations considering just step (i) of the first Sonine approximation, that is, we truncate the expansion for the scaled distribution (36) after the k = 2 term. Henceforth, we call this approximation, nonlinear first Sonine approximation. Afterwards, in the numerical results, we will discuss how both approximations, nonlinear and standard, give almost indistinguishable results.
In the nonlinear first Sonine approximation, it is readily obtained the evolution equation for the temperature [29] where ζ 0 = π − 1 /2 2 1+β Γ 3+β 2 . Unless β = 0 (Maxwell molecules), the equation for the temperature is not closed. Then, we write down the equation for a 2 : again, after a lengthy but straightforward calculation we derive The evolution equations in the standard first Sonine approximation are easily reached just neglecting nonlinear terms in a 2 in Equations (37), that is, For ξ = 0, the steady solution of these equations is (39) Note that (i) 0 ≤ |a s 2 | ≤ 0.133 for 0 ≤ β ≤ 2, which makes it reasonable to use the first Sonine approximation, and (ii) a s 2 is independent of the driving intensity ξ. This will be useful in the linear response analysis, to be developed below, because a sudden change in the driving only changes the stationary value of the temperature but not that of the excess kurtosis. If ξ = 0, the system evolves towards the homogeneous cooling state, in which the excess kurtosis tends to the value as predicted by Equation (38b), and the temperature decays following Haff's law, dT/dt ∝ −T 1+ β 2 . From now on, we use reduced temperature and time, The steady temperature T s plays the role of a natural energy (or granular temperature) unit. In reduced variables, the evolution equations are where we have introduced two parameters of the order of unity, The evolution equations in the first Sonine approximation, (38) or (42), are the particularisation of the equations for the moments (22) to our model: J = 2, and z 1 = T (or θ), z 2 = a 2 . Consistently, they are nonlinear although here, due to the simplifications introduced in the first Sonine approximation, only nonlinear in θ. When the system is close to the NESS, Equations (42) can be linearised around it by writing θ = 1 + ∆θ, a 2 = a s 2 + ∆a 2 , d ds Of course, the general solution of this linear system for arbitrary initial conditions ∆θ(0) and ∆a 2 (0) can be immediately written, but we omit it here.

Kovacs hump in linear response
Now we look into the Kovacs hump in the linear response approximation. Following the discussion leading to Equation (20), first we have to calculate the relaxation function φ T for the granular temperature. The system is at the steady state corresponding to a driving ξ 0 for t < 0, at t = 0 the driving is instantaneously changed to ξ and only the linear terms in ∆ξ = ξ − ξ 0 are retained. We choose the normalisation of φ T (s) in such a way that φ T (0) = 1, that is Since T s changes with ξ but a 2 does not, we have to solve Equation (44) for ∆a 2 (0) = 0 and arbitrary (small enough) ∆θ(0). The solution is where M ij is the (i, j) element of the matrix M and λ ± its eigenvalues, Both eigenvalues λ ± are negative, since Tr M < 0 and det M > 0 for all β > 0. Therefore, |λ + | < |λ − | and it is λ + that dominates the relaxation of the granular temperature for long times. Moreover, c ± > 0 and thus the linear relaxation function φ T (s) is always positive and decays monotonically to zero. Next we consider a Kovacs-like experiment: the system was at the NESS corresponding to a driving ξ 0 , with granular temperature T s,0 for t < 0, the driving is suddenly changed to ξ 1 at t = 0 so that the system starts to relax towards a new steady temperature T s,1 for 0 ≤ t ≤ t w , and this relaxation is interrupted at t = t w , because the driving is again suddenly changed to the value ξ such that the stationary granular temperature T s equals its instantaneous value at t w . The time evolution of the granular temperature for t ≥ t w is given by the particularisation of Equations (20) and (21) to our situation, that is, where we have made use of the normalisation φ T (0) = 1. In the linear response approximation, the jumps in the driving values can be substituted by the corresponding jumps in the stationary values of the granular temperature.
It is important to stress that the positivity and monotonic decay to zero of the linear relaxation function φ T (s) assures that the Kovacs behaviour is normal, that is, (i) K T (s) is always positive and bounded from above by φ T (s) and (ii) there is only one maximum at a certain time s k > s w [6]. The anomalous behaviour found in the uniformly heated hard-sphere granular for large enough inelasticity is thus not present here. This is consistent with the quasi-elastic limit we have introduced to simplify the collision operator.

Nonlinear Kovacs hump
Here we consider the Kovacs hump for arbitrary large driving jumps. In our model, we can make use of the smallness of a 2 , which is assumed in the first Sonine approximation, in order to introduce a perturbative expansion of Equations (42) in powers of a s 2 . The procedure is completely analogous to that performed in [12,13] for a dilute gas of inelastic hard spheres and thus we omit the details here.
We start by writing a 2 = a s 2 A 2 , with A 2 of the order of unity, and θ(s) = θ 0 (s) + a s 2 θ 1 (s) + . . . , These expansions are inserted into Equations (42), which have to be solved with initial conditions To the lowest order, θ 0 (s) = 1 whereas A 20 (s) decays exponentially to one, In order to describe the Kovacs hump, we compute θ 1 (s) that verifies the differential equation which can be immediately integrated to give The structure of this result is completely analogous to those in [12,13] and thus the conclusions can also be drawn in a similar way. In particular, we want to highlight that (i) the factor that controls the size of the hump is proportional to a ini 2 − a s 2 , and (ii) the shape of the hump is codified in the factor between brackets that only depends on β. Note that (a ini 2 − a s 2 ) > 0 for the cooling protocols (ξ 1 < ξ < ξ 0 ) considered here and thus no anomalous Kovacs hump is expected in the nonlinear regime either.

Numerical results
Here, we compare the theoretical approach above to numerical results of our model. Specifically, we focus on the case β = 1 that gives a collision rate similar to that of hard-spheres. All simulations have been carried out with a restitution coefficient α = 0.999, which corresponds to the quasi-elastic limit in which our simplified kinetic description holds. Also, we set ω = 1 without loss of generality.

Validation of the evolution equations and linear relaxation
First of all, we check the validity of the first Sonine approximation, as given by Equations (38), to describe the time evolution of our system. In order to do so, we compare several relaxation curves between the NESS corresponding to two different noise strengths. In particular, we always depart from the stationary state corresponding to χ 0 = 1 and afterwards let the system evolve with χ = {0.2, 0.6, 0.8, 1} for t > 0. In Figure 2, we compare the Monte-Carlo simulation (see appendix A for details) of the system with the numerical solution of the evolution equations in the first Sonine approximation (38). In addition, we also have plotted the analytical solution of the linear response system, Equation (44). The agreement is complete between simulation and theory, and it can be observed how the linear response result becomes more accurate as the temperature jump decreases.
In order not to clutter the plot in Figure 2, we have not shown the results for the nonlinear first Sonine approximation, Equations (37). The relative error between their numerical solution and that of the standard first Sonine approximation (38) is at most of order 0.1%, for all the cases we have considered. Henceforth, we always use the latter, which is the usual approach in kinetic theory.

Kovacs hump
Since the numerical integration of the first Sonine approximation perfectly agrees with Monte Carlo simulations, we compare the analytical results for the Kovacs hump with the former. Specifically, we work in reduced variables and therefore we integrate numerically Equations (42).

Linear response
It is convenient to rewrite the expression for the Kovacs hump in an alternative form to compare our theory to numerical results. We take advantage of the simple structure of the relaxation function in the first Sonine approximation, which is the sum of two exponentials, to introduce the factorisation [6] K T (s) = K 0 (s w )K 1 (s − s w ), where Firstly, this factorisation property shows that the position s k of the maximum relative to the waiting time s w , that is, s k − s w , is controlled by the function K 1 . Thus, s k − s w does not depend on the waiting time but only on the two eigenvalues λ ± . Namely, Secondly, the height of the maximum K max does depend on the waiting time s w due to the factor K 0 (s w ). Specifically, it can be shown that K max is a monotonically decreasing function of the waiting time s w that vanishes in the limit as s w → ∞. In order to check the above results, we have fixed the initial and final drivings in the Kovacs protocol χ 0 and χ and changed the intermediate driving value χ 1 . We do so to simplify the comparison, because the time scale s involves the steady value of the temperature, see Equation (41). Note that the smaller χ 1 is, the shorter the waiting time becomes. Therefore, one expects to get a Kovacs hump whose maximum remains at s − s w 0.44 but raises as χ 1 decreases. This is shown in Figure 3, where the numerical solution of the first Sonine approximation equations (42) and the analytical result (53) are compared. Their agreement is almost perfect for the two lowest curves, corresponding to χ 1 = 0.99 and χ 1 = 0.95, as expected, but is still remarkably good for the two topmost ones, corresponding to the not-so-small jumps for χ 1 = 0.8 and χ 1 = 0.5.

Nonlinear regime
Also we explore the Kovacs effect out of the linear regime. Figure 4 is similar to Figure 3, but for larger temperature (or driving) jumps. We have also fixed the initial and final values of the driving, χ 0 = 10 and χ = 1. The intermediate values of the driving are the same as in the linear case except for the largest one, χ 1 = 0.99, which we have omitted for the sake of clarity (its hump is too small in the scale of the figure). Now, the linear response theory results just provide the qualitative behaviour of the hump, correctly predicting the position of the maximum but not its height. On the one hand, and consistently with the numerical results in an active matter model [15], the Kovacs hump out the linear response regime is larger than the prediction of linear response theory. On the other hand, the position of the maximum remains basically unchanged and its height still increases as χ 1 decreases.
We also compare our analytical expansion in a s 2 with the numerical solutions of Equations (42) for large jumps. Specifically, in order to make things as simple as possible, we choose χ 1 = 0. If the waiting time is long enough, the system reaches the homogeneous cooling state and a 2 (s w ) = a HCS 2 , which is then the initial condition for the final stage of the Kovacs protocol. Moreover, we can compute the location of the maximum of the hump from Equation (52), obtaining This result is numerically indistinguishable from that of linear response, as given by Equation (54), since the relative error is around 1%. In Figure 5, we put forward a comparison between the Kovacs hump obtained from the numerical solution of the first Sonine approximation equations and our theoretical expression for the nonlinear regime, Equation (52). Fixing ξ 1 = 0 and ξ = 1, as ξ 0 increases (as the waiting time is increased), the hump approaches Equation (52) with a ini 2 = a HCS 2 . Moreover, our theory perfectly reproduces all the numerical curves when we substitute the actual values of a ini 2 into Equation (52) .

Monotonicity of an H-functional
The non-monotonicity in the relaxation of the granular temperature that is brought about by the Kovacs protocol is not automatically transferred to other relevant physical magnitudes. Specifically, here we deal with H-functional   There is strong numerical evidence about H(t) being a Lyapunov functional for granular fluids, thus allowing it to be considered as an out-of-equilibrium entropy relative to that of the NESS, in this context [19,20]. Moreover, it has been analytically proven that H(t) is a Lyapunov functional in our system for the Maxwell collision rule β = 0 [21].
We have computed H(t) numerically from Equation (56) within the first Sonine approximation, 2 for the Kovacs protocols considered in Figures 3 and 4. Once more, we have taken β = 1, for which the analytical proof in [21] does not hold. The results are shown in Figure 6 and make it clear that H(t) still monotonically decreases for the Kovacs-like protocols, in both the linear (left panel) and nonlinear (right panel) regimes. At the time of the maximum in the hump, s − s w 0.44, no special signature is observed in the entropy.

Discussion
One of the main results in our paper is the extension of the linear response expression for the Kovacs hump in thermal systems, as given by Equations (1) and (2), to the realm of athermal systems, Equations (20) and (21). This extension is (i) not trivial, since athermal systems tend in the long time limit to a NESS, not to an equilibrium state and (ii) very general, since it can be done starting from the evolution equations at either the mesoscopic or the macroscopic level of description. Therefore, it means that a Kovacs-like effect is to be expected for a very wide class of physical properties Y in a very wide class of systems, basically as long as the relaxation function φ Y (t) to the NESS is monotonic.
This theoretical result has been checked in a class of systems that mimic the dynamics of the shear component of the velocity of a granular fluid. In the linear response regime, we have found a good agreement between the theoretical prediction and the numerical results. Furthermore, we have investigated how the linear response result extends to the nonlinear regime. In this region, the linear response theory results remain useful at a qualitative level, but the actual values for the hump lie well above the linear prediction. Interestingly, this kind of giant Kovacs hump has been recently reported in active systems [15]. In addition, the nonlinear Kovacs hump can be theoretically explained by an expansion in the excess kurtosis. 2 That is, we have substituted both f (v, t) and f s (v) by their expressions in the first Sonine approximation and calculated the integral numerically.
The work presented here also opens new interesting perspectives for future research. First, our analysis of the Kovacs effect in the model, being restricted to the quasi-elastic limit, has not found the anomalous behaviour shown by a gas of inelastic hard spheres for large enough inelasticity in the nonlinear regime [12,13]. The possibility of such a behaviour in linear response, either in the model or in the granular gas, deserves further investigation. Second, our work clearly shows the compatibility of the non-monotonic decay of the granular temperature (or the corresponding relevant physical variable) and the monotonic decay of the nonequilibrium entropy or H-functional [19][20][21]31]. In this respect, to elucidate if the hump leaves some signature in the decay of the nonequilibrium entropy is compelling.