Kovacs-Like Memory Effect in Athermal Systems: Linear Response Analysis

Abstract

We analyze 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. Furthermore, we prove that the observed non-monotonic relaxation is consistent with the monotonic decay of the non-equilibrium entropy.

1. Introduction

The equilibrium state of physical systems is characterized by the value of a few macroscopic variables, for example pressure, volume and temperature in fluids. This characterization 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 characterized 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 characterize 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 equaled 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 $\varphi \left(t\right)$ of the energy E is recorded. Specifically, $\varphi \left(t\right)=\langle E\left(t\right)\rangle -{\langle E\rangle }_{\mathrm{e}}$, where ${\langle E\rangle }_{\mathrm{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 $\langle E\rangle (t=t_w)$ equals ${\langle E\rangle }_{\mathrm{e}}$. At this time $t_w$, the temperature is increased to T, but the energy does not remain constant: it displays the behavior that is qualitatively shown by $K\left(t\right)$. At first, $K\left(t\right)$ increases from zero until a maximum is attained for $t=t_k$, and only afterwards, it goes back to zero. Similarly to the relaxation function, we have defined $K\left(t\right)=\langle E\left(t\right)\rangle -{\langle E\rangle }_{\mathrm{e}}$, for $t\ge t_w$. Note that $K\left(t\right)\le \varphi \left(t\right)$ 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]

$$K\left(t\right)=\frac{T_0-T_1}{T_0-T}\varphi \left(t\right)-\frac{T-T_1}{T_0-T}\varphi (t-t_w),$$

(1)

where the final temperature T and the waiting time $t_w$ are related by:

$$\frac{T-T_1}{T_0-T_1}=\frac{\varphi \left(t_w\right)}{\varphi \left(0\right)}.$$

(2)

In linear response, the relaxation function $\varphi \left(t\right)$ decays monotonically in time because it is proportional to the equilibrium time correlation function (fluctuation-dissipation theorem) ${\langle E\left(0\right)E\left(t\right)\rangle }_{\mathrm{e}}-{\langle E\rangle }_{\mathrm{e}}^2={\sum }_i{c}_{i}exp\left({\lambda }_{i}t\right)$, with $c_i>0$ and ${\lambda }_i<0$ for all i [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\le K\left(t\right)\le \varphi \left(t\right)$, which assures that the hump always has a positive sign (from now on, “normal” behavior), (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 on 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 $\xi$ 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\to \xi$, $E\to T_g$. One of the main differences found is the emergence of “anomalous” Kovacs behavior for large enough dissipation, when $K\left(t\right)$ 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 ${\xi }_0-\xi$, ${\xi }_0-{\xi }_1$ that are not small.

More recently, Kovacs-like behavior 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 two-fold. 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 behavior 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].

2. Linear Theory for Kovacs-Like Memory Effects

2.1. General Markovian Dynamics

Here, we consider a general system, whose state is completely characterized by a vector $\mathit{x}$ with M components, $\mathit{x}=\{x_1,x_2,\dots ,x_M\}$. For example, in a one-dimensional Ising chain of N spins, $x_i={\sigma }_i=\pm 1$, and $M=N$; for a gas comprising N particles with positions ${\mathit{r}}_i$ and velocities ${\mathit{v}}_i$, $M=6N$, and $\mathit{x}=\{{\mathit{r}}_1,{\mathit{v}}_1,\dots ,{\mathit{r}}_N,{\mathit{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 labeled with a discrete index $\alpha$, $1\le \alpha \le \Omega$. For example, this is the case of the Ising system above, where $\Omega =2^N$. The generalization for a continuous index is straightforward, by changing sums into integrals and the Kronecker delta by the Dirac delta [9].

At the mesoscopic level of description, we assume that $\mathit{x}$ is a Markov process, and its dynamics is governed by a master equation for the probabilities $P({\mathit{x}}_{\alpha },t)$,

$${\partial }_{t}P({\mathit{x}}_{\alpha },t)=\sum _{\tilde{\alpha }}\left[W\left({\mathit{x}}_{\alpha }\right|{\mathit{x}}_{\tilde{\alpha }};\xi )P({\mathit{x}}_{\tilde{\alpha }},t)-W\left({\mathit{x}}_{\tilde{\alpha }}\right|{\mathit{x}}_{\alpha };\xi )P({\mathit{x}}_{\alpha },t)\right],$$

(3)

where $W\left({\mathit{x}}_{\alpha }\right|{\mathit{x}}_{\tilde{\alpha }};\xi )$ are the transition rates from state ${\mathit{x}}_{\tilde{\alpha }}$ to state ${\mathit{x}}_{\alpha }$, and our notation explicitly marks their dependence on some control parameter $\xi$. Equation (3) can be formally written as:

$${\partial }_t\left|\mathcal{P}\right(t)\rangle =\mathbb{W}\left(\xi \right)\left|\mathcal{P}\right(t)\rangle ,$$

(4)

where $\left|\mathcal{P}\right(t)\rangle$ is a vector (column matrix) whose components are the probabilities $P({\mathit{x}}_{\alpha },t)$ and $\mathbb{W}\left(\xi \right)$ is the linear operator (square matrix) that generates the dynamical evolution of $\left|\mathcal{P}\right(t)\rangle$,

$$\mathbb{W}\left({\mathit{x}}_{\alpha }\right|{\mathit{x}}_{\tilde{\alpha }};\xi )=W\left({\mathit{x}}_{\alpha }\right|{\mathit{x}}_{\tilde{\alpha }};\xi )-{\delta }_{\alpha ,\tilde{\alpha }}\sum _{{\alpha }^{\prime }}W\left({\mathit{x}}_{{\alpha }^{\prime }}\right|{\mathit{x}}_{\alpha };\xi ).$$

(5)

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 $|{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle$, which verifies:

$$\mathbb{W}\left(\xi \right)|{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle =0,$$

(6)

and depends on the parameter $\xi$ 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 an NESS in the long time limit, not an equilibrium state.

Now, we consider the system evolving from a certain initial state at time $t_0$, characterized by the distribution $\left|\mathcal{P}\right(t_0)\rangle$. The formal solution of the master equation can be written as:

$$|\mathcal{P}\left(t\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle =e^{(t-t_0)\mathbb{W}\left(\xi \right)}|\mathcal{P}\left(t_0\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle .$$

(7)

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 $\mathit{x}$ by $\mathit{Y}\left(\mathit{x}\right)$: its expected or average value is given by :

$$\langle \mathit{Y}\left(t\right)\rangle =\sum _{\alpha }\mathit{Y}\left({\mathit{x}}_{\alpha }\right)P({\mathit{x}}_{\alpha },t)=\langle \mathcal{Y}\left|\mathcal{P}\right(t)\rangle ,$$

(8)

where $|\mathcal{Y}\rangle$ is a ket whose components are $\mathit{Y}\left({\mathit{x}}_{\alpha }\right)$, and $\langle \mathcal{Y}|$ its corresponding bra (row matrix with the same components. Note that we are assuming that Y is a real quantity for all the configurations. By substituting (7) into (8), it is obtained that:

$$\Delta \mathit{Y}(t;\xi )\equiv \langle \mathit{Y}\left(t\right)\rangle -{\langle \mathit{Y}\rangle }_{\mathrm{s}}\left(\xi \right)=\langle \mathcal{Y}|e^{(t-t_0)\mathbb{W}\left(\xi \right)}|\mathcal{P}\left(t_0\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle ,$$

(9)

in which ${\langle \mathit{Y}\rangle }_{\mathrm{s}}\left(\xi \right)$ is the average value at the steady state of the system corresponding to $\xi$.

Now, we investigate the relaxation of the system from the steady state for ${\xi }_0=\xi +\Delta \xi$ to the steady state for $\xi$. We do so in linear response, that is $\Delta \xi$ is considered to be small, and we neglect all terms beyond those linear in $\Delta \xi$. Thus, at $t=0$ we have that the probability distribution is:

$$\left|\mathcal{P}\right(t=0)\rangle =|P_{\mathrm{s}}(\xi +\Delta \xi )\rangle =|P_{\mathrm{s}}\left(\xi \right)\rangle +\Delta \xi |\frac{d{P}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle +O{(\Delta \xi )}^2.$$

(10)

Substitution of (10) into (9) yields the formal expression for the relaxation of Y in linear response,

$$\Delta \mathit{Y}(t;\xi )=\Delta \xi \langle \mathcal{Y}\left|e^{(t-t_0)\mathbb{W}\left(\xi \right)}\right|\frac{d{P}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle .$$

(11)

In order to have an order of unity function, one may define a normalized relaxation function:

$${\varphi }_{\mathit{Y}}(t;\xi )\equiv \underset{\Delta \xi \to 0}{lim}\frac{\Delta \mathit{Y}(t;\xi )}{\Delta \xi }=\langle \mathcal{Y}\left|e^{(t-t_0)\mathbb{W}\left(\xi \right)}\right|\frac{d{P}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle .$$

(12)

Sometimes, the relaxation function is further normalized by considering ${\varphi }_{\mathit{Y}}\left(t\right)/{\varphi }_{\mathit{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.

2.2. 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 $\xi$ controlling the dynamics is changed in the following stepwise manner:

$$\xi \left(t\right)=\left\{\begin{array}{cc}{\xi }_0,\hfill & -\infty <t<0,\hfill \\ {\xi }_1,\hfill & 0<t<t_w,\hfill \\ \xi ,\hfill & t>t_w.\hfill \end{array}\right.$$

(13)

Therefore, since ${\xi }_0$ is kept for an infinite time, at $t=0$, the system is prepared in the corresponding steady state, $\mathcal{P}(t=0)={\mathcal{P}}_{\mathrm{s}}\left({\xi }_0\right)$. Our idea is to consider that the jumps ${\xi }_1-{\xi }_0$ and $\xi -{\xi }_1$ are small, in the sense that all expressions can be linearized 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 $\xi$ playing the role of the temperature.

We start by analyzing the relaxation in the first time window, $0<t<t_w$. Therein, we apply (7) with the substitutions $t_0\to 0$ and $\xi \to {\xi }_1$, that is,

$$|\mathcal{P}\left(t\right)-{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)\rangle =e^{t\mathbb{W}\left({\xi }_1\right)}|{\mathcal{P}}_{\mathrm{s}}\left({\xi }_0\right)-{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)\rangle ,0\le t\le t_w.$$

(14)

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 $\xi$. Making use again of (7) with $t_0\to t_w$,

$$\begin{array}{ccc}\hfill |\mathcal{P}\left(t\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle & =& e^{(t-t_w)\mathbb{W}\left(\xi \right)}|\mathcal{P}\left(t_w\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle \hfill \\ & =& e^{(t-t_w)\mathbb{W}\left(\xi \right)}\left[e^{t_w\mathbb{W}\left({\xi }_1\right)}|{\mathcal{P}}_{\mathrm{s}}\left({\xi }_0\right)-{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)\rangle +|{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle \right],t\ge t_w.\hfill \end{array}$$

(15)

It must be stressed that the above expressions, Equations (14) and (15), are exact, and no approximation has been made.

The linear response approximation is introduced now: we assume that both jumps ${\xi }_0-{\xi }_1$ and $\xi -{\xi }_1$ are small, so we can expand both $|{\mathcal{P}}_{\mathrm{s}}\left({\xi }_0\right)-{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)\rangle$ and $|{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle$, similarly to what was done in Equation (10). Namely,

$$|{\mathcal{P}}_{\mathrm{s}}\left({\xi }_0\right)-{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)\rangle =({\xi }_0-{\xi }_1)|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle +O{({\xi }_0-{\xi }_1)}^2,$$

(16a)

$$|{\mathcal{P}}_{\mathrm{s}}\left({\xi }_1\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle =({\xi }_1-\xi )|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle +O{({\xi }_1-\xi )}^2.$$

(16b)

Note that in both (16a) and (16b), the derivatives are evaluated at $\xi$; the difference introduced by evaluating them at either ${\xi }_1$ or ${\xi }_0$ is second order in the deviations. Then,

$$|\mathcal{P}\left(t\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle =({\xi }_0-{\xi }_1)e^{(t-t_w)\mathbb{W}\left(\xi \right)}{e}^{t_w\mathbb{W}\left({\xi }_1\right)}|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle +({\xi }_1-\xi )e^{(t-t_w)\mathbb{W}\left(\xi \right)}|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle .$$

(17)

This equation can be further simplified: since the two terms on its rhs are first order in the jumps, we can substitute $\mathbb{W}\left({\xi }_1\right)$ with $\mathbb{W}\left(\xi \right)$, with the result:

$$|\mathcal{P}\left(t\right)-{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)\rangle =({\xi }_0-{\xi }_1)e^{t\mathbb{W}\left(\xi \right)}|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle -(\xi -{\xi }_1)e^{(t-t_w)\mathbb{W}\left(\xi \right)}|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle .$$

(18)

This is the superposition of two responses: the first term on the rhs gives the relaxation from ${\xi }_0$ to ${\xi }_1$, which starts at $t=0$, whereas the second term stands for the relaxation from ${\xi }_1$ to $\xi$, which starts at $t=t_w$. We have chosen to write $-(\xi -{\xi }_1)$ in the second term because $\xi >{\xi }_1$ in the Kovacs protocol.

The same structure in Equation (18) is transferred to the average value. Taking into account Equation (9):

$$\Delta \mathit{Y}\left(t\right)=({\xi }_0-{\xi }_1)\langle \mathcal{Y}\left|e^{t\mathbb{W}\left(\xi \right)}\right|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle -(\xi -{\xi }_1)\langle \mathcal{Y}\left|e^{(t-t_w)\mathbb{W}\left(\xi \right)}\right|\frac{d{\mathcal{P}}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle ,t\ge t_w,$$

(19)

in which we recognize the relaxation function in linear response, defined in Equation (12). We can also normalize the response in this experiment, by defining a function $K\left(t\right)$ as follows,

$$K_{\mathit{Y}}\left(t\right)\equiv \underset{{\xi }_0\to \xi }{lim}\frac{\Delta \mathit{Y}\left(t\right)}{{\xi }_0-\xi }=\frac{{\xi }_0-{\xi }_1}{{\xi }_0-\xi }{\varphi }_{\mathit{Y}}\left(t\right)-\frac{\xi -{\xi }_1}{{\xi }_0-\xi }{\varphi }_{\mathit{Y}}(t-t_w).$$

(20)

It is understood that, as ${\xi }_0-\xi \to 0$, both prefactors $\frac{{\xi }_0-{\xi }_1}{{\xi }_0-\xi }$ and $\frac{\xi -{\xi }_1}{{\xi }_0-\xi }$ 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 $({\xi }_0,{\xi }_1,\xi )$, provided that the jumps are small enough and the ratio of the jumps is of the order of unity. The function $K\left(t\right)$ corresponds to a Kovacs-like experiment when $\xi$ is chosen as a function of $t_w$ in such a way that $\langle \mathit{Y}\left(t_w\right)\rangle ={\langle \mathit{Y}\rangle }_{\mathrm{s}}\left(\xi \right)$ or $K_{\mathit{Y}}\left(t_w\right)=0$, that is,

$$\frac{\xi -{\xi }_1}{{\xi }_0-{\xi }_1}=\frac{{\varphi }_{\mathit{Y}}\left(t_w\right)}{{\varphi }_{\mathit{Y}}\left(0\right)}.$$

(21)

Alternatively, one may consider that (21) defines $t_w$ as a function of $\xi$.

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. Furthermore, it must be noted that it can be easily extended to the Fokker–Planck, or the equivalent Langevin, equation.

2.3. 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,\dots ,J$, where J is the number of relevant moments. The equations for the moments have the general form:

$$\frac{d}{dt}{z}_i=f_i(z_1,\dots ,z_J;\xi ),$$

(22)

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. This 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_i^{\mathrm{s}}\left(\xi \right)$. Linearization of the dynamical system around the steady state gives:

$$\frac{d}{dt}|\Delta z\left(t\right)\rangle =\mathbb{M}\left(\xi \right)|\Delta z\left(t\right)\rangle ,|\Delta z\left(t\right)\rangle \equiv |z\left(t\right)-z_{\mathrm{s}}\left(\xi \right)\rangle .$$

(23)

We are using a notation completely similar to that in the previous section: $|z\rangle$ is a vector, represented by a column matrix with components $z_i$, and $\mathbb{M}\left(\xi \right)$ is a linear operator, represented by a square matrix with elements:

$$M_{ij}\left(\xi \right)={\left.{\partial }_{z_j}{f}_i\right|}_{|z\rangle =|z_{\mathrm{s}}\left(\xi \right)\rangle }.$$

(24)

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}\ne M_{ji}$, and the operator $\mathbb{M}$ is not Hermitian. However, we do not need $\mathbb{M}$ to be Hermitian to solve the linearized system in a formal way, as shown below.

Analogously to what was done for the master equation, the formal solution of Equation (23) is:

$$|\Delta z\left(t\right)\rangle =e^{(t-t_0)\mathbb{M}\left(\xi \right)}|\Delta z\left(t_0\right)\rangle .$$

(25)

In particular, if the initial condition is chosen to correspond to the steady state for ${\xi }_0=\xi +\Delta \xi$, one has:

$$|\Delta z\left(t\right)\rangle =\Delta \xi e^{(t-t_0)\mathbb{M}\left(\xi \right)}|\frac{d{z}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle .$$

(26)

The response for any of the relevant moments can be extracted by projecting the above result onto the “natural” basis $|u_i\rangle$, whose components are $u_{ij}={\delta }_{ij}$. Then, the normalized linear response function for $z_i$ can be defined by:

$${\varphi }_{z_i}\left(t\right)=\underset{\Delta \xi \to 0}{lim}\frac{\langle u_i|\Delta z\left(t\right)\rangle }{\Delta \xi }=\langle u_i\left|e^{(t-t_0)\mathbb{M}\left(\xi \right)}\right|\frac{d{z}_{\mathrm{s}}\left(\xi \right)}{d\xi }\rangle .$$

(27)

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); thus, it is not repeated here.

3. A Lattice Model with Conserved Momentum and Non-Conserved Energy

3.1. 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 $\mathit{v}\equiv \{v_1,...,v_N\}$. The dynamics proceeds through inelastic nearest-neighbor binary collisions: each pair $(l,l+1)$ collides inelastically with a characteristic rate proportional to $|v_l-v_{l+1}{|}^{\beta }$, with $\beta \ge 0$. For $\beta =0$, nearest-neighbor particles collide independently of their relative velocity (the so-called Maxwell-molecule model [24]), whereas for $\beta =1$ and $\beta =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 operator ${\widehat{b}}_l$,

$$\begin{array}{ccc}\hfill {\widehat{b}}_l{v}_l& =& v_l-\frac{1+\alpha }{2}\left(v_l-v_{l+1}\right),\hfill \end{array}$$

(28a)

$$\begin{array}{ccc}\hfill {\widehat{b}}_l{v}_{l+1}& =& v_{l+1}+\frac{1+\alpha }{2}\left(v_l-v_{l+1}\right),\hfill \end{array}$$

(28b)

where $0<\alpha \le 1$ is the normal restitution coefficient. Momentum is always conserved in collisions, $({\widehat{b}}_l-1)(v_l+v_{l+1})=0$, but energy is not; in fact, $({\widehat{b}}_l-1)(v_l^2+v_{l+1}^2)=({\alpha }^2-1){(v_l-v_{l+1})}^2/2\le 0$, with equality holding only in the elastic limit $\alpha =1$.

In order to have a steady state, a mechanism that injects energy into the system, and thus compensates on 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:

$${\left.{\partial }_{\tau }{v}_i\left(\tau \right)\right|}_{\mathrm{noise}}={\xi }_i\left(\tau \right)-\frac{1}{N}\sum _{j=1}^N{\xi }_j\left(\tau \right),$$

(29)

and ${\xi }_i\left(\tau \right)$ are Gaussian white noises of amplitude $\chi$:

$${\langle {\xi }_i\left(\tau \right)\rangle }_{\mathrm{noise}}=0,{\langle {\xi }_i\left(\tau \right){\xi }_j\left({\tau }^{\prime }\right)\rangle }_{\mathrm{noise}}=\chi {\delta }_{ij}\delta (\tau -{\tau }^{\prime }),$$

(30)

for $i,j=1,\dots ,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(\mathit{v},\tau )$ probability density of finding the system in state $\mathit{v}$ at time $\tau$. Instead, we directly write the “kinetic” equation for the one-particle distribution function, namely $P_1(v;l,\tau )=\int d\mathit{v}{P}_N(\mathit{v},\tau )\delta (v_l-v)$. Therefrom, all the one-site velocity moments, which are the relevant physical quantities for our present purposes, can be derived, $\langle v_l^n\left(\tau \right)\rangle \equiv {\int }_{-\infty }^{+\infty }dv{v}^n{P}_1(v;l,\tau )$. 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 analyze 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,\tau )\equiv P_1(v,l;\tau )$ and consider the quasi-elastic limit, that is $ϵ\equiv 1-{\alpha }^2\ll 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]:

$${\partial }_{\tau }f(v,\tau )=\frac{\omega ϵ}{2}{\partial }_v{\int }_{-\infty }^{+\infty }d{v}^{\prime }(v-v^{\prime }){|v-v^{\prime }|}^{\beta }f(v,\tau )f(v^{\prime },\tau )+\frac{\chi }{2}{\partial }_v^{2}f(v,\tau ).$$

(31)

We can define a dimensionless time scale, $t=\omega ϵ\tau$, over which:

$${\partial }_{t}f(v,t)=\frac{1}{2}{\partial }_v{\int }_{-\infty }^{+\infty }d{v}^{\prime }(v-v^{\prime }){|v-v^{\prime }|}^{\beta }f(v,t)f(v^{\prime },t)+\frac{\xi }{2}{\partial }_v^{2}f(v,t),$$

(32)

where $\xi$ is the rescaled strength of the noise, $\xi =\frac{\chi }{\omega ϵ}$. 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:

$$T\equiv \langle v^2\rangle ={\int }_{-\infty }^{+\infty }dv{v}^{2}f(v,t).$$

(33)

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 $v_0=\sqrt{2T}$ to make v dimensionless, with the change:

$$v=v_{0}c,f(v,t)dv=\phi (c,t)dc\iff \phi (c,t)=v_{0}f(v,t).$$

(34)

Taking moments in (32) and making the change of variables above, one gets:

$$\frac{d}{dt}T=-\zeta T^{1+\frac{\beta }{2}}+\xi ,\zeta =2^{\frac{\beta }{2}}{\int }_{-\infty }^{+\infty }dc{\int }_{-\infty }^{+\infty }d{c}^{\prime }{|c-c^{\prime }|}^{2+\beta }\phi (c,t)\phi (c^{\prime },t)$$

(35)

The first term on the rhs stems from collisions and cools the system, in the sense that it always makes the granular temperature decrease. The second term stems from the stochastic thermostat and heats the system, and thus, in the long time limit, an NESS is attained in which both terms counterbalance each other.

3.2. 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,

$$\phi (c,t)=\frac{e^{-c^2}}{\sqrt{\mathsf{\pi }}}\left[1+\sum _{k=2}^{\infty }{a}_k\left(t\right)L_k^{\left(-\frac{1}{2}\right)}\left(c^2\right)\right],$$

(36)

where $L_k^{\left(m\right)}\left(x\right)$ are the associated Laguerre polynomials [30]. In kinetic theory, $m=\frac{d}{2}-1$, with d being the spatial dimension, and often, the notation $S_k\left(x\right)\equiv L_k^{\left(\frac{d}{2}-1\right)}$ 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, $\langle c^4\rangle =3(1+a_2)/4$.

Although the linearization 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 the 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, the evolution equation for the temperature is readily obtained [29]:

$$\frac{d}{dt}T=-{\zeta }_0{T}^{1+\frac{\beta }{2}}\left[1+\frac{\beta (2+\beta )}{16}{a}_2+\frac{\beta (2+\beta )(2-\beta )(4-\beta )}{1024}{a}_2^2\right]+\xi ,$$

(37a)

where ${\zeta }_0={\pi }^{-1/2}{2}^{1+\beta }\mathrm{\Gamma }\left(\frac{3+\beta }{2}\right)$. Unless $\beta =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:

$$\begin{array}{cc}\hfill T\frac{d}{dt}{a}_2=& -2\xi a_2-\frac{{\zeta }_0}{3}\beta T^{1+\frac{\beta }{2}}\hfill \\ & \times \left[1+\frac{56+\beta (6+\beta )}{16}{a}_2-\frac{(2+\beta )[384+(2-\beta \left)\beta \right(4+\beta \left)\right]}{1024}{a}_2^2-\frac{3(4-\beta )(2-\beta )(2+\beta )}{512}{a}_2^3\right].\hfill \end{array}$$

(37b)

The evolution equations in the standard first Sonine approximation are easily reached just neglecting nonlinear terms in $a_2$ in Equations (37), that is:

$$\frac{d}{dt}T=-{\zeta }_0{T}^{1+\frac{\beta }{2}}\left[1+\frac{\beta (2+\beta )}{16}{a}_2\right]+\xi ,$$

(38a)

$$T\frac{d}{dt}{a}_2=-\frac{{\zeta }_0}{3}\beta T^{1+\frac{\beta }{2}}\left[1+\frac{56+\beta (6+\beta )}{16}{a}_2\right]-2\xi a_2.$$

(38b)

For $\xi \ne 0$, the steady solution of these equations is:

$$T_s={\left(\frac{\xi }{{\zeta }_0\left[1+\frac{\beta (2+\beta )}{16}{a}_2^{\mathrm{s}}\right]}\right)}^{\frac{2}{2+\beta }},a_2^{\mathrm{s}}=-\frac{16\beta }{96+56\beta +6{\beta }^2+{\beta }^3}.$$

(39)

Note that (i) $0\le |a_2^{\mathrm{s}}|\le 0.133$ for $0\le \beta \le 2$, which makes it reasonable to use the first Sonine approximation, and (ii) $a_2^{\mathrm{s}}$ is independent of the driving intensity $\xi$. 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 $\xi =0$, the system evolves towards the homogeneous cooling state, in which the excess kurtosis tends to the value:

$$a_2^{\mathrm{HCS}}=-\frac{16}{56+\beta (6+\beta )},$$

(40)

as predicted by Equation (38b) and the temperature decays following Haff’s law, $dT/dt\propto -T^{1+\frac{\beta }{2}}$.

From now on, we use reduced temperature and time,

$$\theta =\frac{T}{T_{\mathrm{s}}},s={\zeta }_0{T}_s^{\beta /2}t.$$

(41)

The steady temperature $T_{\mathrm{s}}$ plays the role of a natural energy (or granular temperature) unit. In reduced variables, the evolution equations are:

$$\begin{array}{ccc}\hfill \frac{d}{ds}\theta & =& 1-{\theta }^{1+\frac{\beta }{2}}+\frac{\beta (2+\beta )}{16}\left(a_2^{\mathrm{s}}-a_2{\theta }^{1+\frac{\beta }{2}}\right),\hfill \end{array}$$

(42a)

$$\begin{array}{ccc}\hfill \theta \frac{d}{ds}{a}_2& =& {\kappa }_1\left(a_2-a_2^{\mathrm{HCS}}\right)\left(1-{\theta }^{1+\frac{\beta }{2}}\right)-{\kappa }_2\left(a_2-a_2^{\mathrm{s}}\right),\hfill \end{array}$$

(42b)

where we have introduced two parameters of the order of unity,

$${\kappa }_1=-\frac{\beta }{3a_2^{\mathrm{HCS}}},{\kappa }_2=-\frac{\beta }{3a_2^{\mathrm{s}}},$$

(43)

$0\le {\kappa }_1\le 3$ and $2\le {\kappa }_2\le 5$ for $0\le \beta \le 2$.

The evolution equations in the first Sonine approximation, (38) or (42), are the particularization of the equations for the moments (22) to our model: $J=2$, and $z_1=T$ (or $\theta$), $z_2=a_2$. Consistently, they are nonlinear, although here, due to the simplifications introduced in the first Sonine approximation, only nonlinear in $\theta$. When the system is close to the NESS, Equation (42) can be linearized around it by writing $\theta =1+\Delta \theta$, $a_2=a_2^{\mathrm{s}}+\Delta a_2$,

$$\frac{d}{ds}\left(\begin{array}{c}\Delta \theta \\ \Delta a_2\end{array}\right)=\mathit{M}·\left(\begin{array}{c}\Delta \theta \\ \Delta a_2\end{array}\right),\mathit{M}=\left(\begin{array}{cc}-\frac{2(2+\beta )(12+\beta )}{48+4\beta +{\beta }^2}& -\frac{\beta (2+\beta )}{16}\\ -{\kappa }_1\left(1+\frac{\beta }{2}\right)\left(a_2^{\mathrm{s}}-a_2^{\mathrm{HCS}}\right)& -{\kappa }_2\end{array}\right).$$

(44)

Of course, the general solution of this linear system for arbitrary initial conditions $\Delta \theta \left(0\right)$ and $\Delta a_2\left(0\right)$ can be immediately written, but we omit it here.

3.3. 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 ${\varphi }_T$ for the granular temperature. The system is at the steady state corresponding to a driving ${\xi }_0$ for $t<0$; at $t=0$, the driving is instantaneously changed to $\xi$, and only the linear terms in $\Delta \xi =\xi -{\xi }_0$ are retained. We choose the normalization of ${\varphi }_T\left(s\right)$ in such a way that ${\varphi }_T\left(0\right)=1$, that is:

$${\varphi }_T\left(s\right)\equiv \underset{\Delta T\left(0\right)\to 0}{lim}\frac{\Delta T\left(s\right)}{\Delta T\left(0\right)}=\underset{\Delta \theta \left(0\right)\to 0}{lim}\frac{\Delta \theta \left(s\right)}{\Delta \theta \left(0\right)}.$$

(45)

Since $T_{\mathrm{s}}$ changes with $\xi$, but $a_2$ does not, we have to solve Equation (44) for $\Delta a_2\left(0\right)=0$ and arbitrary (small enough) $\Delta \theta \left(0\right)$. The solution is:

$${\varphi }_T\left(s\right)=c_{+}{e}^{{\lambda }_{+}s}+c_{-}{e}^{{\lambda }_{-}s},$$

(46a)

$$c_{+}=\frac{M_{11}-{\lambda }_{-}}{{\lambda }_{+}-{\lambda }_{-}},c_{-}=\frac{{\lambda }_{+}-M_{11}}{{\lambda }_{+}-{\lambda }_{-}},$$

(46b)

where $M_{ij}$ is the $(i,j)$ element of the matrix $\mathit{M}$ and ${\lambda }_{\pm }$ its eigenvalues,

$${\lambda }_{\pm }=\frac{\mathrm{T}\mathrm{r}\mathit{M}\pm \sqrt{{\left(\mathrm{T}\mathrm{r}\mathit{M}\right)}^2-4det\mathit{M}}}{2}=\frac{\mathrm{T}\mathrm{r}\mathit{M}\pm \sqrt{{(M_{11}-M_{22})}^2+4M_{12}{M}_{21}}}{2}.$$

(47)

Both eigenvalues ${\lambda }_{\pm }$ are negative, since $\mathrm{T}\mathrm{r}\mathit{M}<0$ and $det\mathit{M}>0$ for all $\beta >0$. Therefore, $|{\lambda }_{+}|<|{\lambda }_{-}|$, and it is ${\lambda }_{+}$ that dominates the relaxation of the granular temperature for long times. Moreover, $c_{\pm }>0$, and thus, the linear relaxation function ${\varphi }_T\left(s\right)$ 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 ${\xi }_0$, with granular temperature $T_{\mathrm{s},0}$ for $t<0$; the driving is suddenly changed to ${\xi }_1$ at $t=0$ so that the system starts to relax towards a new steady temperature $T_{\mathrm{s},1}$ for $0\le t\le t_w$, and this relaxation is interrupted at $t=t_w$, because the driving is again suddenly changed to the value $\xi$ such that the stationary granular temperature $T_{\mathrm{s}}$ equals its instantaneous value at $t_w$. The time evolution of the granular temperature for $t\ge t_w$ is given by the particularization of Equations (20) and (21) to our situation, that is,

$$K_T\left(s\right)=\frac{{\xi }_0-{\xi }_1}{{\xi }_0-\xi }{\varphi }_T\left(s\right)-\frac{\xi -{\xi }_1}{{\xi }_0-\xi }{\varphi }_T(s-s_w),\frac{\xi -{\xi }_1}{{\xi }_0-{\xi }_1}={\varphi }_T\left(s_w\right),$$

(48)

where we have made use of the normalization ${\varphi }_T\left(0\right)=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.

The structure of the linear relaxation function ${\varphi }_T\left(s\right)$, as a linear combination of decreasing exponentials $exp\left({\lambda }_{\pm }t\right)$, ${\lambda }_{\pm }<0$, with positive weights $c_{\pm }$, assures that the Kovacs behavior is normal: (i) $K_T\left(s\right)$ is always positive and bounded from above by ${\varphi }_T\left(s\right)$ and (ii) there is only one maximum at a certain time $s_k>s_w$ [6]. The anomalous behavior 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.

3.4. 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_2^{\mathrm{s}}$. 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_2^{\mathrm{s}}{A}_2$, with $A_2$ of the order of unity, and:

$$\theta \left(s\right)={\theta }_0\left(s\right)+a_2^{\mathrm{s}}{\theta }_1\left(s\right)+\dots ,A_2\left(s\right)=A_{20}\left(s\right)+a_2^{\mathrm{s}}{A}_{21}\left(s\right)+\dots .$$

(49)

These expansions are inserted into Equations (42), which have to be solved with initial conditions $\theta \left(s_w\right)=1$, $A_2\left(s_w\right)=A_2^{\mathrm{ini}}$. To the lowest order, ${\theta }_0\left(s\right)=1$, whereas $A_{20}\left(s\right)$ decays exponentially to one,

$$A_{20}\left(s\right)-1\sim \left(A_2^{\mathrm{ini}}-1\right)e^{-{\kappa }_2\left(s-s_w\right)}.$$

(50)

In order to describe the Kovacs hump, we compute ${\theta }_1\left(s\right)$ that verifies the differential equation:

$$\frac{d{\theta }_1}{ds}=-\left(1+\frac{\beta }{2}\right){\theta }_1+\frac{\beta (2+\beta )}{16}\left(A_{20}-1\right),$$

(51)

which can be immediately integrated to give:

$$\theta \left(s\right)-1\sim \left(a_2^{\mathrm{ini}}-a_2^{\mathrm{s}}\right)\frac{\beta (2+\beta )}{8(2+\beta -2{\kappa }_2)}\left[e^{-{\kappa }_2\left(s-s_w\right)}-e^{-\left(1+\frac{\beta }{2}\right)\left(s-s_w\right)}\right],s\ge s_w,$$

(52)

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_2^{\mathrm{ini}}-a_2^{\mathrm{s}}$ and (ii) the shape of the hump is codified in the factor between brackets that only depends on $\beta$. Note that $(a_2^{\mathrm{ini}}-a_2^{\mathrm{s}})>0$ for the cooling protocols (${\xi }_1<\xi <{\xi }_0$) considered here, and thus, no anomalous Kovacs hump is expected in the nonlinear regime either.

4. Numerical Results

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

4.1. 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 ${\chi }_0=1$ and afterwards let the system evolve with $\chi =\{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, Equation (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.

4.2. 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).

4.2.1. 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 factorization [6]:

$$K_T\left(s\right)=K_0\left(s_w\right)K_1(s-s_w),$$

(53a)

where:

$$K_0\left(s_w\right)=c_{+}{c}_{-}\frac{e^{{\lambda }_{+}{s}_w}-e^{{\lambda }_{-}{s}_w}}{1-{\varphi }_T\left(s_w\right)},K_1(s-s_w)=e^{{\lambda }_{+}(s-s_w)}-e^{{\lambda }_{-}(s-s_w)}.$$

(53b)

Firstly, this factorization 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 ${\lambda }_{\pm }$. Namely,

$$s_k-s_w=\frac{1}{{\lambda }_{+}-{\lambda }_{-}}ln\left(\frac{{\lambda }_{-}}{{\lambda }_{+}}\right)\underset{\beta =1}{\simeq }0.442.$$

(54)

Secondly, the height of the maximum $K_{max}$ does depend on the waiting time $s_w$ due to the factor $K_0\left(s_w\right)$. 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\to \infty$.

In order to check the above results, we have fixed the initial and final drivings in the Kovacs protocol ${\chi }_0$ and $\chi$ and changed the intermediate driving value ${\chi }_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 ${\chi }_1$ is, the shorter the waiting time becomes. Therefore, one expects to get a Kovacs hump whose maximum remains at $s-s_w\simeq 0.44$, but raises as ${\chi }_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 ${\chi }_1=0.99$ and ${\chi }_1=0.95$, as expected, but is still remarkably good for the two topmost ones, corresponding to the not-so-small jumps for ${\chi }_1=0.8$ and ${\chi }_1=0.5$.

4.2.2. Nonlinear Regime

Furthermore, 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, ${\chi }_0=10$ and $\chi =1$. The intermediate values of the driving are the same as in the linear case except for the largest one, ${\chi }_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 behavior of the hump, correctly predicting the position of the maximum, but not its height. On the one hand, and consistent with the numerical results in an active matter model [15], the Kovacs hump out of 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 ${\chi }_1$ decreases.

We also compare our analytical expansion in $a_2^{\mathrm{s}}$ with the numerical solutions of Equations (42) for large jumps. Specifically, in order to make things as simple as possible, we choose ${\chi }_1=0$. If the waiting time is long enough, the system reaches the homogeneous cooling state and $a_2\left(s_w\right)=a_2^{\mathrm{HCS}}$, 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:

$$s_k-s_w=\frac{2}{{\kappa }_2-2-\beta }ln\left(\frac{2{\kappa }_2}{2+\beta }\right)\underset{\beta =1}{\simeq }0.437.$$

(55)

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 ${\xi }_1=0$ and $\xi =1$, as ${\xi }_0$ increases (as the waiting time is increased), the hump approaches Equation (52) with $a_2^{\mathrm{ini}}=a_2^{\mathrm{HCS}}$. Moreover, our theory perfectly reproduces all the numerical curves when we substitute the actual values of $a_2^{\mathrm{ini}}$ into Equation (52).

4.3. 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 the H-functional:

$$H\left(t\right)={\int }_{-\infty }^{+\infty }dvf(v,t)log\left[\frac{f(v,t)}{f_{\mathrm{s}}\left(v\right)}\right],$$

(56)

There is strong numerical evidence about $H\left(t\right)$ 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\left(t\right)$ is a Lyapunov functional in our system for the Maxwell collision rule $\beta =0$ [21].

We have computed $H\left(t\right)$ numerically from Equation (56) within the first Sonine approximation (that is, we have substituted both $f(v,t)$ and $f_{\mathrm{s}}\left(v\right)$ by their expressions in the first Sonine approximation and calculated the integral numerically) for the Kovacs protocols considered in Figure 3 and Figure 4. Once more, we have taken $\beta =1$, for which the analytical proof in [21] does not hold. The results are shown in Figure 6 and make it clear that $H\left(t\right)$ 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\simeq 0.44$, no special signature is observed in the entropy.

5. 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 an 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 ${\varphi }_{\mathit{Y}}\left(t\right)$ 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.

Our work also opens new interesting perspectives for future research. Remarkably, the Kovacs-like memory effects analyzed here involve only jump protocols for the control parameter, which is time-independent both in the waiting time window $0<t<t_w$ (value ${\xi }_1$) and in the return to stationarity region $t>t_w$ (value $\xi$). Very recently, the response of the system to continuous perturbations [31,32] has been studied, which raises the interesting question of building up a theory, similar to ours, in which the control is continuously varied up to $t_w$ in a first stage and then abruptly changed to the value that corresponds to the stationary value of the relaxing physical observable. To the best of our knowledge, a thorough study of this kind of response is lacking. Moreover, the analysis of the nonlinear response in light of generalized fluctuation-dissipation relations [33] is worthy of consideration.

On another note, our analysis of the Kovacs effect in the model, being restricted to the quasi-elastic limit, has not found the anomalous behavior shown by a gas of inelastic hard spheres for large enough inelasticity in the nonlinear regime [12,13]. The possibility of such a behavior in linear response, either in the model or in the granular gas, deserves further investigation. In addition, 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,34]. In this respect, to elucidate if the hump leaves some signature in the decay of the nonequilibrium entropy is compelling.