When the Hotter Cools More Quickly: Mpemba Effect in Granular Fluids

Under certain conditions, two samples of fluid at different initial temperatures present a counterintuitive behavior known as the Mpemba effect: it is the hotter system that cools sooner. Here, we show that the Mpemba effect is present in granular fluids, both in uniformly heated and in freely cooling systems. In both cases, the system remains homogeneous, and no phase transition is present. Analytical quantitative predictions are given for how differently the system must be initially prepared to observe the Mpemba effect, the theoretical predictions being confirmed by both molecular dynamics and Monte Carlo simulations. Possible implications of our analysis for other systems are also discussed.

Let us consider two identical beakers of water, initially at two different temperatures, put in contact with a thermal reservoir at subzero (on the Celsius scale) temperature.While one may intuitively expect that the initially cooler sample would freeze first, it has been observed that this is not always the case [1].This paradoxical behavior named the Mpemba effect (ME) has been known since antiquity and discussed by philosophers like Aristotle, Roger Bacon, Francis Bacon, and Descartes [2,3].Nevertheless, physicists only started to analyze it in the second part of the past century, mainly in popular science or education journals .
There is no consensus on the underlying physical mechanisms that bring about the ME.Specifically, water evaporation [4,5,9,24], differences in the gas composition of water [11,17,25], natural convection [6,23,26], or the influence of supercooling, either alone [14,27] or combined with other causes [28][29][30][31], have been claimed to have an impact on the ME.Conversely, the own existence of the ME in water has been recently put in question [32].Notwithstanding, Mpemba-like effects have also been observed in different physical systems, such as carbon nanotube resonators [33] or clathrate hydrates [34].
The ME requires the evolution equation for the temperature to involve other variables, which may facilitate or hinder the temperature relaxation rate.The initial values of those additional variables depend on the way the system has been prepared, i.e., "aged," before starting the relaxation process.Typically, aging and memory effects are associated with slowly evolving systems with a complex energy landscape, such as glassy [35][36][37][38][39][40][41][42][43] or dense granular systems [44][45][46].However, these effects have also been observed in simpler systems, like granular gases [47][48][49][50] or, very recently, crumpled thin sheets and elastic foams [51].
In a general physical system, the study of the ME implies finding those additional variables that control the temperature relaxation and determining how different they have to be initially in order to facilitate its emergence.In order to quantify the effect with the tools of nonequilibrium statistical mechanics, a precise definition thereof is mandatory.An option is to look at the relaxation time to the final temperature as a function of the initial temperature [1,2,4,5,9,11,24,32,34].Alternatively, one can analyze the relaxation curves of the temperature: if the curve for the initially hotter system crosses that of the initially cooler one and remains below it for longer times, the ME is present [3,9,18,20,[27][28][29][30]33].
In this Letter, we combine both alternatives above and investigate the ME in a prototypical case of intrinsically out-of-equilibrium system: a granular fluid [52][53][54][55], i.e., a (dilute or moderately dense) set of mesoscopic particles that do not preserve energy upon collision.As a consequence, the mean kinetic energy, or granular temperature T (t) [53], decays monotonically in time unless an external energy input is applied.The simplicity of the granular fluid makes it an ideal benchmark for other, more complex, nonequilibrium systems.We analyze the time evolution of the granular fluid starting from different initial preparations and quantitatively investigate how the ME appears.This is done for both the homogeneous heated and freely cooling cases.
Our granular fluid is composed of smooth inelastic hard spheres.Therefore, the component of the relative velocity along the line joining the centers of the two colliding particles is reversed and shrunk by a constant factor α [52], the so-called coefficient of normal restitution.In addition, the particles are assumed to be subject to random forces in the form of a white-noise thermostat with variance m 2 ξ 2 , where m is the mass of a particle.Thus, the velocity distribution function (VDF) f (v, t) obeys an Enskog-Fokker-Planck kinetic equation [56][57][58].
The granular temperature and the excess kurtosis (or second Sonine coefficient) are defined as respectively, where n = dvf (v, t) is the number density.From the kinetic equation for the VDF, one readily finds [56] where κ ≡ 2ng(n)σ 2 π/m, σ and g(n) are the sphere diameter and the pair correlation function at contact [59], respectively, χ ≡ 3m 2κ ξ 2 , and µ 2 and µ 4 are dimensionless collisional rates.
Note that Eqs.(1) are formally exact, but (a) T and a 2 are coupled, and (b) the equations are not closed in those two variables since µ n are functionals of the whole VDF.However, if inelasticity is not too large, the nonlinear contributions of a 2 and the complete contributions of higher order cumulants can be neglected.This is the so-called first Sonine approximation [56,60], which yields µ n ≃ µ . Using the first Sonine approximation above in Eqs. ( 1), they become a closed set, but the T -a 2 coupling still remains.Taking into account this coupling, and since µ 2 is an increasing function of a 2 , it turns out that the relaxation of the granular temperature T from an initially "cooler" (smaller T ) sample could possibly be overtaken by that of an initially "hotter" one, if the initial excess kurtosis of the latter is larger enough.We build on and quantify the implications of this physical idea in the following.
First, we consider the uniformly heated case (i.e., χ = 0) and prepare the granular fluid in an initial state that is close to the steady one, in the sense that Eqs.(1) can be linearized around the stationary values [56,57] 4 − 5µ 2 ], where where the matrix Λ has elements 2 (1+a st 2 )].Thus, the relaxation of the temperature reads where Let us imagine two initial states f 0 = f A and f B , with (θ 0 , a 2,0 ) = (θ A , a 2A ) and (θ B , a 2B ), respectively.We assume that θ A > θ B and a 2A > a 2B .Both cooling (θ A > θ B > 1) and heating (θ B < θ A < 1) processes may be considered.From Eq. ( 3), the time τ c for the possible crossing of the two relaxation curves satisfies where ∆θ 0 = θ A − θ B and ∆a 2,0 = a 2A − a 2B .For a given α, in this simplified description the crossover time τ c depends on (θ A , a 2A ) and (θ B , a 2B ) (or, more generally, on the details of the two initial VDFs f A and f B ) only through the single control parameter ∆θ 0 /∆a 2,0 .Figure 1(a) displays τ c as a function of the ratio ∆θ 0 /∆a 2,0 for α = 0.9.Equation (4) implies that there is a maximum of the control parameter ∆θ 0 /∆a 2,0 for which the ME can be observed, namely, This quantity determines the phase diagram for the occurrence of the ME, as shown in Fig. 1(b).Equation ( 5) can be read in two alternative ways.First, it means that, for a given difference ∆a 2,0 of the initial kurtosis, the ME appears when the difference ∆θ 0 of the scaled initial temperatures is below a maximum value (∆θ 0 ) max , proportional to ∆a 2,0 .Second, for a given value of ∆θ 0 , the ME is observed only for a large enough difference of the initial kurtosis, i.e., ∆a 2,0 > (∆a 2,0 ) min , with (∆a 2,0 ) min proportional to ∆θ 0 .This quantitatively measures how different the initial conditions of the system must be in order to have the ME.In order to check the accuracy of our theoretical results, we compare them in Fig. 2 with MD simulations (at a density nσ 3 = 0.02) and with the numerical integration of the Enskog-Fokker-Planck equation by means of the DSMC method [61].In all our simulations, α = 0.9 and the initial VDF is assumed to have a gamma-distribution form [62] in the variable v 2 with parameters adjusted to reproduce the chosen values of θ 0 and a 2,0 .First, three different initial conditions (A, B, and C ) with temperatures above the stationary, θ A = 1.04, θ B = 1.035, and θ C = 1.03, and excess kurtosis a 2A = 0.5, a 2B = 0, and a 2C = −0.35,are considered.The ME is clearly observed as a crossover of the relaxation curves of the temperature [see, also, Fig. 1(a)].Second, we analyze a "heating" protocol by choosing initial temperatures below the steady value, namely, θ ′ A = 0.97, θ ′ B = 0.965, and θ ′ C = 0.96, with the same values of the excess kurtosis as in the "cooling" case.Again, a crossover in the temperature relaxation curves appears, signaling the granular analog of the inverse ME proposed in a recent work [63].It is interesting to note that the evolution curves corresponding to θ C = 1.03 and θ ′ A = 0.97 are nonmonotonic.This peculiar behavior is predicted by Eq. (3) to take place if 2 /µ st 2 < δθ 0 /δa 2,0 < 0. Alternatively, we can characterize the system celerity for cooling (or heating) by defining a relaxation time τ ǫ as the time after which |θ(τ ǫ ) − 1| < ǫ, with ǫ ≪ 1.From Eq. (3), Figure 3 shows τ ǫ as a function of the initial temperature θ 0 for ǫ = 10 −4 and the same values of the initial excess kurtosis as considered in Fig. 2. In this diagram, for a given pair of a 2,0 , the range of initial temperatures for which the ME emerges is clearly visualized.Note that this range does not change if the value of the bound ǫ is changed to ǫ ′ , since the diagram is only shifted vertically by an amount 1 λ− ln(ǫ/ǫ ′ ).A relevant question is whether or not the ME keeps appearing in the zero driving limit.In the undriven case (χ = 0), the granular fluid relaxes to the so-called homogeneous cooling state (HCS) [52], which is the reference state for deriving the granular hydrodynamics [64].If the linear relaxation picture developed above remained valid in the nonlinear relaxation regime, at least qualitatively, the answer would be negative.Note that the maximum temperature difference (∆T 0 ) max would vanish in the limit as χ → 0 (T st ∝ χ 2/3 → 0), as a consequence of (∆θ 0 ) max being independent of χ.Interestingly, we show below that this simple scenario does not hold, and the ME is also observed for very small driving: indeed, (∆T 0 ) max remains finite in this limit.
For very small driving, there is a wide initial time region inside which the system evolves as if it were cooling freely.Therefore, for the sake of simplicity, we now take χ = 0 in the evolution equations (1).While the system freely cools for all times (lim t→∞ T = 0), the excess kurtosis tends to a constant value [56,57]  laxed to its stationary value a HCS

