Large Mpemba-like effect in a gas of inelastic rough hard spheres

TL;DR

This paper demonstrates a giant Mpemba effect in a gas of inelastic rough spheres, where hotter samples can cool faster than colder ones, explained by coupled rotational and translational temperature dynamics.

Contribution

It reveals a large Mpemba effect in a gas of inelastic rough spheres and explains it using a simple Maxwellian approximation for velocity distribution.

Findings

The effect occurs even with temperature differences over an order of magnitude.

Theoretical predictions match molecular dynamics and Monte Carlo simulations.

Rotational and translational temperatures are key to the effect.

Abstract

We report the emergence of a giant Mpemba effect in the uniformly heated gas of inelastic rough hard spheres: The initially hotter sample may cool sooner than the colder one, even when the initial temperatures differ by more than one order of magnitude. In order to understand this behavior, it suffices to consider the simplest Maxwellian approximation for the velocity distribution in a kinetic approach. The largeness of the effect stems from the fact that the rotational and translational temperatures, which obey two coupled evolution equations, are comparable. Our theoretical predictions agree very well with molecular dynamics and direct simulation Monte Carlo data.

Equations

$$\theta (t)=\frac{T_r(t)}{T_t(t)},T(t)=\frac{T_t(t)+T_r(t)}{2}.$$

$${\theta }^{\text{st}}=$$

$${\gamma }^{\text{st}}\equiv$$

$${\partial }_{t^{\ast }}\ln T^{\ast }=\Phi (T^{\ast },\theta ),{\partial }_{t^{\ast }}\ln \theta =\Psi (T^{\ast },\theta ),$$

$$\Phi (T^{\ast },\theta )=$$

$$\Psi (T^{\ast },\theta )=$$

$${\Phi }_1$$

$${\Phi }_3$$

$$\Phi (T_{0A}^{\ast },{\theta }_{0A})<\Phi (T_{0B}^{\ast },{\theta }_{0B}).$$

$$\tilde{\theta }(\alpha ,\beta )=2-3\kappa -3(1+\kappa )\frac{1-{\alpha }^2}{1-{\beta }^2}$$

Full text

Let us consider two beakers of water at different temperatures. Mpemba and Osborne showed that the initially hotter sample cools sooner under certain conditions Mpemba and Osborne (1969), i.e., the curve giving the time evolution of its temperature crosses that of the initially cooler sample and stays below it for longer times. This is called the Mpemba memory effect, which is known since antiquity in cultures for which water in the form of ice and snow is common Ross (1931). Later, the Mpemba effect has been clearly identified in different physical systems Greaney et al. (2011); Ahn et al. (2016); Lu and Raz (2017); Lasanta et al. (2017); Baity-Jesi et al. (2018), although there is still some debate about its existence in water Burridge and Linden (2016).

From a physical point of view, one would like to answer how different the initial preparation of two samples of the system under study must be so that the Mpemba effect arises. This is the main—currently unresolved in general—question, although there has been some recent progress in this respect Lu and Raz (2017); Lasanta et al. (2017). Lu and Raz Lu and Raz (2017) analyzed the Mpemba effect in a generic Markovian system by monitoring the relaxation of an entropy-like variable that measures the distance to the steady state. This makes it possible to define and investigate Mpemba-like effects in systems for which there is not an obvious definition of a nonequilibrium temperature, but makes the comparison with the usual experimental setup described above difficult.

A different approach was carried out by some of us in the study of the Mpemba effect for a granular fluid of smooth hard spheres Lasanta et al. (2017). Therein, the granular temperature—basically the average kinetic energy per particle—is the physical quantity monitored to investigate the Mpemba effect. In the smooth-sphere case, the angular velocities play no role since there is no energy transfer between the translational and the rotational degrees of freedom, and the kinetic energy is thus purely translational. We showed that the Mpemba effect stems from the coupling of the granular temperature and the kurtosis, which measures the deviation of the velocity distribution function from the Maxwellian shape at the lowest order. More specifically, it is the difference between the initial values of the kurtosis of the two samples that controls the appearance of the Mpemba effect.

