Linear response in the uniformly heated granular gas

PHYSICAL REVIEW E 104, 024903 (2021) Linear response in the uniformly heated granular gas Bernardo Sánchez-Rey Departamento de Física Aplicada I, E.P.S., Universidad de Sevilla, Virgen de África 7, E-41011 Sevilla, Spain Antonio Prados Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain (Received 21 October 2020; revised 7 June 2021; accepted 29 July 2021; published 23 August 2021) We analyze the linear response properties of the uniformly heated granular gas. The intensity of the stochastic driving fixes the value of the granular temperature in the nonequilibrium steady state reached by the system. Here, we investigate two specific situations. First, we look into the “direct” relaxation of the system after a single (small) jump of the driving intensity. This study is carried out by two different methods. Not only do we linearize the evolution equations around the steady state, but we also derive generalized out-of-equilibrium fluctuation-dissipation relations for the relevant response functions. Second, we investigate the behavior of the system in a more complex experiment, specifically a Kovacs-like protocol with two jumps in the driving. The emergence of an anomalous Kovacs response is explained in terms of the properties of the direct relaxation function: it is the second mode changing sign at the critical value of the inelasticity that demarcates anomalous from normal behavior. The analytical results are compared with numerical simulations of the kinetic equation, and a good agreement is found. DOI: 10.1103/PhysRevE.104.024903 I. INTRODUCTION Linear response is at the root of many important results in physics. In this respect, the fluctuation-dissipation theorem is a milestone in the development of (nonequilibrium) statistical physics [1–3]. In its original form, it relates the linear relaxation of a system to equilibrium, from an initial nonequilibrium state, with certain equilibrium time correlation functions. Very recently, this result has been extended to more general situations, which include the relaxation to nonequilibrium steady states (NESSs) [4]. The above general picture makes it relevant to investigate the linear relaxation of intrinsically out-of-equilibrium systems, such as granular gases [5–8], to their NESS. Due to the energy dissipation in collisions, an external energy input mechanism is needed to drive the system to a NESS. One of the simplest physical situations is that of the uniformly heated granular gas, in which all the particles of the gas are submitted to independent white noise forces of a given amplitude [9,10]. Therein, the granular gas remains homogeneous for all times, if it was initially so, and the main physical property of the granular gas is the granular temperature—basically, the average kinetic energy per degree of freedom. It is the amplitude of the white noise force, i.e., the intensity of this stochastic thermostat, that determines the stationary value of the granular temperature. The main purpose of this work is to analyze the linear relaxation of the granular temperature in the uniformly heated granular gas, in several different physical situations. Due to the non-Gaussian character of the velocity distribution function, the velocity moments obey an infinite hierarchy of coupled ordinary differential equations (ODEs), which must be supplemented with a closure assumption. Here, we work in the first Sonine approximation [11], in which the state of the granular system is characterized by the granular temperature T, and the second Sonine coefficient a2—the excess kurtosis. The dynamics of the gas is described by a set of two ODEs for the temperature and the excess kurtosis [12], first derived in Ref. [13] for the free cooling case. These evolution equations are nonlinear and have been shown to describe accurately the granular gas in many different physical situations; see, for example, [9,10,12,14–17]. First, we would like to investigate the response of the gas to an instantaneous perturbation of the driving intensity. This is a relevant problem: different memory effects have been recently reported in the uniformly heated granular gas [16–18]. In principle, the existence of these memory effects suggests that the linear relaxation of the granular temperature is nonexponential. Indeed, nonexponential relaxation [19–25] is a key ingredient for the emergence of aging and memory effects [26–37] in many different physical contexts. In light of the above, it is worth elucidating the relaxation behavior of the granular temperature in this single-jump experiment. More specifically, we would like to clarify its exponential or nonexponential character, and also its monotonicity properties. In this respect, it is important to clarify the role played by the inelasticity. For example, a monotonic decay of the direct relaxation function ensures that the Kovacs effect is “normal,” i.e., the hump has the same sign as in molecular systems [38]. For small enough perturbations, the evolution equations for Tand a2are linearized and analytically solved. Moreover, the 2470-0045/2021/104(2)/024903(13) 024903-1 ©2021 American Physical Society BERNARDO SÁNCHEZ-REY AND ANTONIO PRADOS PHYSICAL REVIEW E 104, 024903 (2021) relaxation of both quantities to their steady values is shown to be directly related to some time correlation functions (calculated in the NESS), by extending the ideas in Ref. [4]for an out-of-equilibrium fluctuation-dissipation relation (FDR). The correlation functions involve the derivative of the N- particle velocity distribution function, and thus we introduce a factorization assumption—sometimes called propagation of chaos [39]—to get an explicit expression for the relevant time correlations. Our analytical predictions are compared with direct simulation Monte Carlo (DSMC) results—i.e., the numerical integration of the kinetic equation. This is done for the two procedures described above. In the “direct” route, the system is initially put in the NESS corresponding to a certain value ξ0of the driving, the driving is instantaneously changed to ξat t=0+, and subsequently the time evolution of the granular temperature is recorded. In the FDR route, the system is initially put in the NESS corresponding to ξ, and the corresponding time correlation function is evaluated in this NESS—there is no need to change the driving. More detailed insight into the dynamics of the system can be acquired with more complicated driving protocols. A particularly relevant one is the so-called Kovacs experiment, which has been extensively analyzed in the realm of glassy systems [40–52]. The Kovacs protocol involves two jumps of the parameter ξ—temperature, driving intensity, etc.—controlling the relaxation of the system. First, it is changed from ξ0to ξ1<ξ 0at t=0+and the system relaxes towards the steady—either equilibrium or NESS—state for ξ1. This relaxation lasts for a waiting time tw: it is interrupted when the value of the “thermodynamic” property of interest equals the steady-state value for some intermediate value of ξ,att=t+ wthe value of the control is changed to ξ(ξ1<ξ<ξ 0). If the subsequent behavior is nonmonotonic, with the thermodynamic property departing from its steady value before returning thereto, additional variables are needed to completely characterize the steady state of the system. On the one hand, the vast majority of the studies of the Kovacs effect have been done in the nonlinear regime, i.e., for large jumps of the driving, and thus the results are mainly numerical. On the other hand, a general theory of the Kovacs effect is only available within the limits of applicability of linear response theory. For molecular systems—such that the steady distribution is the canonical one—with Markovian dynamics, it has been shown that the Kovacs hump has some general characteristic features that stem from the standard version of the FDR [47]. Specifically, the Kovacs hump has only one extremum and cannot change sign, i.e., it is always normal. In athermal systems, there appear Kovacs responses that deviate from the normal behavior just described. The Kovacs effect is said to be anomalous when only one extremum is present but the sign of the hump changes with the system parameters. This behavior has been observed both in the uniformly heated granular gas [16,18] and in active matter models [53]. Also, more than one extremum has been reported in the rough granular gas [54]. It must be stressed that all these observations correspond to the nonlinear regime. In particular, in the uniformly heated granular case, the extreme case ξ1=0 (i.e., no driving in the waiting time window) was considered [16,18]. The emergence of these more complex Kovacs responses, including the anomalous behavior, is not well understood on a general basis. Although the existence of anomalous behavior is supported by analytical calculations in specific systems, the general physical reason behind it is not known. One of the main objectives of this work is to shed light on this point by considering the linear response limit in the uniformly heated granular gas. To the best of our knowledge, the connection between the behavior of the direct relaxation function and the Kovacs hump has not been analyzed in the context of granular gases. Here, we show that the anomalous behavior survives in linear response and thus we can explain its emergence in terms of the behaviour of the one-jump relaxation function. This is done by bringing to bear recent results on the generalization of the mathematical structure of the linear Kovacs hump to athermal systems [38,53]. The organization of the paper is as follows. In Sec. II, we put forward the evolution equations for the granular temperature and the excess kurtosis, and we carry out their linearization around the NESS in Sec. IIA. The linear relaxation of the granular temperature to an instantaneous change of the driving is considered in Sec. III for different values of the inelasticity. The generalized FDR for a jump in the driving is the subject of Sec. IV. Specifically, it is derived in Sec. IV A and particularized for the relevant response functions in Sec. IV B. Section Vpresents DSMC results for the relaxation of the granular temperature and the excess kurtosis—or, alternatively, the fourth moment of the velocity. Next, Sec. VI is devoted to analyzing the Kovacs experiment. Finally, conclusions and a brief discussion of perspectives for future work are presented in Sec. VII. II. EVOLUTION EQUATIONS Let us analyze the dynamics of a granular gas of smooth hard spheres [7]. This system comprises Nhard particles of mass mand diameter σin dimension d, which undergo inelastic collisions. In the binary collisions, the tangential component of the relative velocity remains unaltered while the normal component is reversed and shrunk by a factor α, the normal restitution coefficient, 0 ⩽α⩽1. Apart from the elastic limit α=1, kinetic energy is lost in every collision and the undriven system “cools,” in the sense that its granular temperature—basically the average of the kinetic energy— decreases monotonically in time. We consider the uniformly heated granular gas, i.e., the system described above is also submitted to a random forcing that inputs energy into the system. This is modelled as independent white noise forces F(p) i(t) acting over each grain, F(p) i(t)=0 and F(p) i(t)F(q) j(t)=m2ξ2δijδpqδ(t− t), i,j=1,...,d,p,q=1,...,N. As a consequence, the system reaches a steady state in the long time limit, in which the energy input by the thermostat cancels—in average—the energy loss in collisions [9]. At the kinetic level of description, the dynamical evolution of the system is governed by the Boltzmann-Fokker-Planck equation for the velocity distribution function [9,10]. Because 024903-2 LINEAR RESPONSE IN THE UNIFORMLY HEATED … PHYSICAL REVIEW E 104, 024903 (2021) the granular gas is an intrinsically nonequilibrium system, its velocity distribution function is non-Gaussian, even in the steady state. The granular temperature T(t) is defined as d 2T(t)≡1 2mv2(t).(1) To measure the departure from the Maxwellian distribution, the simplest approach is to consider the excess kurtosis a2, a2=d d+2v4 v22−1,(2) which vanishes for the Gaussian case. In kinetic theory, working with the pair (T,a2) is known as the first Sonine approximation, because a2is the first nonzero coefficient of the expansion of the velocity distribution function in a series of Sonine polynomials [9,11,15,55], P1(v;t)=vT(t)−dπ−d/2e−w2[1 +a2(t)S2(w2)],(3) where wis defined by w=v vT(t),v2 T(t)=2T(t) m=2 dv2(t),(4) and S2(x)=1 2x2−d+2 2x+d(d+2) 8.(5) The evolution equations for the granular temperature and the excess kurtosis are derived from the Boltzmann-Fokker- Planck equation [9,12,13,15,16]. They are usually written in terms of the collision rate ζ0, and the stationary values of the granular temperature Tsand the excess kurtosis as 2, Ts=mξ2 ζ01+3 16 as 22/3 ,ζ 0=2nσd−1(1−α2)πd−1 2 √md(d/2) , (6a) as 2=16(1 −α)(1 −2α2) 73 +56d−24dα−105α+30(1 −α)α2.(6b) By defining scaled—order of unity—variables as follows: θ=T Ts ,A2=a2 as 2 ,τ=ζ0√Ts 2t.(7) we get dθ dτ=21−θ3/2+3 16as 2(1 −A2θ3/2),(8a) θdA2 dτ=4[(θ3/2−1)A2+Bθ3/2(1 −A2)],(8b) with the parameter [56] B=73 +8d(7 −3α)+15α[2α(1 −α)−7] 16(1 −α)(3 +2d+2α2).(9) A couple of comments on the evolution equations (8)are appropriate. First, they are nonlinear, so in general it is not possible to write down an analytical closed solution for them [57]. Second, the parameters Band as 2only depend on αand d, being independent of the steady value of the temperature. In fact, all the dependence on Tshas been incorporated into the timescale τintroduced in Eq. (7). It will be useful for our purposes to define the dimensionless velocity c=v vs T ,vs T=2Ts m,(10) the second moment of which is directly related to the dimensionless temperature, c2= v2 vs T2=d 2θ, c2s=d 2.(11) Equation (2) still holds with the change v→c, whereas the fourth moment of cat the steady state is given by c4s=d+2 dc22 s1+as 2=d(d+2) 41+as 2.(12) A. Linearization around the steady state In this paper, we are interested in the linear relaxation of the granular gas to the steady state, in which the temperature and the excess kurtosis are θs=1 and As 2=1, as a consequence of the scaling. Therefore, we write θ=1+δθ, A2=1+δA2,(13) insert them into Eq. (8), and neglect nonlinearities. The following linear system is thus obtained: d dtδz=M·δz,(14a) δz≡δθ δA2,M≡−31+3 16 as 2−3 8as 2 6−4B.(14b) The matrix—or operator [58]—Mdoes not depend on the noise intensity ξbut only on the restitution coefficient αand the dimension d. This is a consequence of the right-hand side of the nonlinear system (8) being independent of Ts, as already pointed out above. In the following, we denote the elements of the matrix M by Mij, i.e., M11 =−31+3 16as 2,M12 =−3 8as 2, M21 =6,M22 =−4B.(15) The eigenvalues of the matrix Mare both negative, λ±=TrM±(TrM)2−4 det M 2<0,(16) and the corresponding eigenvectors read ζ+=M12 λ+−M11,ζ−=M12 λ−−M11.(17) The explicit expression for the trace and the determinant of M are TrM=−4B−31+3 16as 2<0,(18a) det M=12B1+3 16as 2+9 4as 2>0.(18b) 024903-3 BERNARDO SÁNCHEZ-REY AND ANTONIO PRADOS PHYSICAL REVIEW E 104, 024903 (2021) The general solution of the linear system is δz(τ)=C+ζ+eλ+τ+C−ζ−eλ−τ,(19) or, equivalently, δθ(τ) δA2(τ)=C+M12 λ+−M11eλ+τ+C−M12 λ−−M11eλ−τ.(20) The constants C+and C−are determined by the initial conditions through the relation δθ(0) δA2(0)=δz(0) =C+ζ++C−ζ−.(21) III. DIRECT RELAXATION AFTER A SINGLE JUMP IN THE DRIVING Let us consider an experiment in which the system is initially at the steady state corresponding to some value of the driving ξ+δξ. At the initial time, the intensity of the driving is suddenly changed from ξ+δξ to ξ: subsequently, the granular gas relaxes from its nonequilibrium steady state (NESS) at ξ+δξ to the corresponding NESS at ξ. The initial values of the scaled variables are δθ(0) = 0,δA2(0) =0,(22) because the steady value of the temperature depends on the driving intensity, but that of the excess kurtosis does not. Equation (21) predicts that both C+and C−are proportional to δθ(0) in this case, so we write C±=δθ(0) ˜ C±(23) and readily obtain ˜ C+=M11 −λ− M12(λ+−λ−),˜ C−=λ+−M11 M12(λ+−λ−).(24) Therefore, we have that δz(τ) δθ(0) =˜ C+ζ+eλ+τ+˜ C−ζ−eλ−τ.(25) In what follows, we denote the relaxation function for the physical property Yby φY(τ), specifically φY(τ)=δY(τ) δθ(0) .(26) First, we analyze the relaxation of the granular temperature, i.e., φθ(τ)≡δθ(τ) δθ(0) =a+eλ+τ+a−eλ−τ,(27) with a+=M11 −λ− λ+−λ− ,a−=λ+−M11 λ+−λ− .(28) In Fig. 1, we show the two coefficients a+and a−as a function of the restitution coefficient α,ford=2 and 3. Both for the two- and the three-dimensional case, two distinct properties are observed: (i) while a+>0 for all values of α, a−changes sign at the critical value αc=1/√2, specifically sgn(a−)=sgn(α−αc),α c=1/√2,(29) FIG. 1. Coefficients a±as a function of the restitution coefficient α. Both the three-dimensional (solid line) and the two-dimensional (dashed line) cases are plotted. The coefficient a+is always positive, whereas a−changes sign at the critical value α=αc—as given by Eq. (29). Note that |a−|is quite small throughout, with its least small value at α=0, and thus a+=1−a−is always very close to unity. and (ii) the ratio |a−/a+|1 for all values of α, with its larger value taking place at α=0, for which it equals 0.0159 (0.0114) for d=2(d=3). In principle, the change of sign of a−at αcmight bring about a nonmonotonic behavior of the relaxation function, but this is not the case. In the long-time limit, the relaxation is dominated by the largest eigenvalue λ+ and, therefore, φθ(τ)∼a+eλ+τ>0,(λ+−λ−)τ1.(30) To understand the above behavior of a±, we bring to bear the smallness of as 2and expand all quantities in powers of as 2. In particular, the eigenvalues λ±are given by λ+=−3−9(1 +4B) 16(4B−3)as 2+Oas 22,(31) λ−=−4B+9 4(4B−3)as 2+Oas 22.(32) Note that 4B−3>0(λ+>λ −). In turn, the corresponding expansion for the coefficients a±is a+=1+9 4(4B−3)2as 2+Oas 22,(33) a−=− 9 4(4B−3)2as 2+Oas 22.(34) 024903-4 LINEAR RESPONSE IN THE UNIFORMLY HEATED … PHYSICAL REVIEW E 104, 024903 (2021) On the one hand, it is neatly observed that a−is proportional to as 2, thus being rather small across the whole range of inelasticities. Moreover, the change of sign of a−at αcis clearly linked to the vanishing of as 2thereat. On the other hand, a+=1+O(as 2) is very close to unity throughout, consistent with the normalization condition a++a−=1. Additional information is given by the relaxation of the excess kurtosis or, equivalently, of the fourth moment c4. Equation (25) implies φA2(τ)≡δA2(τ) δθ(0) =a+a−(λ+−λ−) M12 (eλ+τ−eλ−τ).(35) The coefficient in front of the difference of exponentials can also be expanded in powers of as 2, making use of Eqs. (15) and (31)–(34), with the result a+a−(λ+−λ−) M12 =6 4B−3+27(5 +4B) 8(4B−3)3as 2+Oas 22. (36) The dominant behavior 6/(4B−3) is positive for all α:this means that δA2(τ)hasthesamesignasδθ(0) for all times, i.e., the same sign of δξ [59]. Note that, in our experiment, the driving is instantaneously decreased from ξ+δξ to ξat the initial time: therefore, δξ > 0 and δA2(τ)>0. The relaxation function φc4(τ) can be written as a linear combination of φθand φA2. Bringing to bear Eqs. (2), (10), (27), and (35), one gets δc4=d(d+2) 21 2as 2δA2+1+as 2δθ(37) for the deviation of c4from its steady-state value, given by Eq. (12). The relaxation function thus reads φc4(τ)=d(d+2) 21 2as 2φA2(τ)+1+as 2φθ(τ).(38) IV. GENERALIZED FDR FOR A JUMP IN THE DRIVING The relaxation functions above can be related to certain time correlations in the NESS of the granular gas by means of a generalized FDR. In the following, we first derive a generalized FDR for the relaxation of an arbitrary function of the velocities, and afterwards we particularize it to get the specific relations for the relaxation functions considered in the previous section. A. Derivation of the generalized FDR Now we investigate the evolution of the granular gas, after an instantaneous jump in the driving like the one considered in the previous section. At the N-particle level, the fluid is described by PN(,t), which is the (one-time) probability density for finding the system with =(v1,v2,...,vN)at time t. For homogeneous situations—like those considered throughout this paper, is a Markov process described by Kac’s equation [39,60], and for any function of the velocities we can write f(t)=dtf(t)PN(t,t) =dtd0f(t)Tt(t|0)PN(0,0),(39) where Tt(t|0) is the transition probability from 0to tin a time interval t. Thus, we have PN(,t,,0) =Tt(|)PN(,0) (40) for the two-time probability density for finding the N-particle granular gas with =(v1,v2,...,vN) at time tand with =(v 1,v 2,...,v N) at time t=0. Now we consider that the initial distribution at t=0 corresponds to the steady state for T(0) =Ts+δT(0), PN(0,0)=Ps N(0)+δT(0) ∂Ps N(0) ∂Ts .(41) Inserting this equation into Eq. (39), one gets δf(t)=δT(0)dtd0f(t)T(t|0)∂Ps N(0) ∂Ts , (42) in which we have defined, consistently with the notation we employ throughout, δf(t)≡f(t)−f()s(43) for the deviation of the average value of ffrom its steady-state value. Equation (42) can be rewritten as δf(t)=δT(0)f(t)∂lnPs N(0) ∂Tss ,(44) from which the time correlation function Cf(t)≡f(t)∂lnPs N(0) ∂Tss =dtd0f(t)∂lnPs N(0) ∂Ts T(t|0)Ps N(0) (45) emerges. Therefore, we can also obtain the relaxation functions by measuring correlations at the NESS, which is the expression of the generalized FDR. The subindex sin ···s stresses that the average in the time correlation function is done in the steady state, because T(t|0)Ps N(0)=Ps N(t,t,0,0) (46) is the two-time probability density at the steady state with temperature Ts[61]. The difficulty of this approach is the necessity of knowing the N-particle distribution Ps N(0): note that one deals with the one-particle distribution function in the kinetic approach. This difficulty is usually circumvented by introducing the following factorization ansatz [4]: Ps N()= N  i=1 P1(vi),(47) which is sometimes called “propagation of chaos” [39]. This is more restrictive than the usual molecular chaos hypothesis employed to derive the Boltzmann equation: therein, only 024903-5 BERNARDO SÁNCHEZ-REY AND ANTONIO PRADOS PHYSICAL REVIEW E 104, 024903 (2021) the two-particle distribution is assumed to factorize. With the assumption (47), Cf(t) simplifies to Cf(t)= N  i=1 Ci f(t),Ci f(t)≡f(t)∂lnPs 1(vi0) ∂Tss ,(48) where Ps 1is the stationary solution of the Boltzmann equation, which is—at least approximately—known. B. Response functions Here, we derive the specific formulas relating the relaxation functions φθand φc4to time correlations. To achieve this goal, we need ∂TslnPs N() and thus we use the propagation of chaos assumption (47). Consistent with our approach, we employ the one-particle velocity distribution function in the first Sonine approximation. Particularizing (3) to the steady state, we have Ps 1(v)=vs T−de−c21+as 2S2(c2),(49) where cis defined in Eq. (10). Taking into account that ∂c2/∂Ts=−c2/Ts, we get ∂lnP1(v) ∂Ts=−d 2Ts+1 Ts c2−1 Ts c2as 2S 2(c2) 1+as 2S2(c2).(50) By defining g(x)=−d 2+x−xas 2S 2(x) 1+as 2S2(x),(51) we can write ∂ln P1(v) ∂Ts=1 Ts g(c2),Ci f(t)=1 Tsf(t)gc2 i0s,(52) and thus the response function is given by δf(t)=δθ(0) N  i=1f(t)gc2 i0s.(53) Now we particularize the above general relation to the specific cases in which we are interested. Note that the dimensionless time τis then used, instead of t. First, we consider the following choice for the function f: f(τ)=1 N N  i=12 dc2 i(τ)−1,f(τ)=δθ(τ).(54) Making use of Eq. (53), we have for the relaxation function of the granular temperature φθ(τ)=R2(τ)≡1 N N  i=1 N  j=12 dc2 i(τ)−1gc2 j(0)s .(55) Second, we consider the relaxation of the fourth moment. Therefore, now we choose fto be f(τ)=1 N N  i=1c4 i(τ)−c4s,f(τ)=δc4(τ), (56) FIG. 2. Relaxation function of the granular temperature as a function of the scaled time τ. Panels (a)–(c) display φθ(τ) for three different values of the restitution coefficient, namely α=0.3, 0.7, and 0.9, respectively—symbols correspond to simulation results in a system of 106particles averaged over 10 runs, and lines to the analytical prediction (27). The three simulation curves are plotted together in panel (d), where it is clearly seen that they superimpose. Therefore, the linear relaxation of the granular temperature to the corresponding NESS is “universal” in the τscale. where c4sis given by Eq. (12). Again, we employ (53)to write φc4(τ)=R4(τ)≡1 N N  i=1 N  j=1c4 i(τ)−c4sgc2 j(0)s. (57) V. NUMERICAL SIMULATIONS Now we compare our analytical predictions for the relaxation function with numerical data obtained from the DSMC numerical integration of the Boltzmann-Fokker-Planck equation. We have employed a system of N=106hard disks (d=2) of unit mass m=1 and unit diameter σ=1, with the binary collision rule between particles iand j, v i=vi−1+α 2(ˆ σ·vij)ˆ σ,vj=vj+1+α 2(ˆ σ·vij)ˆ σ. (58) Above, (vi,vj) are the precollisional velocities, (v i,vj)are the postcollisional ones, vij ≡vi−vjis the relative velocity, and ˆ σis the unit vector pointing from the center of particle jto the center of particle iat the collision. Moreover, in order to simulate the stochastic thermostat, the hard disks are submitted to random kicks every Nc=500 collisions. In the kick, each component of the velocity of every particle is incremented by a random number extracted from a Gaussian distribution of variance ξ2t, where tis the time interval corresponding to the number of collisions Nc. A. Relaxation function of the temperature In particular, we show the relaxation functions corresponding to α=0.3(<αc), α=0.7(≃αc), and α= 0.9(>αc)—see panels (a)–(c) of Fig. 2. The symbols 024903-6 LINEAR RESPONSE IN THE UNIFORMLY HEATED … PHYSICAL REVIEW E 104, 024903 (2021) FIG. 3. Relaxation rates λ±as a function of the restitution coefficient. Again, both d=3 (solid line) and d=2 (dashed) are plotted. The largest rate λ+is always close to −3 throughout. The smaller rate λ−has a stronger dependence on α, diverging in the elastic limit α→1−. represent simulation results from a NESS at (ξ+δξ)2= 0.205 to the corresponding NESS at ξ2=0.2, while the solid lines represent our theoretical prediction (27). The agreement is excellent, which comes as no surprise because the first Sonine approximation is known to give a good description of the granular gas dynamics—although the majority of the studies have been done in the nonlinear regime [12,16,18]. In principle, Eq. (27) tells us that the relaxation of the granular temperature to the NESS is nonexponential. However, the smallness of the coefficient a−throughout the whole αrange entails that, from a practical point of view, the relaxation is very well fitted by a single exponential for all values of the restitution coefficient. Moreover, the relaxation functions for different values of αsuperimpose, when plotted as a function of the scaled time τ[62]. This is neatly observed in panel (d) of Fig. 2, where we have put together the simulation results from the three previous panels. This stems from the weak α-dependence of the eigenvalue λ+, which is close to −3for all values of α, as predicted by Eq. (31) and shown in panel (a) of Fig. 3. Therein, it is observed that the maximum deviation from −3 is of the order of 3% for the largest value of as 2—i.e., in the completely inelastic limit α=0. For completeness, we show λ−in panel (b) of the same figure. Now we proceed to evaluate the relaxation function of the temperature by using the FDR relation (55). The correlation function on the right-hand side has been numerically evaluated from the DSMC data, employing the system with N=104–106hard disks and averaging over a long trajectory of duration τ=400–6000, depending on the value of α. We compare the numerics with Eq. (55)forφθin Fig. 4. On the one hand, we observe in the left panel that there appears a discrepancy at the initial time. The value stemming from the FDR relation deviates from unity in general: it is larger (smaller) for α=0.3(α=0.9), whereas it equals 1 for α=0.7≃αc. On the other hand, the right panel shows that this discrepancy disappears if we normalize the correlation function R2(τ) with its initial value, i.e., if we compare φθ(τ) with R2(τ)/R2(0). Therefore, the correlation R2(τ) correctly predicts the decay, R2(τ)=R2(0)φθ(τ), but not the initial value. In light of the above, it is worth asking the reason behind the discrepancy for the initial time. The correlation function includes both diagonal terms (i=j) and nondiagonal terms (i= j) in its double sum, thus we split it accordingly into R2(0) =Rdia 2(0) +Rndia 2(0),(59a) Rdia 2(0) =1 N N  i=12 dc2 i−1gc2 is =2 dc2 1−1gc2 1s ,(59b) Rndia 2(0) =1 N N  i=1 j=i2 dc2 i−1gc2 js =(N−1)2 dc2 1−1gc2 2s .(59c) Note that the above one-time averages are done in the NESS and are thus time-independent, so we have written c(0) →c. Employing the factorization assumption (47) and the first Sonine approximation (49)forPs 1, we obtain that Rdia 2(0) =2 dc2−1g(c2)s=1,(60a) Rndia 2(0) =(N−1)2 dc2−1sg(c2)s=0.(60b) We have also evaluated Rdia 2(τ) and Rndia 2(τ) from the DSMC data. We have always got Rdia 2(0) ≃1, whereas in general Rndia 2(0) = 0: the discrepancy at the initial time stems from the nondiagonal terms. The physical reason behind this is the O(N−1) correlations that are completely neglected when writing Eq. (47). To account for this behavior, it is necessary to consider the two-particle time correlations, which is outside the scope of this work [63]. As shown in the right panel of Fig. 4,Rdia 2(τ) decays faster than the whole correlation function R2(τ); the relaxation time of the latter approximately doubles that of the former. It must be remarked that the discrepancy between the relaxation function and the correlation is not a consequence of the first Sonine approximation. One can easily incorporate the contribution from higher cumulants into R2(τ), for instance the following term that includes the sixth-cumulant a3and the third Sonine polynomial S3(x). The only change is that now 024903-7 BERNARDO SÁNCHEZ-REY AND ANTONIO PRADOS PHYSICAL REVIEW E 104, 024903 (2021) FIG. 4. Check of the FDR for the relaxation function of the granular temperature. The left panel shows the correlation function R2(τ)for three different values of the restitution coefficient α: from top to bottom, α=0.3 (squares), 0.7 (circles), and 0.9 (triangles). Also shown is the theoretical prediction for the direct relaxation function of the granular temperature φθ(τ) (line). As discussed in the text, the decay of φθ is well described by R2, although the initial values do not exactly match. This is illustrated in the right panel, where we plot R2(τ)/R2(0) for the specific case α=0.3. Therein, R2(τ)/R2(0) is compared with both φθ(τ) (solid line) and the diagonal part of the correlation Rdia 2(τ) (dashed line). It is clearly observed that the diagonal part equals unity for τ=0, but it decays faster than the whole correlation R2: roughly, its relaxation time is divided by 2. g(x) reads g(x)=−d 2+x−xas 2S 2(x)+as 3S 3(x) 1+as 2S2(x)+a3S3(x).(61) We have checked that the DSMC results for the correlation function remain basically unaltered. B. Relaxation function of the fourth moment Simulation curves for the time evolution of A2are compared with the theoretical prediction (35)inFig.5. All curves correspond to a system with restitution coefficient α=0.3 and they have been averaged over 100 runs. Two experiments in which the driving intensity is decreased from ξ+δξ to ξ are considered. In both cases, we have that ξ2=0.2, but with two different values of the jump, (ξ+δξ)2=0.25 (squares) and (ξ+δξ)2=0.35 (solid circles). FIG. 5. Relaxation function of the scaled excess kurtosis A2= a2/as 2as a function of τfor α=0.3. In the simulations, the driving is instantaneously changed from (ξ+δξ)2=0.25 (squares) and (ξ+ δξ)2=0.35 (solid circles) to ξ2=0.2atτ=0. Despite the jumps being quite large, the agreement with the theoretical prediction, as given by Eq. (35) (line), is remarkably good. Fluctuations in A2are larger than those of the granular temperature, and thus it is necessary to consider larger jumps in the driving, roughly one order of magnitude larger than those employed in Fig. 2. Still, linear response theory works pretty well: for (ξ+δξ)2=0.25, δξ/ξ ≃0.12 but the agreement is almost perfect; for (ξ+δξ)2=0.35, ξ/ξ ≃0.3 and the theory only slightly overestimates the relaxation function φA2(τ). For these larger jumps in the driving, we have also looked into the relaxation of the granular temperature. This is done in Fig. 6; the agreement between the simulation curves and the linear response result is even better than that of the excess kurtosis. We have also checked the accuracy of the FDR relation (57) for the relaxation function φc4of the fourth moment. We have numerically computed the time correlation function R4(τ) from the DSMC data; the values for the simulation parameters are the same as for the temperature. In Fig. 7, FIG. 6. Relaxation function of the granular temperature for the same situation considered in Fig. 5. Note the logarithmic scale on the vertical axis. Although the jumps in the driving are quite large— relative size of 10–30%—the linear response prediction (solid line) fits perfectly well the simulation results (symbols). The coding for the different lines is the same as in Fig. 5. 024903-8 LINEAR RESPONSE IN THE UNIFORMLY HEATED … PHYSICAL REVIEW E 104, 024903 (2021) FIG. 7. Check of the FDR relation for the relaxation function of the fourth-moment. The squares represent the correlation function R4(τ)forα=0.3, measured in DSMC simulations, while the solid line is the theoretical prediction given by Eq. (38). The decay of φc4 is well described by R4, although their initial values do not exactly match. On the other hand, the diagonal part of the correlation Rdia 4(τ) (dashed line) correctly gives the initial value φc4(0), but it decays faster than the whole correlation R4. for α=0.3, the numerics (squares) is compared with the theoretical prediction for φc4, as given by Eq. (38) (solid line). Similarly to the case of the temperature relaxation, the correlation function R4(τ) correctly predicts the decay but not the initial value of φc4. Again, this discrepancy comes out from the nondiagonal terms of the correlation function, as the evaluation of the diagonal part Rdia 4(τ) (dashed line) shows. VI. THE KOVACS-LIKE EXPERIMENT WITH TWO JUMPS IN THE DRIVING In this section, we analyze a more complex experiment, which was first done by Kovacs to investigate the glassy response of polymers [40,41]. The main idea is the following: first, the system under study is “aged” by following a certain protocol. After this aging stage, the system has values of its macroscopic variables equal to those in the steady state. However, the macroscopic variables do not remain flat but display a “hump.” This is the Kovacs effect, which tells us that there are other, nonmacroscopic, variables that also have to be taken into account to understand the relaxation properties of the system. Here, we describe the Kovacs-like experiment in terms of the variables of our model [64]. Instead of letting the system relax directly from the NESS for ξ+δξ to that for ξ,att=0 we suddenly change the driving intensity from ξ+δξ to a lower driving ξ−δξ<ξ. Then, the system starts to relax to the NESS for ξ−δξ, at which Ts(ξ−δξ)=Ts(ξ)−δT< Ts(ξ). At some time—which we call henceforth the waiting time τw—the instantaneous value of the granular temperature T(τ) equals Ts(ξ). At τ=τw, we suddenly change the driving intensity from ξ−δξto ξ. This two-jump procedure in the driving characterizes the Kovacs protocol. A qualitative plot of the Kovacs protocol is shown in Fig. 8. Specifically, we have considered the usual case in which δξ and δξare both positive and therefore Ts+δT>Ts>Ts− FIG. 8. Qualitative plot of the Kovacs protocol. The system starts from a NESS with granular temperature T0=Ts+δT, corresponding to a driving ξ0=ξ+δξ. In the waiting time window, the driving is suddenly decreased to ξ1=ξ−δξ<ξ and the granular temperature approaches the corresponding steady-state value T1= Ts−δT<Ts. At the waiting time, the granular temperature equals its steady value for the driving ξ. However, for longer times the granular temperature does not remain flat: it displays a nonmonotonic behavior, the Kovacs hump. The normal case is depicted, in which the temperature displays a maximum, but anomalous behavior—a minimum instead of a maximum—can also be observed. Note that δTand δTare defined in such a way that they are both positive. δT—sometimes termed the cooling protocol. The Kovacs effect is brought to the fore if T(τ) does not remain flat for τ>τw, signaling that, actually, the system is not in the NESS for the driving intensity ξat τ=τw.IfT(τ)−Ts>0, the Kovacs effect is normal—this is the situation that is always found in molecular systems [40,41,47], and the one plotted in the figure. If T(τ)−Ts<0, the Kovacs effect is anomalous, a possibility that has recently been reported in some athermal systems [16,18,53]. In the following, we analyze the Kovacs hump in linear response. Note that the theoretical approach presented here essentially differs from that in Refs. [16,18], where the nonlinear case—large jumps in the driving—was investigated. Therein, a perturbative analysis in powers of a2,smade it possible to derive an analytical expression for the Kovacs hump, but only when the system freely cooled in the waiting time window [65]. The linear response approach developed below, although limited to small jumps in the driving, improves our understanding of the Kovacs effect by connecting the shape of the hump with the properties of the direct relaxation function. The evolution of the temperature for τ⩾τwdirectly follows from the general linear response scheme in Ref. [38]. Namely, it is given by δT(τ)=(T0−T1)[φθ(τ)−φθ(τw)φθ(τ−τw)],(62) where φθis the direct relaxation function (i.e., for the one jump experiment), normalized in the sense that φθ(τ=0) = 024903-9