2
. So as to look for a possible crossover of the cooling curves, we linearize the equations around T * = 1 (by choosing T ref such that the initial temperatures verify |T * 0 − 1| ≪ 1) and a 2 = a HCS 2 .Therefrom, the evolution of T * is obtained as where δT * = T * − 1, δa 2,0 = a 2,0 − a HCS 2 , and λ a = Thus, the ME actually survives in the zero driving limit.Had we considered a small value of the driving χ instead of χ = 0, Eqs. ( 7)-( 9) would characterize the strongly nonlinear regime, in which the initial scaled temperature θ 0 = T 0 /T st ≫ 1.In a first stage of the relaxation, as long as the granular temperature T ≫ T s , the driving can be neglected, the system freely cools, and the ME is observed provided that the condition (9) is fulfilled.Afterwards, the initially hotter system remains below the initially cooler one forever.When approaching the steady state, both the temperature and the excess kurtosis start to evolve towards their stationary values T st and a st 2 , but in both curves one has a 2,0 = a HCS 2 , and Eq. ( 5) tells us that no further crossing of the curves takes place (∆a 2,0 = 0).
In summary, we have shown by means of a simple analytical theory that the ME naturally appears in granular fluids, as a consequence of the relevance of non-Gaussianities in the time evolution of T .Specifically, this allows us to (i) prove that the ME is to be expected on quite a general basis and for a wide range of systems, as long as non-Gaussianities are present and (ii) quantitatively predict the region of parameters within which the ME is present.Moreover, we have also predicted the existence of an inverse ME: when the system is heated instead of cooled, the initially cooler sample may heat sooner [63].In this way, we have provided a general theoretical framework for the understanding of the ME.
The main assumptions in our theory are: (i) the validity of the kinetic description, (ii) the system remaining homogeneous for all times, and (iii) the first Sonine approximation within the kinetic description.All these assumptions have been validated in the paper.First, the numerical integration of the Enskog equation provided by the DSMC simulations has been successfully compared with MD simulations.Second, we have also checked that the system remains homogeneous in the MD simulations, both for the heated and undriven cases.Concretely, in the latter, the system size has been chosen to be well below the clustering instability threshold [55,65].Third, the accuracy of the first Sonine approximation has been confirmed by the excellent agreement between our analytical results and the DSMC simulations, even for the not-so-small values of the excess kurtosis a 2 considered throughout.
Finally, we stress that non-Gaussianities may have a leading role in the emergence of the ME in other systems, even when there is no inelasticity.For example, the temperature of a molecular fluid, which is basically the mean kinetic energy per particle, does not remain constant if the system interacts with a thermal reservoir.Let us assume that the coupling with the reservoir brings about a nonlinear drag, as considered, for instance, in Refs.[66][67][68][69].Then, the evolution equation of the temperature would involve higher moments of the transient nonequilibrium VDF.In this quite general situation, the ME would also stem from those non-Gaussianities [70].