In the granular fluid of smooth hard spheres, the kurtosis is typically small. On the one hand, this facilitates the theoretical analysis, because it makes it possible to linearize the evolution equations and thus give a quantitative prediction of how different the initial kurtoses must be to enable the Mpemba effect. On the other hand, the smallness of the kurtosis limits the magnitude of the Mpemba effect: The initial temperatures must be quite close; see Eq. (5) and Fig. 1(b) of Ref. Lasanta et al. (2017).

It has very recently been shown that a different memory effect, the Kovacs effect Kovacs et al. (1979); Bertin et al. (2003); Mossa and Sciortino (2004); Aquino et al. (2006); Prados and Brey (2010), is much larger and more complex in a granular gas of rough spheres Lasanta et al. (2019) than in the smooth-sphere case Prados and Trizac (2014); Trizac and Prados (2014). What makes it possible to understand the largeness of the Kovacs effect is the coupling between the translational and rotational temperatures, which are of the same order of magnitude. In addition, the basic features of the memory effect can be understood within the Maxwellian (Gaussian) approximation, without having to resort to higher order cumulants.

The above picture prompts us to look into the Mpemba effect in a fluid of rough inelastic hard spheres. Remarkably, we show in the following that the Mpemba effect can be explained within a Gaussian framework, both qualitatively and quantitatively. Physically speaking, the coupling between the rotational and translational degrees of freedom and thus the existence of two comparable but different temperatures suffices to explain the memory effect. Moreover, we give a picture of the underlying physical conditions and discuss the possible relevance of the two-temperature mechanism in other systems.

Therefore, let us consider a dilute gas of inelastic rough hard spheres, with mass $m$, diameter $\sigma$, and moment of inertia $I$. Henceforth, we employ the dimensionless moment of inertia $\kappa\equiv 4I/m\sigma^{2}$, and its specific value for uniform solid spheres ($\kappa=\frac{2}{5}$) whenever a definite value is needed. Translational and angular particle velocities are denoted as $\bm{v}$ and $\bm{\omega}$, respectively.

Collisions between macroscopic particles are inelastic, i.e., energy is not conserved Foerster et al. (1993). For the inelastic rough hard sphere model, a binary collision is characterized by two parameters: the coefficient of normal restitution $0\leq\alpha\leq 1$ and the coefficient of tangential restitution $-1\leq\beta\leq 1$ Kremer (2010); Brilliantov and Pöschel (2004); Garzó (2019). They are intrinsic properties of the material and determine the shrinking of the normal and tangential components of the relative velocity of the two surface points at contact not (a). The collisional model based on these two parameters is sufficiently accurate in a variety of materials Foerster et al. (1993), with most of them presenting experimental values in the intervals $\alpha\in(0.7,0.95)$ and $\beta\in(-0.5,0.5)$ Louge . Consistently, the analysis carried out in this paper focuses on this region of the $(\alpha,\beta)$ plane.

Additionally, energy is homogeneously injected to the translational degrees of freedom of all the grains by a stochastic thermostat $\mathbf{F}$ modeled as a Gaussian white noise, i.e., $\langle{\bf F}_{i}(t)\rangle={\bf 0}$, $\langle{\bf F}_{i}(t){\bf F}_{j}(t^{\prime})\rangle=\mathsf{I}m^{2}\chi_{0}^{2}\delta_{ij}\delta(t-t^{\prime})$, where $i,j$ refer to the particles’ indexes, $\mathsf{I}$ is the $3\times 3$ unit matrix, and $\chi_{0}^{2}$ gives the “strength” of the stochastic forcing Williams (1996); Williams and MacKintosh (1996); Swift et al. (1998); van Noije and Ernst (1998); Montanero and Santos (2000); Gradenigo et al. (2011).

Here, we provide the minimal theoretical framework needed for the understanding of the Mpemba effect in the granular gas (see Refs. Vega Reyes et al. (2014); Vega Reyes and Santos (2015); sup for a detailed account of the kinetic theory calculations). The dynamics of our system is governed by the inelastic Boltzmann–Fokker–Planck equation for the single-particle velocity distribution function $f(\bm{v},\bm{\omega},t)$ Montanero and Santos (2000). From the kinetic equation, the evolution equations for the average quantities of interest are derived Vega Reyes et al. (2014); Vega Reyes and Santos (2015).

