On the emergence of large and complex memory effects in nonequilibrium fluids

Control of cooling and heating processes is essential in many industrial and biological processes. In fact, the time evolution of an observable quantity may differ according to the previous history of the system. For example, a system that is being subject to cooling and then, at a given time $t_{w}$ for which the instantaneous temperature is $T(t_w)=T_{\mathrm{st}}$, is suddenly put in contact with a temperature source at $T_{\mathrm{st}}$ may continue cooling down temporarily or, on the contrary, undergo a temperature rebound. According to current knowledge, there can be only one"spurious"and small peak/low. However, our results prove that, under certain conditions, more than one extremum may appear. Specifically, we have observed regions with two extrema and a critical point with three extrema. We have also detected cases where extraordinarily large extrema are observed, as large as the order of magnitude of the stationary value of the variable of interest. We show this by studying the thermal evolution of a low density set of macroscopic particles that do not preserve kinetic energy upon collision, i.e., a granular gas. We describe the mechanism that signals in this system the emergence of these complex and large memory effects, and explain why similar observations can be expected in a variety of systems.

Control of cooling and heating processes is an essential aspect in industrial 1 and biological processes 2 . In fact, the time evolution of an observable quantity, including temperature, in a system may differ according to its previous history [3][4][5][6][7][8][9] . For example, a system that is being subject to cooling and then is suddenly put in contact with a constant temperature source may continue cooling down temporarily or, on the contrary, undergo a temperature rebound 10 . According to current knowledge, there can be only one "spurious" peak/low. In order to check the validity of this, we study here the thermal evolution of a low density set of macroscopic particles that do not preserve kinetic energy upon collision, i.e., a granular gas 11,12 . Our results prove that several extrema may be observed, with the granular temperature moving above and below its steady value. Furthermore, we describe a general mechanism that signals the emergence of these complex memory effects. We also show a route for the observation of extremely large memory effects, as large as the order of magnitude of the stationary value of the variable of interest.
Experimental observations reveal that the response to an excitation of complex condensed matter systems may depend on the entire system's history, and not just on the instantaneous value of the state variables 3,4,13,14 . This is usually called memory effect. Memory effects signal the breakdown of the thermodynamic (or hydrodynamic or macroscopic, depending on the physical context), description. Some typical memory effects include shape memory in polymers 6 , aging and rejuvenation in spin glasses 15 , active matter 16 , and polymers 17 , and the counterintuitive Mpemba effect 18 .
One of the most relevant memory effects related to thermal processes was originally observed by Kovacs and collaborators 3 in a polymer system, which was subject to quenching to a low temperature T 1 from an equilibrium state at temperature T 0 > T 1 . After a long enough waiting time t w but still relaxing towards equilibrium at T 1 , the temperature was suddenly increased back to an intermediate value T s , T 1 < T s < T 0 , such that the instantaneous value of the volume V(t = t w ) equalled the equilibrium value V s corresponding to T s . Subsequently, the volume V(t) did not remain flat but followed a non-monotonic evolution. This nonmonotonic behavior, denominated later as Kovacs hump, consists in reaching one maximum before returning to its equilibrium value V s . We have described above the typical cooling procedure, but also a heating protocol can be considered Quite recently, Kovacs-like memory effects have been thoroughly investigated in glassy systems 19,20 , granular fluids 10 , active matter 21 and disordered mechanical systems 22 . Also, the Kovacs memory effect has been rigorously studied in the context of linear response theory 23,24 . Therein, it has been shown that the normal Kovacs response, i.e., the one originally observed by Kovacs 3 and described just above, with a single maximum, is the only possibility in molecular systems, for which the fluctuation-dissipation theorem holds. However, in athermal systems the response may be anomalous 10,21 : within a region of the system parameters space, a unique temperature minimum can be observed instead of a maximum in a cooling protocol (conversely for a heating protocol). Moreover, the observed deviation of the lone extremum away from the stationary value is typically several orders of magnitude smaller than the stationary value 3,10,20,21 .
One of the main aims of our work is to show that the actual memory effects landscape is in general far more complex. First, several extrema may appear, instead of only one, in a single heating/cooling protocolà la Kovacs. Second, very large memory humps, of the order of magnitude of the stationary value of the magnitude of interest, can be observed. To the best of our knowledge, both features have not yet been reported in the literature.
We consider a collection of identical solid spheres at low particle density so that collisions are always instantaneous and binary but inelastic, i.e., energy is not conserved and we deal with a granular gas 11 . We consider in this case particles with homogeneous mass density and employ the rough hard sphere collisional model with constant coefficients of normal and tangential restitution, α and β, respectively, which is quite realistic for a variety of materials at low particle density 25 .
Let us discuss first why the granular gas of rough spheres is a good candidate for eventually finding complex memory effects. Memory effects appear always in complex systems that consist of many structural units, for which a continuum description seems in principle appropriate. Within this kind of description, the instantaneous value of the complete set of macroscopic variables completely characterizes the system's time evolution 26 . However, there are states that cannot be completely described only with the system macroscopic variables, and it is precisely for these states where a memory effect can be observed. As a matter of fact, this kind of distinct arXiv:1808.08401v1 [cond-mat.soft] 25 Aug 2018 states for which the macroscopic description fails are theoretically very well understood in the context of the kinetic theory of gases 26 . Furthermore, the granular gas of rough spheres can have extremely long relaxation times before it falls into a state where the macroscopic description is valid 27,28 , giving room to the emergence of eventual long lasting memory effects. And, most importantly, in this kind of system there are always two intrinsic, independent and potentially large temperature scales, the translational and rotational granular temperatures, with a highly non-linear coupling. All these facts open new spaces in the search of new important features in complex memory effects, including eventually multiple extrema.
To keep things simple, we consider the granular gas to be in a spatially homogeneous state at all times. The translational velocities are denoted by v, while the angular (or rotational) velocities are denoted by ω. The system is thermalized by a stochastic but homogeneous volume force 29,30 F wn , characterized by a noise intensity χ 2 0 , see Methods section. The kinetic description of our system starts from the corresponding Boltzmann-Fokker-Planck equation for the granular gas under this kind of forcing 28 , see also Methods. The exact solution to this kinetic equation can be formally expressed by means of an expansion around a Maxwellian distribution with variances T t (translational temperature) and T r (rotational temperature) in the translational and angular velocities, respectively. The total granular temperature is given by T = (T t + T r )/2, which is proportional to the average kinetic, translational plus rotational, energy. By adopting a dimensionless time scale τ , proportional to the number of collisions per particle, see Methods, the evolution equations for the temperatures can be written as Above, θ ≡ T r /T t is the temperature ratio and where 3 is a dimensionless measure of the noise intensity, with n being the particle density. The reduced collisional moments µ 20 and µ 02 , see Methods for more reference, are functionals of the whole velocity distribution and therefore the above system of equations is not closed. In order to solve it, we use the first Sonine approximation, which refers to the first non-trivial truncation of the aforementioned exact infinite expansion 27 . For this, together with equations (1), we need to incorporate the evolution equations for the fourth-order cumulants and the initial values of γ, θ and these cumulants, see Methods.
We generate a common initial state for all the temperature evolution curves we subsequently analyze. At an arbitrary time, which we choose to be the time origin τ = 0, and The granular gas is prepared in an initial state (τ = 0) for which all the energy is concentrated in the translational degrees of freedom, as described in the text. In a first stage, 0 < τ < τw, the granular gas freely cools. Then, at the waiting time τ = τw, the noise intensity is suddenly increased from zero to a value such that the instantaneous temperature T (τw) coincides with the corresponding steady temperature Ts. The curves shown for τ > τw correspond to the so-called normal Kovacs response. Time τ measures the average number of collisions per particle. Note that, in order to visualize the Kovacs effect, the relative deviations of T (τ ) from Ts in the response curves have been magnified by a factor r = 5 for all the protocols, except for the transition one, for which r = 100. All the curves correspond to normal and tangential restitution coefficients α = 0.8 and β = 0, respectively.
over an arbitrary previous microscopic state, we apply an instantaneous thermal pulse to the granular gas. In this way, the rotational modes (T r ) of the granular gas are quenched, whereas the translational modes (T t ) are subject to a large heating. As a result, the initial kinetic energy is concentrated in the translational modes, so that the total initial temperature is T (0) = T t (0)/2 and the temperature ratio is θ(0) = 0. Moreover, all the fourth-order cumulants vanish because the initial distribution that results from the heat pulse is a twovariate (T r , T t ) Maxwellian. By this procedure, the system forgets all the previous thermal history of the system, assuring always the same non-equilibrium initial state. From the initial state we have just characterized, the granular gas is left to cool freely, due to the intrinsically inelastic particle collisions 11 , for a waiting time τ w . At τ = τ w , we suddenly apply the stochastic force, with an intensity such that the corresponding steady temperature T s to be reached equals the instantaneous temperature value at the moment of turning the noise on, i.e., T s = T (τ w ). If T (τ > τ w ) further departs from T s , then a Kovacs-like memory effect is observed. What we call protocol is the thermal procedure that we have just described. Depending on the waiting time τ w for turning the stochastic heating on, the system spans different classes of temperature evolution curves. This is depicted and explained in Figure 1. For the sake of simplicity, we investigate the two limiting cases in Figure 1; i.e., τ w = 0 (heating protocol, HP) and τ w → ∞ (cooling protocol, CP).
In the CP, since the system is left cooling down for a long time, the system is already in the homogeneous cooling state (HCS) 11 at τ w . In the HCS, the temperature T (t) is the only relevant variable and decays in time following Haff's law 11 , whereas the temperature ratio and the fourth-order cumulants are time independent and their values depend only on the parameters α and β 27 . Therefore, the conditions for this protocol at τ w are γ(τ w ) = γ s [(1 + θ s )/(1 + θ HCS )] 3 2 , θ(τ w ) = θ HCS , and the fourth-order cumulants at τ = τ w also equal their HCS values. In the HP, the initial conditions for the Kovacs experiment are different. Since we turn on the stochastic force right after the thermal pulse, the initial conditions are those of a bi-variate Maxwellian. Therefore, we have that γ(τ w ) = γ s (1 + θ s ) 3 2 , θ(τ w ) = 0, and, in addition, all the cumulants vanish at τ = τ w .
We clearly show in Fig. 2 the appearance of very large memory effects. Two data sets from molecular dynamics (MD) simulations of the granular gas for both the HP and the CP, together with their corresponding theoretical predictions, are represented. The temperature humps displayed here, of approximately 100% for the CP and 10% for the HP, are larger by at least two orders of magnitude than previously observed memory effects in athermal systems, which at most range from a few thousandths to a few hundredths of the stationary value of the relevant variable 10,21 . The theoretical curves displayed in Fig.2 have been obtained by means of a bi-variate Maxwellian approximation, in which all the cumulants are assumed to be zero, see Methods. Thus, the essential property driving the giant memory effect here is the existence of two independent temperature scales, translational and rotational, i.e. the breakdown of equipartition as given by the fact that θ = 1.
Let us denote the earliest minimum and maximum in the temperature evolution as T m and T M , respectively. We also define H m ≡ T m /T s − 1 < 0, H M ≡ T M /T s − 1 > 0 accordingly. In Fig. 3 we present contour plots highlighting the regions with large |H m | (HP normal response, CP anomalous response) and H M (HP anomalous response, CP normal response). We also plot the transition line H M = |H m | from normal to anomalous response. Huge Kovacs humps appear, especially in the normal region for the HP, in which the size of reported humps can be as large as 100%, relative to the steady temperature.
In order to characterize and quantify complexity in the thermal response we define the parameter S where H 1 (equal to either H m or H M ) is the magnitude of the earliest extremum. Note that S = 0 if there is only one extremum. Thus, S = 0 is the signature of the emergence of more complex response, i.e. with more than one extremum, in the normal-to-anomalous transition. In the transition region, |S| attains its maximum value, |S| = 1, when both extrema are of the same size and neither dominates. The sign of S ∈ [−1, 1] is equal to that of the earliest extremum, providing further information on the detailed structure of the response. Figure 4 represents S as a function of the restitution coefficients α, β, by solving the system (1). Panels A, C, and D correspond to the HP, whereas panels B and E correspond to the CP. We have highlighted in blue (red) regions with S < 0 (S > 0), whereas all points with "simple" memory behavior, i.e., S = 0, remain white. The complex regions are thin but still occupy noticeable sections of the parameter space, especially taking into account that they fall into ranges of experimental values of α and β commonly present in a variety of materials 25 . In Fig. 4A (HP), we clearly observe two zones rich in complex memory effects. In panel C, the first complex zone is zoomed in. This region is attached to the smooth limit, β ∼ −1, and only displays S > 0 for high inelasticities, up to α = 1/ √ 2. In panel D, the second complex region is zoomed in. Within this region, which is close to the quasielastic limit α ∼ 1, the system displays both S > 0 and S < 0 behavior. In Fig. 4B (CP), only one complex Kovacs region next to the quasielastic limit, inside which S > 0, has been identified. Panel E shows a close-up thereof.
It is important to mention that we have found that all of the details of the complex regions emerge in the theoretical solution only when the cumulants are taken into account. This indicates that the temperatures T t , T r do not explain in full detail by themselves the complexity of memory effects found in the rough granular gas. Let us also point out that we have found for the HP a critical narrow region with a discontinuous transition from S > 0 to S < 0, which is signaled in panel C in Fig. 4 and represented in time evolution curves in Fig. 5C. Three consecutive temperature jumps appear in this critical region there before stabilization in the stationary value. Otherwise, all transitions we have observed are continuous, passing from complex to simple (only one extremum) behavior. The latter can be either the normal behavior of molecular systems 3,15,20 (also present in non-equilibrium systems) or the anomalous behavior exclusive of non-equilibrium systems 10,21,24 . This is appropriately tagged in Fig. 4, in which we have labelled the corresponding normal and anomalous regions. Figure 5 displays the evolution curves of the temperature for the three different Kovacs transitions that we have found, in all cases depicted here for the HP: the S > 0 transition in panel A, the S ≈ 0 transition in panel B (in its inset we show the three consecutive humps), and finally the S < 0 case in panel C. All theoretical curves are compared against the numerical solution of the kinetic equation, obtained by means of the direct simulation Monte Carlo (DSMC) method. The agreement is in general excellent, which once more shows the accuracy of our theoretical approach. Although the size of the humps in the transition regions appear smaller than those in Fig. 2 with simple memory behavior, yet they are of the same order of magnitude as those previously reported in the smooth granular gas 10 .
Our work puts forward a general mechanism for the emergence of significantly large memory effects. Enormous humps can be expected if the time evolution of the system under scrutiny is controlled by at least two independent and comparable in magnitude physical variables (here the translational and rotational temperatures) but with only one (here the total temperature) being relevant for the macroscopic or hydrodynamic description. In addition, complex Kovacs response, with multiple extrema, can be expected if the time evolution of the system depends on several additional relevant variables. Here, these additional variables are the fourth-order cumulants, whose sometimes non-monotonic relaxation 27 probably enhances memory effect complexity.
Therefore, memory effects of this size and complexity can potentially be present in other athermal or molecular systems if there are several variables of comparable magnitude are coupled in their time evolution in non-linear form, even if only a subset thereof is relevant in the macroscopic description. This may be relevant, for instance, in active matter systems, where non-linear effects are important in general 31 . We think our results are also especially significant for future experimental work, since we expect these large memory effects to be measurable in granular dynamics experiments; a thermally homogeneous system may be achieved by means of homogeneous turbulent air fluidization 30 .