Figure 1 .
Figure 1.(a) Crossover time τc as a function of the ratio ∆θ0/∆a2,0 for α = 0.9.(b) Phase diagram in the plane ∆θ0/∆a2,0 vs α.The regions of the plane inside which there appears or does not appear the ME are separated by the curve (∆θ0/∆a2,0)max.

Figure 2 .
Figure 2. Relaxation of the scaled temperature to the steady state for α = 0.9.The upper and the lower curves correspond to the ME for the cooling and heating processes, respectively (see text).The direct simulation Monte Carlo (DSMC) (open symbols) and molecular dynamics (MD) (filled symbols) data show an excellent agreement with the theoretical prediction (lines).

Figure 4 .
Figure 4. Evolution of the temperature in the free cooling case.Again, the agreement between the theory (lines) and both DSMC (open symbols) and MD (filled symbols) simulation data is excellent.The inset shows t * c as a function of ∆T * 0 /∆a2,0.

Fig. 4 .
In turn, a 2 decays exponentially to a HCS 2 with a characteristic time λ −1 a .Similar to the thermostatted case, we consider two initial states (T * 0 , a 2,0 ) = (T * A , a 2A ) and (T * B , a 2B ), with∆T * 0 = T * A − T * B > 0, ∆a 2,0 = a 2A − a 2B > 0.Logically, only the cooling case makes sense.In Fig.4, we plot two relaxation curves of the temperature for α = 0.9, with T * A = 1, T * B = 0.99, a 2A = 0.5, a 2B = −0.35,with the choice T ref = T A .The ME is clearly observed, and the crossover time t * c is Therefore, there is a maximum value of the ratio ∆T * 0 /∆a 2.0 for which the ME appears,