We restrict ourselves to homogeneous and isotropic states, for which $\langle\bm{v}\rangle=\mathbf{0}$ and $\langle\bm{\omega}\rangle=\mathbf{0}$. The basic physical information is thus encoded in the translational and rotational granular temperatures $T_{t}={\frac{m}{3}}\langle v^{2}\rangle$ and $T_{r}={\frac{I}{3}}\langle\omega^{2}\rangle$. Alternatively, the same information is provided by the temperature ratio $\theta$ and the total temperature $T$,

[TABLE]

The granular gas is inherently a nonequilibrium system and, therefore, equipartition is broken, i.e., $\theta\neq 1$. Thus, the simplest description of the rough-sphere granular gas is provided by the Maxwellian approximation, in which the following bivariate Gaussian form is assumed for the velocity distribution function, $f(\bm{v},\bm{\omega},t)\simeq n\!\left[{mI}/{4\pi^{2}T_{t}(t)T_{r}(t)}\right]^{\frac{3}{2}}\exp\left[-\frac{mv^{2}}{2T_{t}(t)}-\frac{I\omega^{2}}{2T_{r}(t)}\right]$, where $n$ is the number density.

In the long-time limit, the granular gas reaches a steady state due to the action of the stochastic force. This steady state is completely characterized by $\theta^{\text{st}}$ and $T^{\text{st}}$ in the Maxwellian approximation. Their expressions in terms of the coefficients of restitution and the stochastic forcing intensity are Vega Reyes and Santos (2015)

[TABLE]

Note that $\theta^{\text{st}}\leq 1$ is independent of $\alpha$ in the Maxwellian approximation but higher order approximations introduce a—rather weak—dependence on $\alpha$ Vega Reyes and Santos (2015).

It is useful to introduce dimensionless variables for temperature and time. Then, we define $T^{*}\equiv T/T^{\text{st}}$ and $t^{*}\equiv 2n\sigma^{2}\sqrt{\pi T_{t}^{\text{st}}/m}t$. In the Maxwellian approximation, $T^{*}$ and $\theta$ evolve according to sup

[TABLE]

with the definitions

[TABLE]

Above, we have introduced the notation

[TABLE]

with $K\equiv\kappa(1+\beta)^{2}/(1+\kappa)^{2}$ not (b). Note that the time evolution of the temperature is governed by the function $\Phi$, which does not only depend on $T^{*}$; this is a necessary condition for the appearance of the Mpemba effect.

Imagine two initial states $(T_{0A}^{*},\theta_{0A})$ and $(T_{0B}^{*},\theta_{0B})$, with $T_{0A}^{*}>T_{0B}^{*}>1$, of the same granular gas, i.e., with the same values of the coefficients of restitution $\alpha$ and $\beta$. Let us denote by $T_{A}^{*}(t^{*})$ and $T_{B}^{*}(t^{*})$ the associated decays of the temperature to the steady state: A Mpemba-like effect is brought about when there exists a crossing time $t^{*}_{\times}$ such that $T_{A}^{*}(t^{*})<T_{B}^{*}(t^{*})$ for $t^{*}>t^{*}_{\times}$.

A necessary—and physically intuitive—condition for having the Mpemba effect is that the initially hotter sample cools faster than the cooler one for short times, when the system still keeps “memory” of its initial conditions and is in the first stage of the so-called kinetic regime. For short enough times, we can consider that the system is exponentially cooling with a characteristic rate roughly equal to the initial value of $-\Phi$ not (c) and then a necessary condition for the Mpemba effect to be present is

[TABLE]