ACKNOWLEDGMENTS
The authors thank Prof. J. S. Urbach for fruitful discussion. This work has been supported by the Spanish Agencia Estatal de Investigación Grants (partially financed by the ERDF) No. MTM2017-84446-C2-2-R, MTM2014-56948-C2-2-P (A. L.) and No. FIS2016-76359-P (F. V. R. and A. S.), and also by Universidad de Sevilla's VI Plan Propio de Investigación Grant PP2018/494 (A. P.). Use of computing facilities from Extremadura Research Centre for Advanced Technologies (CETA-CIEMAT), funded by the ERDF is also acknowledged.

AUTHOR CONTRIBUTIONS
The ideas of this work were originally conceived by A. L. and F. V. R. and further developed by the four authors. The paper was written by all authors. A. L. and F. V. R. obtained all computer simulations and theoretical data. A. S. wrote the symbolic code for theoretical system's trajectories.  The equations for Tt and Tr have been obtained by multiplying both sides of equation (1) by the translational and rotational kinetic energies, respectively, and integrating over all particle velocity values. The parameters ξt and ξr are respectively. In general, neither ξt nor ξr does have a definite sign, whereas the cooling rate is always positive because energy is dissipated in collisions.
To proceed further, it is convenient to go to dimensionless variables. Time is measured in a scale τ , which is roughly proportional to the number of collisions, because ν(t) is the collision frequency. Dimensionless velocities are introduced as a reduced velocity distribution function as and the dimensionless collision kernel as In dimensionless variables, the evolution equations for the temperatures can be written as equations (1) in the main text. Therein, there appear the reduced collisional moments µ 20 ≡ µ Note that, aside from the nondimensionalizing factors, the production rates ξt and ξr are basically identical to µ 20 and µ 02 , respectively. These are functionals of the whole distribution function and thus the evolution equations for the temperatures are not closed. In order to close the dynamical equations, a formally exact expansion in orthogonal polynomials can be performed 27 . For isotropic states, we can expand the velocity distribution around the Maxwellian φ M (c, w) = π −3 e −c 2 −w 2 , where Ψ ( ) jk (c, w) are certain products of Laguerre and Legendre polynomials. By normalization, a which we call the fourth-order cumulants henceforth.
Maxwellian approximation.-The simplest description is obtained by substituting the Maxwellian velocity distribution into the collision integrals (11). Equivalently, one may consider that all the non-trivial cumulant vanish in this approach, which yields where κ ≡ 4(mσ 2 ) −1 I is the dimensionless moment of inertia. Insertion of equation (14) into the evolution equations (1) in the main text gives rise to the Maxwellian approximation. First Sonine approximation.-A more elaborate approximation can be done by incorporating the lowest order cumulants, which we defined in equations (13), as the first corrections to the Maxwellian. A closed set of six coupled differential equations can be obtained for θ(τ ), γ(τ ), a  The resulting set of six differential equations can be numerically solved with appropriate initial conditions for each physical situation, as discussed in the main text. In this way, we obtain the time evolution of the temperatures in the so-called first Sonine approximation, to which we refer throughout this work.
Computer simulations. We use in this work data sets obtained from computer simulations from two independent and different methods: Direct Simulation Monte Carlo method (DSMC), which obtains an exact numerical solution of the relevant kinetic equation, in our case Eq. (1), and molecular dynamics (MD) simulation, which solves particles trajectories. A detailed description of the DSMC method may be found elsewhere 32 . In our DSMC simulations, and in order to reduce statistical noise in the temperature time evolution curves, we have used an average of 100 statistical replicas of a system with 2 × 10 6 particles. In the MD case, we have simulated 1000 inelastic hard spheres at a density nσ 3 = 0.01 and averaged over 500 trajectories.