Let us investigate the behavior of $\Phi(T^{*},\theta)$ as a function of $\theta$, for fixed $T^{*}$, to understand under which conditions the Mpemba effect is expected. There are three distinct terms in $\Phi$: (i) the first one, $\Phi_{1}(T^{*})$, is a heating term that stems from the stochastic forcing and is thus independent of $\theta$, (ii) the second one, $\Phi_{2}(T^{*},\theta)$, is the typical cooling term of granular gases, which is also present for smooth spheres not (d), and (iii) the third one, $\Phi_{3}(T^{*},\theta)$, is a purely roughness term (note that $K=0$ for $\beta=-1$) and heats (cools) the system when $\theta<\theta^{\text{st}}$ ($\theta>\theta^{\text{st}}$). The sign and magnitude of $\Phi(T^{*},\theta)$ results from the competition among those three terms.

In light of the above, we analyze the behavior of $\Phi_{2}$ and $\Phi_{3}$ as a function of $\theta$, for fixed $T^{*}$. While the cooling term $\Phi_{2}$ is a monotonically increasing function of $\theta$ not (e), the roughness term $\Phi_{3}$ shows a more complex behavior. Starting from $\theta=0^{+}$, $\Phi_{3}$ first decreases with increasing $\theta$, vanishes at $\theta=\theta^{\text{st}}$, reaches a (negative) minimum value at $\theta=2+3\theta^{\text{st}}$, and finally tends to zero from below in the limit $\theta\to\infty$.

The overall behavior of $\Phi_{2}+\Phi_{3}$ as a function of $\theta$ depends on the values of the coefficients of restitution $(\alpha,\beta)$. Both $\Phi_{2}$ and $\Phi_{3}$ grow with increasing $\theta$ beyond the minimum of $\Phi_{3}$, i.e., for $\theta>2+3\theta^{\text{st}}$. On the other hand, as $\alpha$ approaches unity, the decay of $\Phi_{3}$ for small $\theta$ dominates over the growth of $\Phi_{2}$, resulting in a nonmonotonic dependence of $\Phi$ on $\theta$. This is illustrated in Fig. 1(a), which puts forward a density plot of $\Phi$ as a function of $(T^{*},\theta)$ for uniform solid spheres in the limiting case $(\alpha=1,\beta=0)$. A nonmonotonic behavior is neatly observed, especially to the right of the locus $\Phi(T^{*},\theta)=0$, i.e., where $\Phi<0$ and the system cools. As $\alpha$ is decreased, the magnitude of the cooling term $\Phi_{2}$ increases, eventually becoming the dominant one for small $\theta$ if $1-\alpha$ is large enough. Therein, a monotonically increasing behavior is observed, as illustrated in Fig. 1(b) for $(\alpha=0.9,\beta=0)$.

To carry out a more quantitative analysis, we impose that $\partial\Phi(T^{*},\theta)/\partial\theta|_{\widetilde{\theta}}=0$, which leads to

[TABLE]

and study the sign of $\widetilde{\theta}$ not (f). If $\widetilde{\theta}\leq 0$, there is no physically meaningful minimum and $\Phi$ is a monotonically increasing function of $\theta$. If $\widetilde{\theta}>0$, $\Phi$ displays a minimum at $\theta=\widetilde{\theta}$. Equation (7) implies that $\widetilde{\theta}>0$ if $\alpha$ is sufficiently close to unity and $|\beta|$ is sufficiently small. This is consistent with the qualitative discussion above, and it is illustrated in Fig. 2. Therein, the locus $\widetilde{\theta}=0$ separates the regions inside which $\Phi$ is monotonic (below it) and nonmonotonic (above it). On the one hand, $\widetilde{\theta}\leq 0$ for all $\beta$ when $\alpha\leq\alpha_{c}=\sqrt{(1+6\kappa)/3(1+\kappa)}$, which gives $\alpha_{c}=\sqrt{17/21}\simeq 0.9$ for uniform solid spheres. On the other hand, $\widetilde{\theta}>0$ for $\alpha>\alpha_{c}$ only if $\beta^{2}<1-3(1-\alpha^{2})(1+\kappa)/(2-3\kappa)$ not (g).

For the sake of conciseness, we restrict ourselves to the simpler monotonic situation $\alpha\leq\alpha_{c}$ in the remainder of the paper. Let us look again at Fig. 1(b), in which the limiting—less favorable—case $(\alpha=\alpha_{c},\beta=0)$ is shown. As already stated before, we consider two points, $A\equiv(T_{0A}^{*},\theta_{0A})$ and $B\equiv(T_{0B}^{*},\theta_{0B})$ with $T_{0A}^{*}>T_{0B}^{*}$, corresponding to different initial conditions, to study the Mpemba effect. Since $\Phi$ is monotonic, it suffices to take $\theta_{0A}<\theta_{0B}$ (i.e., the initially hotter system has its kinetic energy more concentrated in the translational modes than the initially cooler one) to fulfill Eq. (6). Equation (6), however, is not a sufficient condition for the Mpemba effect to appear: The disparity between $\Phi(T_{0A}^{*},\theta_{0A})$ and $\Phi(T_{0B}^{*},\theta_{0B})$ must be large enough, because the relaxation is not purely exponential. This entails that the disparity between $\theta_{0A}$ and $\theta_{0B}$ must be also large enough. This is illustrated in Fig. 3(a), where we analyze the emergence of the Mpemba effect for several pairs of the initial temperatures $(T_{0A}^{*},T_{0B}^{*})$. For each one of these pairs, the lines in the $(\theta_{0B},\theta_{0A})$ plane delimiting the regions with and without the Mpemba effect are plotted.

For the emergence of the Mpemba effect, the most favorable situation is having the kinetic energy completely concentrated in (i) the translational degrees of freedom for the higher temperature $T_{0A}^{*}$, i.e., $\theta_{0A}=0$, and (ii) the rotational degrees of freedom for the lower temperature $T_{0B}^{*}$, i.e., $\theta_{0B}\to\infty$. This limiting case is considered in Fig. 3(b). The range of initial temperature ratios $T_{0A}^{*}/T_{0B}^{*}$ leading to the emergence of the Mpemba effect increases (almost exponentially) with the colder initial temperature $T_{0B}^{*}$. For instance, if $T_{0B}^{*}=15.3$, $T_{0A}^{*}$ can be $100$ times larger than $T_{0B}^{*}$ and still the Mpemba effect is observed. In all the cases plotted in Fig. 3, we have considered that the Mpemba effect is present when the crossing time $t^{*}_{\times}$—if it ever exists—is smaller than $10$ not (h).

To check the accuracy of our theoretical predictions, we have carried out molecular dynamics (MD) and direct Monte Carlo (DSMC) simulations sup . Different systems are considered by varying the values of $\alpha$ and $\beta$. For each particular pair $(\alpha,\beta)$, the system is initialized at three points $A$, $B$, and $C$ with $(T^{*}_{0A},\theta_{0A})=(5,0.01)$, $(T^{*}_{0B},\theta_{0B})=(4,10)$, and $(T^{*}_{0C},\theta_{0C})=(3,100)$, respectively. According to our theoretical framework, the Mpemba effect is expected to emerge in the cases $(A,B)$ and $(A,C)$, but not in the case $(B,C)$.

Figure 4 shows the time evolution of the granular temperature. The curves correspond to $\alpha=0.7$ and $0.9$, and $\beta=-0.5$, [math], and $0.5$, but the qualitative picture explained below remains unaltered for other choices, as long as $\alpha\leq\alpha_{c}\simeq 0.9$. The Mpemba effect is clearly observed as the initially hottest sample ($A$) is the fastest one to reach the steady state. Moreover, in all the cases, the three plotted curves—the MD data (circles), the DSMC data (triangles), and the Maxwellian theory (line)—present a remarkably excellent agreement.

Qualitatively, the strength of the Mpemba effect—the separation between the two relaxation curves for times longer than $t_{\times}^{*}$—increases as $\alpha$ decreases at fixed $\beta$ (higher inelasticity) and as $\beta$ goes to smaller values at fixed $\alpha$ (lower roughness). Thus, quantitatively, we can measure the magnitude of the Mpemba effect by defining a Mpemba parameter $\text{Mp}_{\!AB}$ as the extremum of the difference of dimensionless temperatures $T_{B}^{*}(t)-T_{A}^{*}(t)$ for $t^{*}>t_{\times}^{*}$. The values of $\text{Mp}_{\!AB}$ and $\text{Mp}_{\!AC}$ obtained from the solution of Eq. (3) are shown in each one of the panels of Fig. 4 and are consistent with the qualitative behavior described above. For fixed roughness $\beta$, the more inelastic the system is, the larger Mp becomes. For fixed inelasticity $\alpha$, the smoother the system is, the larger Mp becomes. Note that the values of Mp are typically of the order of unity for our system of rough spheres. In the smooth-sphere case, however, Mp is much smaller, being typically of the order of a few thousandths for the cases reported in Ref. Lasanta et al. (2017). While for smooth spheres the influence of the kurtosis on the evolution of $T^{*}$ is very weak, a strong impact of the rotational-translational partition is present in the case of rough spheres.

We have neatly observed the Mpemba effect in the gas of inelastic rough hard spheres. It stems from the coupling in the evolution of the rotational and translational temperatures, which are, in general, of the same order, $\theta=O(1)$. For low enough $\alpha$, i.e., large enough inelasticity, we have shown that the more concentrated in the translational degrees of freedom the total kinetic energy is, the faster the system cools. Therefore, the initially hotter system must have its kinetic energy more concentrated in the translational modes than the initially cooler one to facilitate the Mpemba effect. Moreover, we have quantified how large this concentration must be. In the rough-sphere case, the initial temperatures of the two samples can be quite different and definitely do not need to be close to each other, as happens in the smooth-sphere granular gas.

The strength of the effect has been quantified by the Mpemba parameter Mp defined as the maximum separation—relative to the stationary temperature—between the relaxation curves after the crossing. The Mpemba effect in the rough-sphere granular gas is really huge, with typical order of unity values of Mp. This has to be contrasted with the Mpemba effect in the smooth-sphere case that, although distinctly observed in Ref. Lasanta et al. (2017), is quite small ($\text{Mp}\sim 10^{-3}$).

Also, it is worth emphasizing that the Mpemba effect can be explained with a minimal Gaussian approximation for the distribution function. The largeness of the effect shows that the system is sweeping far-from-equilibrium states but, interestingly, these states can be described by a sort of “local equilibrium” approximation but with two distinct temperatures. This is at difference with the hydrodynamic regime, in which the rotational temperature becomes enslaved to the total one for long enough times, i.e., $\theta$ becomes time independent over the hydrodynamic scale and $T$ is the relevant hydrodynamic field.

In the rough-sphere case, the corrections introduced by the cumulants are expected to remain small not (i). This expectation is supported by the exceptionally good agreement found in Fig. 4 between the Maxwellian analytical predictions and the DSMC data—which give the numerical integration of the kinetic equation. In addition, the agreement of our theory with MD simulations tells us that our kinetic description holds for the low-density gas.

The results in the present work suggest that the two-temperature mechanism found here may be significant for observation of a huge—and thus experimentally measurable—Mpemba effect in a variety of systems with several distinct temperatures, such as structural glasses Narayanaswami (1971); DeBolt et al. (1976), metals undergoing plastic deformation R. Chowdhury and Roy (2019), and granular mixtures Dahl et al. (2002). This is enough to bring to the fore the Mpemba effect, which may even be enhanced if other variables like higher order cumulants are also relevant.

This work has been supported by the Spanish Ministerio de Ciencia, Innovación y Universidades and the Agencia Estatal de Investigación Grants (partially financed by the ERDF) No. MTM2017-84446-C2-2-R (A.L. and A.T.), No. FIS2016-76359-P (F.V.R. and A.S.), No. FIS2017-84440-C2-2-P (A.T.), and No. PGC2018-093998-B-I00 (A.P.), and by the Junta de Extremadura (Spain) Grants No. IB16087 and No. GR18079 (F.V.R., A.S., and M.L.C.), partially funded by the ERDF. Use of computing facilities from Extremadura Research Center for Advanced Technologies (CETA-CIEMAT), funded by the ERDF, is also acknowledged.