Mpemba-like effect protocol for granular gases of inelastic and rough hard disks

1 Introduction

Since the Antiquity, the fact that water could start freezing earlier for initially hotter samples was observed and commented by very influential people of different epochs like Aristotle [1], Francis Bacon [2], or René Descartes [3]. This counterintuitive phenomenon contradicts Isaac Newton’s formulation of its well-known cooling’s law [4, 5], but otherwise it is part of the popular belief in cold countries. The scientific community started to pay attention to this effect since the late 60s of the last century thanks to its accidental rediscovery by a Tanzanian high-school student, Erasto B. Mpemba. Later, he and Dr. Denis Osborne reported their findings [6, 7] and since then the effect is usually known as Mpemba effect (ME).

Whereas the original tested system for ME has been water [6–37], it is still under discussion and no consensus about its occurrence has been agreed [38–40]. In fact, the statistical physics community is currently paying attention to Mpemba-like effects that have been described in a huge variety of complex systems in the last decades, such as ideal gases [41], molecular gases [42–44], gas mixtures [45], granular gases [46–52], inertial suspensions [53, 54], spin glasses [55], Ising models [56–58], non-Markovian mean-field systems [59, 60], carbon nanotube resonators [61], clathrate hydrates [62], active systems [63], or quantum systems [64]. The theoretical approach to the fundamentals of the problem has been done via different routes like Markovian statistics [65–69] or Landau’s theory of phase transitions [70]. Recently, in the context of a molecular gas under a nonlinear drag force, new interpretations and definitions of ME from thermal and entropic point of views, as well as a classification of the whole possible phenomenology, have been carried out [44]. In addition, ME has been experimentally observed in colloids [71, 72], proving that it is a real effect present in nature.

The very first time that ME was observed theoretically in granular gaseous rapid flows was in [46]. The considered system was a set of inelastic and smooth hard spheres (with constant coefficient of normal restitution) heated by a stochastic thermostat, the effect arising by initially preparing the system in far from Maxwellian states. The same type of initial preparation was applied to the case of molecular gases with nonlinear drag [42, 44]. Essentially, the temperature evolution depends on the whole moment hierarchy of the velocity distribution function (particularly on the excess kurtosis, or fourth cumulant, and, more weakly, on the sixth cumulant), this dependence giving rise to the possible appearance of ME.

On the other hand, there is no need to consider an initial velocity distribution function (VDF) far from the Maxwell–Boltzmann one if the temperature is coupled to other basic variables that can be fine-tuned in the initial preparation of the system. This occurs in the case of a monocomponent granular gas made by inelastic and rough hard spheres thermostatted by a stochastic force [47], as well as in driven binary mixtures of either molecular or inelastic gases [45, 48]. In those systems, one does not need to invoke strong nonGaussianities, since the temperature relaxation essentially depends on the rotational-to-translational temperature ratio (in the case of rough particles) or on the partial component temperatures (in the case of mixtures). A similar situation applies in the presence of anisotropy in either the injection of energy [51] or in the velocity flow [53]. However, there is still a lack of protocol defining a possible nearly realistic preparation of the initial states for a granular or molecular gas in homogeneous and isotropic states. Unlike other memory effects, such as Kovac’s effect [73, 74], ME has not a predefined way to elaborate a protocol.

In this work, we have addressed the latter preparation problem for a specific model of granular gases. We consider a monodisperse granular gas of inelastic and rough hard disks, where inelasticity is parameterized via a constant coefficient of normal restitution, α, and the roughness is accounted for by a coefficient of tangential restitution, β, assumed to be constant as well. Disks are “heated” by a stochastic thermostat which injects energy to both translational and rotational degrees of freedom through a combination of a stochastic force and a stochastic torque, both with properties of a white noise. The relative amount of energy injected to the rotational degree of freedom, relative to that injected to the translational degrees of freedom, can be freely chosen. Therefore, we will denote this thermostat as splitting thermostat (ST). The quantity coupled to the temperature that will monitor the possible occurrence of ME will be the rotational-to-translational temperature ratio, as in [47], where, however, a stochastic torque was absent. This double energy-injection based on ST allows us to fix the initial conditions of the variables that play a role in the evolution process, namely the temperature and its coupling. A side effect of providing energy to the rotational degree of freedom is that it favors the possibility of an overshoot of the temperature with respect to its steady-state value. This might cause ME, even in the absence of a crossing between the temperatures of the two samples [44]. Therefore, the protocol must be adapted to this specific phenomenon.

It is worth saying that our theoretical approach is based on a Maxwellian approximation (MA), that is, we assume that both transient and steady-state VDFs are close to a two-temperature Maxwellian. This approach is founded on previous works for the case of zero stochastic torque [75] and on preliminary results for the system at hand [76]. Moreover, the two-dimensional characterization of the physical system is thought to be plausible for hopefully being reproduced in some experimental setup. As will be seen, the reliability of our theoretical approach is confirmed by computer simulations via the direct simulation Monte Carlo (DSMC) method and event-driven molecular dynamics (EDMD).

The paper is structured as follows. In Section 2, the model system for a granular gas of inelastic and rough hard disks thermostatted by stochastic force and torque is introduced. Also, explicit evolution equations and expressions for the steady-state dynamic variables are shown under the MA, and the theoretical results are compared with DSMC and EDMD. Section 3 collects the definition and necessary conditions for ME to occur taking into account the emergence or not of overshoot during evolution. Subsequently, and based on the analysis of this section, two different protocols are presented for observing ME in cases without and with overshoot, respectively. This discussion is accompanied by its proper comparison with simulation results. Finally, concluding remarks are presented in Section 4.

2 The model

We consider a set of mechanically identical inelastic and rough hard disks of mass m, diameter σ, and reduced moment of inertia κ ≡ 4I/mσ² (I being the moment of inertia). The translational velocities lie on the xy plane, i.e., $\mathbf{v}=v_x\widehat{\mathbf{x}}+v_y\widehat{\mathbf{y}}$, while the angular velocities point along the orthogonal z axis, $\mathit{\omega }=\omega \widehat{\mathbf{z}}$. Inelasticity and roughness are characterized by the coefficients of normal and tangential restitution, α and β, respectively, which are assumed to be constant and defined as [77].

$\widehat{\mathit{\sigma }}\cdot {\mathbf{v}}_{12}^{\prime }=-\alpha \widehat{\mathit{\sigma }}\cdot {\mathbf{v}}_{12},{\widehat{\mathit{\sigma }}}_{\perp }\cdot {\mathbf{v}}_{12}^{\prime }=-\beta {\widehat{\mathit{\sigma }}}_{\perp }\cdot {\mathbf{v}}_{12},(1)$

Where $\widehat{\mathit{\sigma }}\equiv ({\mathbf{r}}_2-{\mathbf{r}}_1)/|{\mathbf{r}}_2-{\mathbf{r}}_1|$ is the unit intercenter vector along the line of centers from particle 1 to particle 2, ${\widehat{\mathit{\sigma }}}_{\perp }\equiv \widehat{\mathit{\sigma }}\times \widehat{\mathbf{z}}$ is orthogonal to $\widehat{\mathit{\sigma }}$, v₁₂ ≡v₁ − v₂ is the relative velocity between particles 1 and 2, and primed quantities account for their postcollisional values. Because of their definitions, the ranges of the coefficients of restitution are 0 ≤ α ≤ 1 and −1 ≤ β ≤ 1, α = 1 corresponding to elastic collisions, β = −1 describing a perfectly smooth disks, and β = 1 standing for completely rough disks. In fact, total kinetic energy is only conserved if α = |β| = 1 [77–80].

In the case α ≠ 1 or |β| ≠ 1, that is, when kinetic energy is dissipated upon collisions, the undriven system will evolve up to a completely frozen state. In order to avoid that quench, we will force the particles to externally receive energy via a homogeneous stochastic force F^wn and a homogeneous stochastic torque τ^wn that inject translational and rotational kinetic energy, respectively, with the properties of a white noise. That is,

$\langle {\mathbf{F}}_i^{\mathrm{w}\mathrm{n}}\left(t\right)\rangle =0,\langle {\mathbf{F}}_i^{\mathrm{w}\mathrm{n}}\left(t\right){\mathbf{F}}_j^{\mathrm{w}\mathrm{n}}\left(t^{\prime }\right)\rangle =\mathsf{I}{m}^2{\chi }_{\text{tr}}^2{\delta }_{ij}\delta \left(t-t^{\prime }\right),(\mathrm{2a})$

$\langle {\tau }_i^{\mathrm{w}\mathrm{n}}\left(t\right)\rangle =0,\langle {\tau }_i^{\mathrm{w}\mathrm{n}}\left(t\right){\tau }_j^{\mathrm{w}\mathrm{n}}\left(t^{\prime }\right)\rangle =mI{\chi }_{\text{rot}}^2{\delta }_{ij}\delta \left(t-t^{\prime }\right),(\mathrm{2b})$

where $\mathsf{I}$ is the 2 × 2 identity matrix, i and j are particle indices, and ${\chi }_{\text{tr}}^2$ and ${\chi }_{\text{rot}}^2$ are the intensities of the noises applied to the translational and rotational degrees of freedom, respectively. The combination of the stochastic force and torque is characterized by the pair of parameters $({\chi }_{\text{tr}}^2,{\chi }_{\text{rot}}^2)$ and defines the ST, as described in Section 1. In Section 2.1 we will introduce a more manageable pair of equivalent parameters. An illustration of the system is represented in Figure 1.

FIGURE 1

FIGURE 1. (A) Illustration of the system, where the green horizontal straight arrows and the blue curved ones represent the stochastic force F^wn and the stochastic torque τ^wn in Eqs. 2, respectively. (B) Mechanism of injection-dissipation of energy in the introduced model. Dissipative collisions produce a cooling effect on the translational and rotational temperatures (T_tr and T_rot), together with a transfer between the translational and rotational energies (equipartition effect). Additionally, the external white-noise force and torque inject energy (heating effect).

To dynamically describe the system, we will work under the assumptions of the homogeneous Boltzmann–Fokker–Planck equation (BFPE). That is, we consider a homogeneous and isotropic gas in a dilute regime, such that the evolution due to the collisional process is determined by just binary collisions, assuming Stosszahlansatz (or molecular chaos). The BFPE for this collisional model together, with the ST, is written as follows

$\frac{\partial }{\partial t}f\left(\mathbf{v},\omega ;t\right)-\frac{{\chi }_{\text{tr}}^2}{2}{\left(\frac{\partial }{\partial \mathbf{v}}\right)}^{2}f\left(\mathbf{v},\omega ;t\right)-\frac{{\chi }_{\text{rot}}^2}{2}{\left(\frac{\partial }{\partial \omega }\right)}^{2}f\left(\mathbf{v},\omega ;t\right)=J\left[\mathbf{v},\omega |f\right],(3)$

where f is the one-particle VDF, J[v, ω|f] is the usual Boltzmann collision operator for hard disks defined as [78–80].

$J\left[{\mathbf{v}}_1,{\omega }_1|f\right]=n\sigma \int \mathrm{d}{\mathbf{v}}_2{\int }_{-\infty }^{\infty }\mathrm{d}{\omega }_2{\int }_{+}\mathrm{d}\widehat{\mathit{\sigma }}|{\mathbf{v}}_{12}\cdot \widehat{\mathit{\sigma }}|\times \left[\frac{f\left({\mathbf{v}}_1^{\prime \prime },{\omega }_1^{\prime \prime };t\right)f\left({\mathbf{v}}_2^{\prime \prime },{\omega }_2^{\prime \prime };t\right)}{\alpha |\beta |}-f\left({\mathbf{v}}_1,{\omega }_1;t\right)f\left({\mathbf{v}}_2,{\omega }_2;t\right)\right],(4)$

where n is the number density, the subscript + in the integral over $\widehat{\mathit{\sigma }}$ means the constraint ${\mathbf{v}}_{12}\cdot \widehat{\mathit{\sigma }}>0$, and double primed quantities are precollisional velocities, which are given by [78–80].

${\mathbf{v}}_{1/2}^{″}={\mathbf{v}}_{1/2}\mp \left\{\frac{1+{\alpha }^{-1}}{2}\left({\mathbf{v}}_{12}\cdot \widehat{\mathit{\sigma }}\right)\widehat{\mathit{\sigma }}+\frac{\kappa }{1+\kappa }\frac{1+{\beta }^{-1}}{2}\left[{\mathbf{v}}_{12}\cdot {\widehat{\mathit{\sigma }}}_{\perp }-\frac{\sigma }{2}\left({\omega }_1+{\omega }_2\right)\right]{\widehat{\mathit{\sigma }}}_{\perp }\right\}.(5)$

2.1 Dynamics

The time evolution of the system is fully described by the BFPE, Eq. 3, which allows one to determine the dynamics of macroscopic quantities. The most important and basic quantities to study the dynamics of the system will be the translational and rotational granular temperatures defined at a certain time t as

$T_{\text{tr}}\left(t\right)=\frac{1}{2}m\langle v^2\rangle ,T_{\text{rot}}\left(t\right)=I\langle {\omega }^2\rangle ,(6)$

where the notation ⟨⋅⟩ means the average over the instantaneous VDF,

$⟨X\left(\mathbf{v},\omega ;t\right)⟩=n^{-1}\int \mathrm{d}\mathbf{v}{\int }_{-\infty }^{\infty }\mathrm{d}\omega X\left(\mathbf{v},\omega ;t\right)f\left(\mathbf{v},\omega ;t\right).(7)$

One should notice that the partial noise intensities ${\chi }_{\text{tr}}^2$ and ${\chi }_{\text{rot}}^2$ affect directly the evolution equations of T_tr(t) and T_rot(t), respectively. On the other hand, those quantities are coupled due to the transfer of energy during collisions (see Figure 1B). The description of the dynamics from T_tr(t) and T_rot(t) is equivalent to consider the mean granular temperature, T(t), and the rotational-to-translational temperature, θ(t), defined as follows

$T\left(t\right)=\frac{2}{3}{T}_{\text{tr}}\left(t\right)+\frac{1}{3}{T}_{\text{rot}}\left(t\right),\theta \left(t\right)=\frac{T_{\text{rot}}\left(t\right)}{T_{\text{tr}}\left(t\right)}.(8)$

In the definition of T(t) we have taken into account that there are two translational and one rotational degrees of freedom. This new pair of variables will be useful to study ME, which will be related to the evolution of the mean granular temperature T(t). In addition, this change of dynamical quantities in Eq. 8 induces a change of parameters describing the ST. Thus, we introduce the total noise intensity, χ², and the rotational-to-total noise intensity ratio, ɛ, as

${\chi }^2={\chi }_{\text{tr}}^2+\frac{I}{2m}{\chi }_{\text{rot}}^2,\epsilon =\frac{I}{2m}\frac{{\chi }_{\text{rot}}^2}{{\chi }^2}.(9)$

Notice that ${\chi }_{\text{tr}}^2=(1-\epsilon ){\chi }^2$. Therefore, 0 ≤ ɛ ≤ 1, ɛ = 0 and ɛ = 1 corresponding to the purely translational and the purely rotational thermostat, respectively. The total noise intensity χ² is unbounded from above and, by dimensional analysis, can be equivalently characterized by a noise temperature

$T^{\mathrm{w}\mathrm{n}}\equiv m{\left(\frac{{\chi }^2}{\sqrt{\pi }n\sigma }\right)}^{2/3}.(10)$

Therefore, from now on the ST will be characterized by the pair (T^wn, ɛ). In terms of the new parameters, the BFPE, Eq. 3, reads

$\frac{\partial }{\partial t}f\left(\mathbf{v},\omega ;t\right)-\frac{{\nu }^{\mathrm{w}\mathrm{n}}{T}^{\mathrm{w}\mathrm{n}}}{4m}\left[\left(1-\epsilon \right){\left(\frac{\partial }{\partial \mathbf{v}}\right)}^2+\epsilon \frac{2m}{I}{\left(\frac{\partial }{\partial \omega }\right)}^2\right]f\left(\mathbf{v},\omega ;t\right)=J\left[\mathbf{v},\omega |f\right],(11)$

${\nu }^{\mathrm{w}\mathrm{n}}=2n\sigma \sqrt{\frac{\pi T^{\mathrm{w}\mathrm{n}}}{m}}(12)$

being a reference noise-induced collisional frequency.

Inserting the definitions of partial granular temperatures, Eq. 6, into Eq. 11, one obtains

${\partial }_t{T}_{\text{tr}}=-{\xi }_{\text{tr}}{T}_{\text{tr}}+\frac{1-\epsilon }{2}{\nu }^{\mathrm{w}\mathrm{n}}{T}^{\mathrm{w}\mathrm{n}},{\partial }_t{T}_{\text{rot}}=-{\xi }_{\text{rot}}{T}_{\text{rot}}+\epsilon {\nu }^{\mathrm{w}\mathrm{n}}{T}^{\mathrm{w}\mathrm{n}},(13)$

where

${\xi }_{\text{tr}}=-\frac{m}{2n{T}_{\text{tr}}}\int \mathrm{d}\mathbf{v}{\int }_{-\infty }^{\infty }\mathrm{d}\omega v^{2}J\left[\mathbf{v},\omega |f\right],(\mathrm{14a})$

${\xi }_{\text{rot}}=-\frac{I}{n{T}_{\text{rot}}}\int \mathrm{d}\mathbf{v}{\int }_{-\infty }^{\infty }\mathrm{d}\omega {\omega }^{2}J\left[\mathbf{v},\omega |f\right](\mathrm{14b})$

are the translational and rotational energy production rates [78, 80].

In terms of the quantities defined in Eq. 8, Eqs. 13 become

${\partial }_t\tilde{T}=-\zeta \tilde{T}+\frac{1}{3}{\nu }^{\mathrm{w}\mathrm{n}},(\mathrm{15a})$

${\partial }_t\theta =\theta \left({\xi }_{\text{tr}}-{\xi }_{\text{rot}}\right)-{\nu }^{\mathrm{w}\mathrm{n}}\frac{2+\theta }{6\tilde{T}}\left[\theta -\epsilon \left(2+\theta \right)\right],(\mathrm{15b})$

where $\tilde{T}\equiv T/T^{\mathrm{w}\mathrm{n}}$ and

$\zeta \equiv \frac{2}{3}\frac{{\xi }_{\text{tr}}{T}_{\text{tr}}}{T}+\frac{1}{3}\frac{{\xi }_{\text{rot}}{T}_{\text{rot}}}{T}=\frac{2{\xi }_{\text{tr}}+{\xi }_{\text{rot}}\theta }{2+\theta }(16)$

is the cooling rate.

According to Eqs. 15, the steady-state quantities ${\tilde{T}}^{\mathrm{s}\mathrm{t}}=T^{\mathrm{s}\mathrm{t}}/T^{\mathrm{w}\mathrm{n}}$ and ${\theta }^{\mathrm{s}\mathrm{t}}=T_{\text{rot}}^{\mathrm{s}\mathrm{t}}/T_{\text{tr}}^{\mathrm{s}\mathrm{t}}$ satisfy the conditions

${\zeta }^{\mathrm{s}\mathrm{t}}{\tilde{T}}^{\mathrm{s}\mathrm{t}}=\frac{1}{3}{\nu }^{\mathrm{w}\mathrm{n}},\epsilon {\xi }_{\text{tr}}^{\mathrm{s}\mathrm{t}}=\frac{1-\epsilon }{2}{\xi }_{\text{rot}}^{\mathrm{s}\mathrm{t}}{\theta }^{\mathrm{s}\mathrm{t}},(17)$

which imply a balance between collisional cooling and external heating. The steady-state temperature can be used to define a reduced temperature $T^{\ast }\equiv \tilde{T}/{\tilde{T}}^{\mathrm{s}\mathrm{t}}=T/T^{\mathrm{s}\mathrm{t}}$ and a reduced time $t^{\ast }=\frac{1}{2}{\nu }^{\mathrm{s}\mathrm{t}}t$, where

${\nu }^{\mathrm{s}\mathrm{t}}=2n\sigma \sqrt{\frac{\pi T_{\text{tr}}^{\mathrm{s}\mathrm{t}}}{m}}={\nu }^{\mathrm{w}\mathrm{n}}\sqrt{\frac{3{\tilde{T}}^{\mathrm{s}\mathrm{t}}}{2+{\theta }^{\mathrm{s}\mathrm{t}}}}(18)$

is the steady-state collision frequency. More in general, the time-dependent collision frequency is

$\nu \left(t\right)=2n\sigma \sqrt{\frac{\pi T_{\text{tr}}\left(t\right)}{m}}={\nu }^{\mathrm{w}\mathrm{n}}\sqrt{\frac{3\tilde{T}\left(t\right)}{2+\theta \left(t\right)}}={\nu }^{\mathrm{s}\mathrm{t}}G\left(T^{\ast }\left(t\right),\theta \left(t\right)\right),(\mathrm{19a})$

$G\left(T^{\ast },\theta \right)\equiv \sqrt{T^{\ast }\frac{2+{\theta }^{\mathrm{s}\mathrm{t}}}{2+\theta }}.(\mathrm{19b})$

The above collision frequency can be used to nondimensionalize the energy production rates as

${\mu }_{20}\equiv \frac{{\xi }_{\text{tr}}}{\nu },{\mu }_{02}\equiv \frac{{\xi }_{\text{rot}}}{2\nu },(20)$

where

${\mu }_{k\ell }=-\int \mathrm{d}\mathbf{c}{\int }_{-\infty }^{\infty }\mathrm{d}w{c}^k{w}^{\ell }\mathcal{J}\left[\mathbf{c},w|f\right](21)$

are the reduced collisional moments. In Eq. 21, $\mathbf{c}=\mathbf{v}/\sqrt{2T_{\text{tr}}/m}$ and $w=\omega /\sqrt{2T_{\text{root}}/I}$ are the reduced translational and angular velocities, respectively, and $\mathcal{J}\equiv (2{\mathit{T}}_{\mathrm{tr}}/m)\sqrt{2{\mathit{T}}_{\text{root}}/I}/(n\nu )J$ is the reduced collision operator. Thus, the steady-state conditions (17) become

$2\left({\mu }_{20}^{\mathrm{s}\mathrm{t}}+{\mu }_{02}^{\mathrm{s}\mathrm{t}}{\theta }^{\mathrm{s}\mathrm{t}}\right)={\left(\frac{2+{\theta }^{\mathrm{s}\mathrm{t}}}{3{\tilde{T}}^{\mathrm{s}\mathrm{t}}}\right)}^{3/2}\equiv {\gamma }^{\mathrm{s}\mathrm{t}},\epsilon {\mu }_{20}^{\mathrm{s}\mathrm{t}}=\left(1-\epsilon \right){\mu }_{02}^{\mathrm{s}\mathrm{t}}{\theta }^{\mathrm{s}\mathrm{t}}.(22)$

Using these dimensionless definitions, Eqs. 15 yield

$\frac{1}{2}{\partial }_{t^{\ast }}\ln T^{\ast }=-2G\left(T^{\ast },\theta \right)\frac{{\mu }_{20}+{\mu }_{02}\theta }{2+\theta }+\frac{2}{T^{\ast }}\frac{{\mu }_{20}^{\mathrm{s}\mathrm{t}}+{\mu }_{02}^{\mathrm{s}\mathrm{t}}{\theta }^{\mathrm{s}\mathrm{t}}}{2+{\theta }^{\mathrm{s}\mathrm{t}}},(\mathrm{23a})$

$\frac{1}{2}{\partial }_{t^{\ast }}\ln \theta =G\left(T^{\ast },\theta \right)\left({\mu }_{20}-2{\mu }_{02}\right)-\frac{1}{T^{\ast }}\frac{2+\theta }{2+{\theta }^{\mathrm{s}\mathrm{t}}}\frac{{\theta }^{\mathrm{s}\mathrm{t}}}{\theta }\frac{\theta -\epsilon \left(2+\theta \right)}{{\theta }^{\mathrm{s}\mathrm{t}}-\epsilon \left(2+{\theta }^{\mathrm{s}\mathrm{t}}\right)}\left({\mu }_{20}^{\mathrm{s}\mathrm{t}}-2{\mu }_{02}^{\mathrm{s}\mathrm{t}}\right),(\mathrm{23b})$

where use has been made of Eqs. 17, 22.

According to the definition of collisional moments, Eq. 21, they depend on the whole VDF. This implies that Eqs. 23 do not make a closed set of equations. The same applies to the steady-state solution, Eqs. 22. This shortcoming, however, can be circumvented if an approximate closure is applied. This is the subject of section 2.2.

2.2 Maxwellian approximation

In order to get explicit results from Eqs. 22, 23 by using the simplest possible closure, we resort to the two-temperature MA

$f\left(\mathbf{v},\omega \right)\to n\frac{m}{2\pi T_{\text{tr}}}\sqrt{\frac{I}{2\pi T_{\text{rot}}}}\exp \left(-\frac{m{v}^2}{2T_{\text{tr}}}-\frac{I{v}^2}{2T_{\text{rot}}}\right).(24)$

This approximation does a very good job in the three-dimensional case with ɛ = 0 [75] and it is reasonably expected to perform also well in the case of disks with ɛ ≠ 0.

Within this approximation, the relevant collisional moments can be evaluated with the result [78–80].

${\mu }_{20}=\frac{1-{\alpha }^2}{2}+\frac{\mathcal{K}}{2}\left[\frac{\kappa \left(1-\beta \right)}{2}\left(1+\frac{\theta }{\kappa }\right)+1-\theta \right],(\mathrm{25a})$

${\mu }_{02}=\frac{\mathcal{K}}{2}\left[\frac{1-\beta }{2\kappa }\left(1+\frac{\kappa }{\theta }\right)+1-{\theta }^{-1}\right],\mathcal{K}\equiv \kappa \frac{1+\beta }{{\left(1+\kappa \right)}^2}.(\mathrm{25b})$

Solving Eqs. 22, we get the steady-state expressions

${\theta }^{\mathrm{s}\mathrm{t}}=\kappa \left[\frac{2}{\mathcal{K}}\frac{\left(1-{\alpha }^2\right)\epsilon +\mathcal{K}\left(1+\kappa \right)}{\left(1-\beta \right)\left[1-\epsilon \left(1+\kappa \right)\right]+2\kappa }-1\right],{\tilde{T}}^{\mathrm{s}\mathrm{t}}=\frac{2+{\theta }^{\text{st}}}{3{\left({\gamma }^{\mathrm{s}\mathrm{t}}\right)}^{2/3}},(26)$

with

${\gamma }^{\mathrm{s}\mathrm{t}}=1-{\alpha }^2+\mathcal{K}\frac{1+\kappa }{2\kappa }\left(1-\beta \right)\left(\kappa +{\theta }^{\mathrm{s}\mathrm{t}}\right).(27)$

In particular,

$\epsilon =0\Rightarrow {\theta }^{\mathrm{s}\mathrm{t}}=\frac{1+\beta }{2+{\kappa }^{-1}\left(1-\beta \right)},(\mathrm{28a})$

$\epsilon =1\Rightarrow {\theta }^{\mathrm{s}\mathrm{t}}=2\frac{1+{\mathcal{K}}^{-1}\left(1-{\alpha }^2\right)+\frac{\kappa }{2}\left(1-\beta \right)}{1+\beta }.(\mathrm{28b})$

Eq. 28a agrees with a previous result [78]. Notice that, for the special value ɛ = 0, θ^st is independent of the coefficient of normal restitution α both for disks and spheres [47, 75, 78–80]. However, this property is broken down when energy is injected into the rotational degree of freedom (ɛ ≠ 0).

From Eqs. 26 one can observe that θ^st is independent of T^wn and, at given α and β, it is a monotonically increasing function of ɛ. This is physically expected since, by growing ɛ, we are increasing the relative amount of rotational energy injected with respect to the total energy; therefore, it is presumed that the stationary value of T_rot rises with respect to T_tr at fixed T^wn. Thus, the most disparate values of θ^st correspond to ɛ = 1 and 0, their difference being plotted in Figure 2 as a function of α and β.

FIGURE 2

FIGURE 2. Difference θ^st(ɛ = 1) − θ^st(ɛ = 0) [Eq. 28] in the plane (α, β).

It is interesting to note that, in the MA, Eq. 24, one simply has.

$2\left({\mu }_{20}+{\mu }_{02}\theta \right)={\gamma }^{\mathrm{s}\mathrm{t}}+\mathcal{K}\left(1-\beta \right)\frac{1+\kappa }{2\kappa }\left(\theta -{\theta }^{\mathrm{s}\mathrm{t}}\right),(\mathrm{29a})$

${\mu }_{20}-2{\mu }_{02}=\left(1-\epsilon \frac{2+{\theta }^{\mathrm{s}\mathrm{t}}}{{\theta }^{\mathrm{s}\mathrm{t}}}\right)\frac{{\gamma }^{\mathrm{s}\mathrm{t}}}{2}-\mathcal{K}\frac{1+\beta }{4}\left(1+\frac{2}{\theta {\theta }^{\mathrm{s}\mathrm{t}}}\right)\left(\theta -{\theta }^{\mathrm{s}\mathrm{t}}\right).(\mathrm{29b})$

As a consequence, Eqs. 23 can be recast as

$\frac{1}{2}{\partial }_{t^{\ast }}\ln T^{\ast }=\mathrm{\Phi }\left(T^{\ast },\theta \right)\equiv -{\gamma }^{\mathrm{s}\mathrm{t}}\left[\frac{G\left(T^{\ast },\theta \right)}{2+\theta }-\frac{1/T^{\ast }}{2+{\theta }^{\mathrm{s}\mathrm{t}}}\right]-G\left(T^{\ast },\theta \right)\mathcal{K}\left(1-\beta \right)\frac{1+\kappa }{2\kappa }\frac{\theta -{\theta }^{\mathrm{s}\mathrm{t}}}{2+\theta },(\mathrm{30a})$

$\begin{array}{ll}\hfill \frac{1}{2}{\partial }_{t^{\ast }}\ln \theta =& \left(1-\epsilon \frac{2+{\theta }^{\mathrm{s}\mathrm{t}}}{{\theta }^{\mathrm{s}\mathrm{t}}}\right)\frac{{\gamma }^{\mathrm{s}\mathrm{t}}}{2}\left[G\left(T^{\ast },\theta \right)-\frac{1}{T^{\ast }}\frac{2+\theta }{2+{\theta }^{\mathrm{s}\mathrm{t}}}\frac{{\theta }^{\mathrm{s}\mathrm{t}}}{\theta }\frac{\theta -\epsilon \left(2+\theta \right)}{{\theta }^{\mathrm{s}\mathrm{t}}-\epsilon \left(2+{\theta }^{\mathrm{s}\mathrm{t}}\right)}\right]\hfill \\ \hfill & -G\left(T^{\ast },\theta \right)\mathcal{K}\frac{1+\beta }{4}\left(1+\frac{2}{\theta {\theta }^{\mathrm{s}\mathrm{t}}}\right)\left(\theta -{\theta }^{\mathrm{s}\mathrm{t}}\right).\text{(30b)}\hfill \end{array}$

Eqs. 30 make a closed set of two ordinary differential equations that can be (numerically) solved with arbitrary initial conditions $T^{\ast }(0)\equiv T_0^{\ast }$ and θ(0) ≡θ₀. Although the theoretical results have been derived for arbitrary values of the reduced moment of inertia κ, henceforth all the graphs are obtained for the conventional case of a uniform mass distribution of the disks, i.e., $\kappa =\frac{1}{2}$.

Figures 3, 4 show density plots of Φ(T^∗, θ) for α = 0.9 and 0.7, respectively. In each case, two values of β (0 and −0.7) and three values of ɛ (0, 0.5, and 1) are considered. We observe that, typically, Φ(T^∗, θ) increases with increasing θ at fixed T^∗, while it decreases with increasing T^∗ at fixed θ.

FIGURE 3

FIGURE 3. Density plots of Φ(T*, θ) for α = 0.9 and (A) β = 0, ɛ = 0; (B) β = 0, ɛ = 0.5; (C) β = 0, ɛ = 1; (D) β = −0.7, ɛ = 0; (E) β = −0.7, ɛ = 0.5; and (F) β = −0.7, ɛ = 1. The contour lines (solid for Φ < 0, dashed for Φ > 0) are separated by an amount ΔΦ = 0.1. The thick solid line is the locus Φ(T*, θ) = 0. It intercepts the (brown dash-dotted) vertical line T* = 1 at θ = θ^st.

FIGURE 4

FIGURE 4. Same as Figure 3 but for α = 0.7.

2.3 Comparison with simulation results

In order to check the validity of the MA, we have compared our theoretical predictions against DSMC and EDMD simulation results both for transient and steady-state values.

The DSMC algorithm used is based on the one presented in, e.g., [81, 82], and adapted for the model presented in this work. For our DSMC simulations we have dealt with N = 10⁴ particles and chosen a time step Δt = 4 × 10^–5/ν^wn. In addition, the way of implementing the stochastic force and torque in the EDMD code is based on the approximate Green function algorithm [83], as applied to the ST. We have chosen N = 3.6 × 10³ particles, a density nσ² = 5 × 10^–4 (implying a box length of L/σ ≈ 1897.37), and a time step Δt ≈ 4 × 10^–4/ν^wn. No instabilities were observed.

In Figure 5, results for the steady-state values ${\tilde{T}}^{\mathrm{s}\mathrm{t}}$ and θ^st as functions of ɛ and different values of α and β are presented. Simulation results for DSMC and EDMD come from averages over 100 replicas and over 50 data points, once the steady state is ensured to be reached. A very good agreement between DSMC and EDMD with expressions in Eqs. 26 is observed. Whereas θ^st is an increasing function of ɛ, this is not the case, in general, with ${\tilde{T}}^{\mathrm{s}\mathrm{t}}$, as can be observed in Figure 5E.

FIGURE 5

FIGURE 5. Steady-state values of the mean temperature ${\tilde{T}}^{\mathrm{s}\mathrm{t}}$ and the rotational-to-translational temperature ratio θ^st as functions of ɛ. The values of the coefficients of restitution are (A,B): (α, β) = (0.7, 0); (C,D): (α, β) = (0.7, − 0.7); (E,F): (α, β) = (0.9, 0); and (G,H): (α, β) = (0.9, − 0.7). Thick black lines correspond to the theoretical prediction in Eqs. 26, and symbols refer to DSMC (o) and EDMD (×) simulation results for ɛ = 0, 0.1, 0.25, 0.5, 0.8, 0.9, and 1.

As a test of the transient stage, we present in Figure 6 the evolution of T^∗ and θ (starting from a Gaussian-generated VDF with $T_0^{\ast }=1.5$ and θ₀ = 1) for the same choices of α and β as in Figure 5 and for the representative values ɛ = 0 (translational noise only), 0.8 (both translational and rotational noise), and ɛ = 1 (rotational noise only). We observe again an excellent agreement of the MA, Eqs. 30, with simulation results.

FIGURE 6

FIGURE 6. Time evolution of the reduced mean temperature T* and the rotational-to-translational temperature ratio θ for ɛ = 0, 0.8, and 1. The initial conditions are $T_0^{*}=1.5$ and θ₀ = 1 in all the cases. The values of the coefficients of restitution are (A,B): (α, β) = (0.7, 0); (C,D): (α, β) = (0.7, − 0.7); (E,F): (α, β) = (0.9, 0); and (G,H): (α, β) = (0.9, − 0.7). Lines correspond to the theoretical prediction from Eq. 30, and symbols refer to DSMC and EDMD simulation results.

2.4 Temperature overshoot

As illustrated in Figure 6, the evolution of T^∗(t^∗) for certain initial states might experiment an overshoot $T^{\ast }(t_O^{\ast })=1$ at a finite time $t_O^{\ast }$, followed by a minimum, and then relax to the steady state from below. This overshoot effect becomes more pronounced as $\mathrm{\Phi }(T_0^{\ast },{\theta }_0)$ takes more negative values, i.e., as ɛ increases and/or θ₀ decreases.

In general, at a given initial condition $(T_0^{\ast },{\theta }_0)$, there exists a critical value ${\epsilon }_{\mathrm{c}\mathrm{r}}(T_0^{\ast },{\theta }_0)$, such that the T^∗(t^∗) exhibits overshoot if $\epsilon >{\epsilon }_{\mathrm{c}\mathrm{r}}(T_0^{\ast },{\theta }_0)$. The determination of ${\epsilon }_{\mathrm{c}\mathrm{r}}(T_0^{\ast },{\theta }_0)$ within the MA requires the numerical solution of Eqs. 30. We have analyzed those numerical solutions up to t^∗ = 15 since the overshoot typically takes place in the first stage of the evolution. Thus, the numerical value ɛ = ɛ_cr corresponds to $t_O^{\ast }=15$.

Figure 7 shows ${\epsilon }_{\mathrm{c}\mathrm{r}}(T_0^{\ast },{\theta }_0)$ as a function of $T_0^{\ast }$ at a specific value of θ₀, namely the one given by Eq. 28a and for the same four pairs of coefficients of restitution as in Figures 3–6. In each case, the curve ${\epsilon }_{\mathrm{c}\mathrm{r}}(T_0^{\ast },{\theta }_0)$ splits the plane ɛ vs. $T_0^{\ast }$ into two regions: the region above the curve, where the overshoot effect is present, and the one below the curve, where temperature relaxes to the steady-state value from above. Note that the shape of the curve ɛ_cr in Figure 7 associated with the pair (α, β) = (0.7, −0.7) differs in curvature from the curves associated with the other three pairs.

FIGURE 7

FIGURE 7. Phase diagram for the emergence of overshoot. The numerical critical value ${\epsilon }_{\mathrm{c}\mathrm{r}}(T_0^{*},{\theta }_0)$ is presented as a function of $T_0^{*}$ for θ₀ given by Eq. 28a and four different pairs of the coefficients of restitution: (α, β) = (0.7, 0), (0.7, − 0.7), (0.9, 0), and (0.9, − 0.7).

3 Mpemba effect

As already said in Section 1, ME refers to the counterintuitive phenomenon according to which an initially hotter sample of a given fluid relaxes earlier to the steady state than an initially colder one. In a recent paper [44], we distinguished between the thermal ME—where the relaxation process is described by the temperature of the system (second moment of the VDF)— and the entropic ME—where the deviation from the final steady state is monitored by the Kullback–Leibler divergence (thus involving the full VDF). Whereas this distinction is interesting and the relationship between the thermal ME and the entropic ME is not always biunivocal [44], we focus this paper on the thermal version due to its simpler characterization and its relationship with the original results [6]. Morover, only cooling processes will be considered throughout this work.

Let us assume two samples—denoted by A and B—of the same gas, subject to the same noise temperature $T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ and the same splitting parameter ɛ_ref, so that the final steady-state values ${\tilde{T}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}$ and ${\theta }_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}$ will also be the same. Both samples differ in the initial conditions $(T_{0A}^{\ast },{\theta }_{0A})$ and $(T_{0B}^{\ast },{\theta }_{0B})$, respectively, where $T_{0A}^{\ast }>T_{0B}^{\ast }>1$, that is, A refers to the initially hotter sample and we are considering a cooling experiment.

3.1 Standard Mpemba effect

Let us first consider the standard form of thermal ME [42, 44–47], where both $T_A^{\ast }(t^{\ast })$ and $T_B^{\ast }(t^{\ast })$ cross over at a certain crossing time $t_c^{\ast }$ and then relax from above, i.e., $T_A^{\ast }(t_c^{\ast })=T_B^{\ast }(t_c^{\ast })$ and $T_B^{\ast }(t^{\ast })>T_A^{\ast }(t^{\ast })>1$ for $t^{\ast }>t_c^{\ast }$. The emergence of ME can be subdued to the appearance of a crossing time (or an odd number of them) in the absence of any overshoot effect in the thermal evolution [44]. We will refer to this situation as the standard ME (SME). According to the discussion in Section 2.4, the SME implies that ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}<\min \{{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0A}^{\ast },{\theta }_{0A}),{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})\}$. The case when a temperature overshoot takes place will be discussed in Section 3.2.

It can be reasonably expected that a necessary condition for the occurrence of the SME is that the initial slope is smaller in sample A than in sample B, i.e.,

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

Note here that the usual situation is that both slopes are negative, in which case $|\mathrm{\Phi }(T_{0A}^{\ast },{\theta }_{0A})|>|\mathrm{\Phi }(T_{0B}^{\ast },{\theta }_{0B})|$. Of course, Eq. 31 is not sufficient for the SME since the latter also depends on how close $T_{0A}^{\ast }$ and $T_{0B}^{\ast }$ are and how far both initial temperatures are from unity. From Figures 3, 4 we can conclude that, in general, the inequality (31) is best satisfied if θ_0A ≪θ_0B.

3.2 Overshoot Mpemba effect

The emergence of the temperature overshoot described in Section 2.4 makes the crossover criterion employed in the SME become meaningless. Imagine that such a crossover takes place with $T_A^{\ast }(t_c^{\ast })=T_B^{\ast }(t_c^{\ast })>1$, but then $T_B^{\ast }(t^{\ast })$ relaxes from above while $T_A^{\ast }(t^{\ast })$ overshoots the steady-state value, $T_A^{\ast }(t_O^{\ast })=1$. It is then possible that $T_A^{\ast }(t^{\ast })$ relaxes (from below) later than $T_B^{\ast }(t^{\ast })$. In that case, the initially hotter sample (A) would reach the steady state later than the initially colder sample (B), thus contradicting the existence of a ME, despite the crossover.

Reciprocally, imagine that $T_A^{\ast }(t^{\ast })$ and $T_B^{\ast }(t^{\ast })$ never cross each other but $T_B^{\ast }(t^{\ast })$ overshoots the steady-state value and then relaxes from below. It is now possible that the initially hotter sample (A) reaches the steady state earlier than the initially colder sample (B), thus qualifying as a ME, despite the absence of any crossover. We will refer to this phenomenon as overshoot ME (OME) [44]. From the discussion in Section 2.4 we conclude that the OME requires ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}>{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})$.

To characterize the existence of OME without a thermal crossover, we adopt the quantity [44].

$\mathfrak{D}\left(T^{\ast }\left(t\right)\right)\equiv T^{\ast }\left(t\right)-1-\ln T^{\ast }\left(t\right).(32)$

This quantity is (except for a factor) the Kullback–Leibler divergence of the Maxwellian VDF given by Eq. 24 (with T_rot/T_tr → θ^st) with respect to the steady-state Maxwellian. Note that $\mathfrak{D}(T^{\ast })$ is a positive-definite convex function of T^∗. Therefore, we can define the OME by the crossover of $\mathfrak{D}(T_A^{\ast })$ and $\mathfrak{D}(T_B^{\ast })$ with, however, $T_A^{\ast }<1$ and $T_A^{\ast }<T_B^{\ast }$.

In order to look for OME, the necessary condition for SME [given by Eq. 31] must be reversed. That is,

$\mathrm{\Phi }\left(T_{0A}^{\ast },{\theta }_{0A}\right)>\mathrm{\Phi }\left(T_{0B}^{\ast },{\theta }_{0B}\right).(33)$

Establishing the most favorable conditions for OME is not as simple as just declaring the opposite of the SME condition. Firstly, we want for the colder sample to overshoot as much as possible the steady state, so that the relaxation from below is retarded maximally. This reasoning is translated into the condition of highly negative initial slope $\mathrm{\Phi }(T_{0B}^{\ast },{\theta }_{0B})$, which implies small θ_0B (Figures 3, 4). On the other hand, two competing phenomena exist for the hotter sample: we want to either avoid any overshoot or force it to be as weak as possible, but we also want the relaxation to be faster than in the colder sample. Therefore, one needs θ_0A > θ_0B but, for very large values of θ_0A, one might not find OME due to a slower relaxation of the hotter sample.

3.3 Initial preparation protocols

In order to study the absence or existence of ME, one needs to specify the initial conditions $(T_{0A}^{\ast },{\theta }_{0A})$ and $(T_{0B}^{\ast },{\theta }_{0B})$ of both samples. In previous studies [42, 44, 47] the values of θ_0A and θ_0B were freely chosen, without a specific reference to a previous protocol to initially prepare the samples.

Most of the interest of the present work resides in the proposal of protocols to generate the initial states of the samples involved in a ME experiment. The protocols are based on the proposed ST, and the initial states will be generated by assuming prior thermostat values $(T_i^{\mathrm{w}\mathrm{n}},{\epsilon }_i)$, i = A, B, and allowing both samples to reach their respective steady states $(T_i^{\mathrm{s}\mathrm{t}},{\theta }_i^{\mathrm{s}\mathrm{t}})$ before switching to the common posterior thermostat $(T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}},{\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}})$ at t = 0. The values of the prior thermostats will be chosen to optimize the necessary conditions (31) and (33) for SME and OME, respectively.

According to Eqs. 26, 27, the ratios between the prior and posterior noise temperatures for desired values of $T_{0i}^{\ast }=T_i^{\mathrm{s}\mathrm{t}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}$, i = A, B, where $T_i^{\mathrm{s}\mathrm{t}}=T_i^{\mathrm{w}\mathrm{n}}{\tilde{T}}^{\mathrm{s}\mathrm{t}}({\epsilon }_i)$ and $T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}=T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}{\tilde{T}}^{\mathrm{s}\mathrm{t}}({\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}})$, are

$\frac{T_i^{\mathrm{w}\mathrm{n}}}{T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}}=T_{0i}^{\ast }\frac{2+{\theta }_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}}{2+{\theta }_i^{\mathrm{s}\mathrm{t}}}{\left(\frac{{\gamma }_i^{\mathrm{s}\mathrm{t}}}{{\gamma }_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}}\right)}^{2/3},i=A,B.(34)$

3.3.1 Protocol for the standard Mpemba effect

In this case, we want to have θ_0A << θ_0B. According to Eqs. 26, and as observed in Figure 5, θ^st is an increasing function of ɛ and independent of T^wn. Therefore, the most disparate values of θ_0A and θ_0B are obtained if the prior thermostats of samples A and B have ɛ_A = 0 and ɛ_B = 1, respectively. According to Figure 2, SME would be stronger and/or easier to find for lower values of β at fixed α and for lower values of α at fixed β. In addition, in order to define a cooling process, we need to choose proper values of $T_A^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ and $T_B^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ [Eq. 34], such that $T_{0A}^{\ast }=T_A^{\mathrm{s}\mathrm{t}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}>T_{0B}^{\ast }=T_B^{\mathrm{s}\mathrm{t}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}>1$. Finally, one must fix ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}<\min \{{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0A}^{\ast },{\theta }_{0A}),{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})\}$ to prevent any possible overshoot. This minimum of the critical rotational-to-total noise intensity parameter is expected to be ${\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0A}^{\ast },{\theta }_{0A})$ because it corresponds to a more negative initial slope.

Thus, the designed protocol for SME reads as follows (Figure 8 for an illustrative scheme):

1) Start by fixing ɛ_A = 0 and ɛ_B = 1, in order to ensure θ_0A < θ_0B.

2) Choose $T_A^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ and $T_B^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$, such that $T_{0A}>T_{0B}>T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}$.

3) Let both samples evolve and reach the steady states corresponding to their respective prior thermostats. These steady states will play the role of the initial conditions for our ME experiment.

4) Switch the values of the thermostats of both samples to a common reference pair of values $(T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}},{\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}})$, such that no overshoot is expected, that is, ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}<\min \{{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0A}^{\ast },{\theta }_{0A}),{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})\}$. This thermostat switch fixes the origin of time, t = 0.

5) Finally, let both samples evolve and reach a common steady state.

FIGURE 8

FIGURE 8. Scheme of the protocol for SME. (A) Choice of the splitting parameters ɛ. (B) Choice of the noise temperatures T^wn.

For given (α, β), the numerical solutions of Eqs. 30 for different values of $T_{0A}^{\ast }$ and $T_{0B}^{\ast }$—and with θ_0A and θ_0B given by Eqs. 28a, 28b, respectively—can be analyzed to determine whether SME is present or not. This provides the phase diagram presented in Figure 10A for ɛ_ref = 0.1 and some pairs of coefficients of restitution.

3.3.2 Protocol for the overshoot Mpemba effect

In the OME case, it is convenient to have θ_0A ≫θ_0B, so that the adopted choices of ɛ for the prior thermostats are the reverse of those of SME, i.e., ɛ_A = 1 and ɛ_B = 0. Whereas in Section 3.2 we commented that the best situation for the OME is not always the opposite to that of the SME, the above choice helps us avoid or weaken a possible overshoot for the hotter sample. Again, in order to define a cooling process, we need to choose proper values of $T_A^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ and $T_B^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ [Eq. 34], and such that $T_{0A}>T_{0B}>T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}$. Finally, ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}>{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})$ to ensure overshoot of $T_B^{\ast }(t^{\ast })$.

In analogy with the SME case, a protocol for observing OME is designed as follows (Figure 9):

1) Start by fixing ɛ_A = 1 and ɛ_B = 0, in order to ensure θ_0A > θ_0B.

2) Choose $T_A^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ and $T_B^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$, such that $T_{0A}>T_{0B}>T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{s}\mathrm{t}}$.

3) Let both samples evolve and reach the steady states corresponding to their respective prior thermostats. These steady states will play the role of the initial conditions for our ME experiment.

4) Switch the values of the thermostats of both samples to a common reference pair of values $(T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}},{\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}})$, such that overshoot is ensured, that is, ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}>{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})$. This thermostat switch fixes the origin of time, t = 0.

5) Finally, let both samples evolve and reach a common steady state.

FIGURE 9

FIGURE 9. Same as in Figure 8, but for OME.

Figure 10B shows a phase diagram for the occurrence of OME (with ɛ_ref = 0.9) for the same pairs of coefficients of restitution as before.

FIGURE 10

FIGURE 10. Phase diagrams in the plane $T_{0A}^{*}/T_{0B}^{*}$ vs. $T_{0B}^{*}$ for the emergence of (A) SME (with ɛ_ref = 0.1) and (B) OME (with ɛ_ref = 0.9). Four different pairs of the coefficients of restitution are considered: (α, β) = (0.7, − 0.7), (0.7, 0), (0.9, − 0.7), and (0.9, 0). The insets show magnified views of the indicated regions.

3.4 Comparison with simulation results

In order to check the initialization protocols for detecting SME and OME, we have run DSMC and EDMD simulations. The simulation details are the same as introduced in Section 2.3. Simulation points correspond to the average over ensembles of 100 replicas. Again, no instabilities were observed.

3.4.1 Standard Mpemba effect

Figure 11 presents results for the SME protocol introduced in Section 3.3.1. As we can observe, the theoretical predictions agree very well with DSMC and EDMD simulation data.

FIGURE 11

FIGURE 11. Time evolution of the reduced mean temperature T* and the rotational-to-translational temperature ratio θ for SME initialization protocol, as described in Section 3.3.1. (A,B): (α, β) = (0.7, 0) and ɛ_ref = 0.1; (C,D): (α, β) = (0.7, − 0.7) and ɛ_ref = 0.1; (E,F); (α, β) = (0.9, 0) and ɛ_ref = 0.6; and (G,H): (α, β) = (0.9, − 0.7) and ɛ_ref = 0.1. Thick and dashed lines correspond to the theoretical prediction from Eq. 30, dotted lines represented the steady state value, and symbols refer to DSMC and EDMD simulation results. The inset in panel (E) shows the evolution of the temperature difference $T_A^{*}-T_B^{*}$.

The chosen initial temperature conditions are $(T_{0A}^{\ast },T_{0B}^{\ast })=(\mathrm{3,2})$ for Figure 11A, $(T_{0A}^{\ast },T_{0B}^{\ast })=(\mathrm{4,2})$ for Figures 11C,G, and $(T_{0A}^{\ast },T_{0B}^{\ast })=(\mathrm{3,2.92})$ for Figure 11E. Moreover, ɛ_ref = 0.1 for all cases, except for the case of Figures 11E,F, that is chosen to be ɛ_ref = 0.6. The values of $T_A^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ and $T_B^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ are given in each case by Eq. 34 to ensure the desired initial temperatures. All these cases avoid overshoot, in agreement with the case θ₀ = θ^st(ɛ = 0) shown in Figure 7. Moreover, the cases with ɛ_ref = 0.1 are inside the region where SME is predicted to be present in Figure 10A.

The value ɛ_ref = 0.6 for the system (α = 0.9, β = 0) was chosen instead of ɛ_ref = 0.1 to avoid the need of taking very high initial temperatures (Figure 10A) and also to prevent overshoot (Figure 7). The price paid for this choice of coefficients of restitution is that the difference θ_0A− θ_0B is not too high, see Figures 2, 11F. Then, initial temperature values are restricted to be very similar and the crossing characterizing the SME is much less pronounced than in the other cases. However, the inset in Figure 11E shows a well defined change of sign of the difference $T_A^{\ast }-T_B^{\ast }$.

3.4.2 Overshoot Mpemba effect

The theoretical results stemming from the OME protocol introduced in Section 3.3.2 are compared with simulations in Figure 12, again with an excellent agreement.

FIGURE 12

FIGURE 12. Time evolution of the reduced mean temperature T*, the entropy-like quantity $\mathfrak{D}$, and the rotational-to-translational temperature ratio θ for OME initialization protocol described in Section. 3.3.2 with ɛ_ref = 0.9. (A–C): (α, β) = (0.7, 0); (D–F): (α, β) = (0.7, − 0.7); (G–I); (α, β) = (0.9, 0); and (J–L): (α, β) = (0.9, − 0.7). Thick and dashed lines correspond to the theoretical prediction from Eqs. 30, dotted lines represented the steady state value, and symbols refer to DSMC and EDMD simulation results.

The initial temperatures are $(T_{0A}^{\ast },T_{0B}^{\ast })=(1.22,1.1)$ for Figure 12A, $(T_{0A}^{\ast },T_{0B}^{\ast })=(\mathrm{2,1.5})$ for Figures 12D,J, and $(T_{0A}^{\ast },T_{0B}^{\ast })=(1.2,1.1)$ for Figure 12G. Moreover, a common value ɛ_ref = 0.9 is chosen to ensure that overshoot appears for the colder samples, as depicted in Figure 7, and that OME is expected from the phase diagram shown in Figure 10B.

Whereas no crossing between $T_A^{\ast }$ and $T_B^{\ast }$ is observed in Figure 12, such a crossing is present between the entropy-like quantities ${\mathfrak{D}}_A$ and ${\mathfrak{D}}_B$ [Eq. 32]. As explained in Section 3.2, this characterizes the OME, where the initially colder system relaxes later to the steady state.

4 Concluding remarks

In this work, we have studied the homogeneous states of a dilute granular gas made of inelastic and rough hard disks lying on a two-dimensional plane. The inelasticity and roughness are mathematically described by constant coefficients of normal (α) and tangential (β) restitution, respectively. In order to avoid frozen long-time limiting states, the disks are assumed to be heated by the stochastic force- and torque-based ST. This novel stochastic thermostat injects energy to both translational and rotational degrees of freedom. Each specific thermostat of this type is univocally determined by its associated noise temperature, T^wn ≥ 0 [Eq. 10], and the rotational-to-total noise intensity, 0 ≤ ɛ ≤ 1 [Eq. 9].

The system is assumed to be fully described by the instantaneous one-particle VDF, its dynamics being then given by the BFPE, Eq. (11). It is known that the steady-state VDF is not a Maxwellian, as occurs for ɛ = 0 [75], and this is even more the case with the instantaneous transient VDF. However, these nonGaussianities are expected to be small enough as to approximate the VDF by a Maxwellian, as previously done for ɛ = 0 [47]. Therefore, we have worked under the two-temperature MA introduced in Section 2.2, which allows us to account for the dynamics of the systems just in terms of the reduced total granular temperature, T^∗, and the rotational-to-translational temperature ratio, θ, according to Eqs. 30. The steady-state values can be explicitly expressed in terms of the mechanical properties of the disks, Eqs. 26. Steady and transient states predicted by the MA are tested via DSMC ans EDMD, with very good agreement, as observed in Figures 5, 6. This reinforces the validity of our approach. As expected, θ^st is an increasing function of ɛ and independent of T^wn, whereas ${\tilde{T}}^{\mathrm{s}\mathrm{t}}\equiv T^{\mathrm{s}\mathrm{t}}/T^{\mathrm{w}\mathrm{n}}$ can be a nonmonotonic function of ɛ (Figure 5).

The main core of this work has been the description of ME in cooling processes in this system, with special emphasis on the elaboration of preparation protocols for the generation of the initial states. We noted that, if T^∗ − 1 does not change its sign during the evolution, the usual form of ME, SME, can emerge if θ_0A ≪θ_0B, A being the initially hotter sample. This SME is characterized by a single crossing (or, in general, an odd number of them) between the temperature curves, thus inducing that the initially colder system, B, relaxes more slowly toward the final steady state. However, we have realized that an overshoot (or change of sign of T^∗ − 1) might appear if the rotational-to-total noise intensity is larger than a certain critical value ${\epsilon }_{\mathrm{c}\mathrm{r}}(T_0^{\ast },{\theta }_0)$, as shown for some cases in Figure 7. Therefore, the SME is no longer ensured and, following a recent work [44], ME might appear without a crossing between $T_A^{\ast }$ and $T_B^{\ast }$ but with a crossing between the entropy-like quantities ${\mathfrak{D}}_A$ and ${\mathfrak{D}}_B$ [see Eq. 32]. This peculiar type of ME is termed OME and is described in Section 3.2. According to the MA, the most favorable condition for the OME changes from the SME one to θ_0A ≫θ_0B.

Protocols for generating initial conditions to observe both SME and OME have been presented in Sections 3.3.1, 3.3.2, respectively. While reminiscent of the protocols previously considered in the case of sheared inertial suspensions [53], the protocols proposed here represent novel instructions to elaborate a ME experiment in homogeneous states of granular gaseous systems. We have based those protocols on the steady states of the ST, taking advantage of the increase of θ^st as ɛ increases. Thus, to guarantee the biggest possible difference |θ_0B− θ_0A| we have fixed prior thermalization processes with (ɛ_A, ɛ_B) = (0, 1) for SME and (ɛ_A, ɛ_B) = (1, 0) for OME. Moreover, in this prior thermalization stage, $T_A^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$ and $T_B^{\mathrm{w}\mathrm{n}}/T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}}$, where $(T_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{w}\mathrm{n}},{\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}})$ characterizes the posterior thermostat, are chosen to fine-tune both $T_{0A}^{\ast }$ and $T_{0B}^{\ast }$. It is crucial to choose ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}<\min \{{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0A}^{\ast },{\theta }_{0A}),{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})\}$ for SME to avoid overshoot, and ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}>{\epsilon }_{\mathrm{c}\mathrm{r}}(T_{0B}^{\ast },{\theta }_{0B})$ for OME to ensure the overshoot of the initially colder sample. These protocols, together with the MA predictions, allowed us to elaborate phase diagrams for SME and OME emergences, as depicted in Figure 10.

The theoretical descriptions of the SME and OME initialization protocols have been tested by DMSC and EDMD simulations, finding a very good agreement between theory and simulations, as observed in Figures 11, 12, respectively.

One can then conclude that, despite its simplicity, the MA captures very well the dynamics for this system. In turn, this implies that the SME and OME protocols for the initial-state preparation described in this paper are trustworthy.

We expect that this work will be useful to the ME community in the search for practical protocols able to generate adequate initial states. Whereas this work is restricted to a specific physical model, the phenomenology that arises might be general in Mpemba-like effects in out-of-equilibrium statistical mechanical systems. In fact, the overshoot mechanism introduced has been previously observed in other gaseous models [44], and it is somehow reminiscent of the supercooling hypothesis for the ME in water [19]. Furthermore, the two-parameter-based protocols proposed in this paper might inspire the development of protocols for Mpemba-like effects in other systems.

To sum up, the designs of the protocols can be based on a very obvious statement: generate the initial conditions from steady states for specific values of the parameters of the thermostat; then, change the values of those thermostat parameters and let the system relax to the new steady state. In general, if the temperature (or the energy) of the physical system in a certain model is coupled to n independent control variables (i.e., the dynamics is described by n + 1 differential equations), one should look for a proper choice of n + 1 parameters of the model that will determine directly the value of each variable at the steady state.

Finally, given the simplicity of the hard-disk system studied in this paper, we hope it can be experimentally realizable, thus opening up the possibility of reproducing the ME by an adequate control of the external forcing mechanisms.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/amegiasf/MpembaSplitting.

Author contributions

AM and AS contributed to the conception and design of the study. AM performed the computer simulations and wrote the first draft of the manuscript. Both authors contributed to manuscript revision, read, and approved the submitted version.

Funding

The authors acknowledge financial support from Grant No. PID2020-112936GB-I00 funded by MCIN/AEI/10.13039/501100011033, and from Grants No. IB20079 and No. GR21014 funded by Junta de Extremadura (Spain) and by ERDF “A way of making Europe.” AM is grateful to the Spanish Ministerio de Ciencia, Innovación y Universidades for a predoctoral fellowship FPU2018-3503.

Acknowledgments

The authors are grateful to the computing facilities of the Instituto de Computación Científica Avanzada of the University of Extremadura (ICCAEx), where the simulations were run.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1.Ross WD, editor. The works of Aristotle (translated into English under the editorship of W.D. Ross), III. London: Oxford Clarendon Press (1931).

Google Scholar

2. Bacon F. Novum organum scientiarum (in Latin). In: Translated into English in: The new organon. Cambridge: Cambridge Univ Press (2000). Under the editorship of L. Jardine, M. Silverthorne.

Google Scholar

3. Descartes R. Discours de la méthode pour bien conduire sa raison, et chercher la vérité dans les sciences (in French). In: Translated into English in: Discourse on method, optics, geometry, and meteorology. Hackett Publishing Company (2001). Under the editorship of P. J. Olscamp.

Google Scholar

4. Newton I. VII. Scala graduum caloris. Phil Trans R Soc (1701) 22:824–9. doi:10.1098/rstl.1700.0082

CrossRef Full Text | Google Scholar

5. Newton I. Isaaci Newtoni Opera quae exstant omnia, 4 (1782). p. 403–7. Londini: excudebat Joannes Nichols.

Google Scholar

6. Mpemba EB, Osborne DG. Cool? Phys Educ (1969) 4:312–175. doi:10.1088/0031-9120/4/3/312

CrossRef Full Text | Google Scholar

7. Osborne DG. Mind on ice. Phys Educ (1979) 14:414–7. doi:10.1088/0031-9120/14/7/313

CrossRef Full Text | Google Scholar

8. Elkin S (Editor). The 100 greatest unsolved mysteries. New York: Cavendish Square (2018):14.

9. Kell GS. The freezing of hot and cold water. Am J Phys (1969) 37:564–5. doi:10.1119/1.1975687

CrossRef Full Text | Google Scholar

10. Firth I. Cooler? Phys Educ (1971) 6:310–41. doi:10.1088/0031-9120/6/1/310

CrossRef Full Text | Google Scholar

11. Deeson E. Cooler-lower down. Phys Educ (1971) 6:311–44. doi:10.1088/0031-9120/6/1/311

CrossRef Full Text | Google Scholar

12. Frank FC. The Descartes–Mpemba phenomenon. Phys Educ (1974) 9:121–284. doi:10.1088/0031-9120/9/4/121

CrossRef Full Text | Google Scholar

13. Gallear R. The Bacon–Descartes–Mpemba phenomenon. Phys Educ (1974) 9:114–490. doi:10.1088/0031-9120/9/7/114

CrossRef Full Text | Google Scholar

14. Walker J. Hot water freezes faster than cold water. Why does it do so? Sci Am (1977) 237:246–57. doi:10.1038/scientificamerican0977-246

CrossRef Full Text | Google Scholar

15. Freeman M. Cooler still—an answer? Phys Educ (1979) 14:417–21. doi:10.1088/0031-9120/14/7/314

CrossRef Full Text | Google Scholar

16. Kumar K. Mpemba effect and 18th century ice-cream. Phys Educ (1980) 15:101–268. doi:10.1088/0031-9120/15/5/101

CrossRef Full Text | Google Scholar

17. Hanneken JW. Mpemba effect and cooling by radiation to the sky. Phys Educ (1981) 16:102–7. doi:10.1088/0031-9120/16/1/102

CrossRef Full Text | Google Scholar

18. Wojciechowski B, Owczarek I, Bednarz G. Freezing of aqueous solutions containing gases. Cryst Res Technol (1988) 23:843–8. doi:10.1002/crat.2170230702

CrossRef Full Text | Google Scholar

19. Auerbach D. Supercooling and the Mpemba effect: When hot water freezes quicker than cold. Am J Phys (1995) 63:882–5. doi:10.1119/1.18059

CrossRef Full Text | Google Scholar

20. Knight CA. The Mpemba effect: The freezing times of hot and cold water. Am J Phys (1996) 64:524. doi:10.1119/1.18275

CrossRef Full Text | Google Scholar

21. Maciejewski PK. Evidence of a convective instability allowing warm water to freeze in less time than cold water. J Heat Transfer (1996) 118:65–72. doi:10.1115/1.2824069

CrossRef Full Text | Google Scholar

22. Jeng M. The Mpemba effect: When can hot water freeze faster than cold? Am J Phys (2006) 74:514–22. doi:10.1119/1.2186331

CrossRef Full Text | Google Scholar

23. Esposito S, De Risi R, Somma L. Mpemba effect and phase transitions in the adiabatic cooling of water before freezing. Physica A: Stat Mech its Appl (2008) 387:757–63. doi:10.1016/j.physa.2007.10.029

CrossRef Full Text | Google Scholar

24. Katz JI. When hot water freezes before cold. Am J Phys (2009) 77:27–9. doi:10.1119/1.2996187

CrossRef Full Text | Google Scholar

25. Vynnycky M, Mitchell SL. Evaporative cooling and the Mpemba effect. Heat Mass Transfer (2010) 46:881–90. doi:10.1007/s00231-010-0637-z

CrossRef Full Text | Google Scholar

26. Brownridge JD. When does hot water freeze faster then cold water? A search for the Mpemba effect. Am J Phys (2011) 79:78–84. doi:10.1119/1.3490015

CrossRef Full Text | Google Scholar

27. Vynnycky M, Maeno N. Axisymmetric natural convection-driven evaporation of hot water and the Mpemba effect. Int J Heat Mass Transf (2012) 55:7297–311. doi:10.1016/j.ijheatmasstransfer.2012.07.060

CrossRef Full Text | Google Scholar

28. Balážovič M, Tomášik B. The Mpemba effect, Shechtman’s quasicrystals and student exploration activities. Phys Educ (2012) 47:568–73. doi:10.1088/0031-9120/47/5/568

CrossRef Full Text | Google Scholar

29. Zhang X, Huang Y, Ma Z, Zhou Y, Zhou J, Zheng W, et al. Hydrogen-bond memory and water-skin supersolidity resolving the Mpemba paradox. Phys Chem Chem Phys (2014) 16:22995–3002. doi:10.1039/C4CP03669G

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Vynnycky M, Kimura S. Can natural convection alone explain the Mpemba effect? Int J Heat Mass Transf (2015) 80:243–55. doi:10.1016/j.ijheatmasstransfer.2014.09.015

CrossRef Full Text | Google Scholar

31. Sun CQ. Behind the Mpemba paradox. Temperature (2015) 2:38–9. doi:10.4161/23328940.2014.974441

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Balážovič M, Tomášik B. Paradox of temperature decreasing without unique explanation. Temperature (2015) 2:61–2. doi:10.4161/23328940.2014.975576

CrossRef Full Text | Google Scholar

33. Romanovsky AA. Which is the correct answer to the Mpemba puzzle? Temperature (2015) 2:63–4. doi:10.1080/23328940.2015.1009800

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Jin J, Goddard WA. Mechanisms underlying the Mpemba effect in water from molecular dynamics simulations. J Phys Chem C (2015) 119:2622–9. doi:10.1021/jp511752n

CrossRef Full Text | Google Scholar

35. Ibekwe RT, Cullerne JP. Investigating the Mpemba effect: When hot water freezes faster than cold water. Phys Educ (2016) 51:025011. doi:10.1088/0031-9120/51/2/025011

CrossRef Full Text | Google Scholar

36. Gijón A, Lasanta A, Hernández ER. Paths towards equilibrium in molecular systems: The case of water. Phys Rev E (2019) 100:032103. doi:10.1103/PhysRevE.100.032103

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bechhoefer J, Kumar A, Chétrite R. A fresh understanding of the Mpemba effect. Nat Rev Phys (2021) 3:534–5. doi:10.1038/s42254-021-00349-8

CrossRef Full Text | Google Scholar

38. Burridge HC, Linden PF. Questioning the Mpemba effect: hot water does not cool more quickly than cold. Sci Rep (2016) 6:37665. doi:10.1038/srep37665

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Burridge HC, Hallstadius O. Observing the Mpemba effect with minimal bias and the value of the Mpemba effect to scientific outreach and engagement. Proc R Soc A (2020) 476:20190829. doi:10.1098/rspa.2019.0829

CrossRef Full Text | Google Scholar

40. Elton DC, Spencer PD. Pathological water science — four examples and what they have in common. In: Gadomski A, editor. Biomechanical and related systems. Biologically-inspired systems, 17. Cham: Springer (2021). p. 155–69.

CrossRef Full Text | Google Scholar

41. Żuk PJ, Makuch K, Hłyst R, Maciołek A. Transient dynamics in the outflow of energy from a system in a nonequilibrium stationary state. Phys Rev E (2022) 105:054133. doi:10.1103/PhysRevE.105.054133

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Santos A, Prados A. Mpemba effect in molecular gases under nonlinear drag. Phys Fluids (2020) 32:072010. doi:10.1063/5.0016243

CrossRef Full Text | Google Scholar

43. Patrón A, Sánchez-Rey B, Prados A. Strong nonexponential relaxation and memory effects in a fluid with nonlinear drag. Phys Rev E (2021) 104:064127. doi:10.1103/PhysRevE.104.064127

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Megías A, Santos A, Prados A. Thermal versus entropic Mpemba effect in molecular gases with nonlinear drag. Phys Rev E (2022) 105:054140. doi:10.1103/PhysRevE.105.054140

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Gómez González R, Khalil N, Garzó V. Mpemba-like effect in driven binary mixtures. Phys Fluids (2021) 33:053301. doi:10.1063/5.0050530

CrossRef Full Text | Google Scholar

46. Lasanta A, Vega Reyes F, Prados A, Santos A. When the hotter cools more quickly: Mpemba effect in granular fluids. Phys Rev Lett (2017) 119:148001. doi:10.1103/physrevlett.119.148001

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Torrente A, López-Castaño MA, Lasanta A, Vega Reyes F, Prados A, Santos A. Large Mpemba-like effect in a gas of inelastic rough hard spheres. Phys Rev E (2019) 99:060901. doi:10.1103/PhysRevE.99.060901

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Biswas A, Prasad VV, Raz O, Rajesh R. Mpemba effect in driven granular Maxwell gases. Phys Rev E (2020) 102:012906. doi:10.1103/PhysRevE.102.012906

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Mompó E, López Castaño MA, Torrente A, Vega Reyes F, Lasanta A. Memory effects in a gas of viscoelastic particles. Phys Fluids (2021) 33:062005. doi:10.1063/5.0050804

CrossRef Full Text | Google Scholar

50. Gómez González R, Garzó V. Time-dependent homogeneous states of binary granular suspensions. Phys Fluids (2021) 33:093315. doi:10.1063/5.0062425

CrossRef Full Text | Google Scholar

51. Biswas A, Prasad VV, Rajesh R. Mpemba effect in an anisotropically driven granular gas. EPL (2021) 136:46001. doi:10.1209/0295-5075/ac2d54

CrossRef Full Text | Google Scholar

52. Biswas A, Prasad VV, Rajesh R. Mpemba effect in anisotropically driven inelastic Maxwell gases. J Stat Phys (2022) 186:45. doi:10.1007/s10955-022-02891-w

CrossRef Full Text | Google Scholar

53. Takada S, Hayakawa H, Santos A. Mpemba effect in inertial suspensions. Phys Rev E (2021) 103:032901. doi:10.1103/PhysRevE.103.032901

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Takada S. Homogeneous cooling and heating states of dilute soft-core gases undernonlinear drag. EPJ Web Conf (2021) 249:04001. doi:10.1051/epjconf/202124904001

CrossRef Full Text | Google Scholar

55. Baity-Jesi M, Calore E, Cruz A, Fernandez LA, Gil-Narvión JM, Gordillo-Guerrero A, et al. The Mpemba effect in spin glasses is a persistent memory effect. Proc Natl Acad Sci U S A (2019) 116:15350–5. doi:10.1073/pnas.1819803116

PubMed Abstract | CrossRef Full Text | Google Scholar

56. González-Adalid Pemartín I, Mompó E, Lasanta A, Martín-Mayor V, Salas J. Slow growth of magnetic domains helps fast evolution routes for out-of-equilibrium dynamics. Phys Rev E (2021) 104:044114. doi:10.1103/PhysRevE.104.044114

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Teza G, Yaacoby R, Raz O. Relaxation shortcuts through boundary coupling. arXiv [Preprint] (2021). doi:10.48550/arXiv.2112.10187

CrossRef Full Text | Google Scholar

58. Vadakkayila N, Das SK. Should a hotter paramagnet transform quicker to a ferromagnet? Monte Carlo simulation results for ising model. Phys Chem Chem Phys (2021) 23:11186–90. doi:10.1039/d1cp00879j

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Yang ZY, Hou JX. Non-markovian Mpemba effect in mean-field systems. Phys Rev E (2020) 101:052106. doi:10.1103/PhysRevE.101.052106

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Yang ZY, Hou JX. Mpemba effect of a mean-field system: The phase transition time. Phys Rev E (2022) 105:014119. doi:10.1103/PhysRevE.105.014119

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Greaney PA, Lani G, Cicero G, Grossman JC. Mpemba-like behavior in carbon nanotube resonators. Metall Mat Trans A (2011) 42:3907–12. doi:10.1007/s11661-011-0843-4

CrossRef Full Text | Google Scholar

62. Ahn YH, Kang H, Koh DY, Lee H. Experimental verifications of Mpemba-like behaviors of clathrate hydrates. Korean J Chem Eng (2016) 33:1903–7. doi:10.1007/s11814-016-0029-2

CrossRef Full Text | Google Scholar

63. Schwarzendahl FJ, Löwen H. Anomalous cooling and overcooling of active systems. Phys. Rev. Lett. (Forthcoming 2022).

Google Scholar

64. Carollo F, Lasanta A, Lesanovsky I. Exponentially accelerated approach to stationarity in Markovian open quantum systems through the Mpemba effect. Phys Rev Lett (2021) 127:060401. doi:10.1103/PhysRevLett.127.060401

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Lu Z, Raz O. Nonequilibrium thermodynamics of the Markovian Mpemba effect and its inverse. Proc Natl Acad Sci U S A (2017) 114:5083–8. doi:10.1073/pnas.1701264114

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Klich I, Raz O, Hirschberg O, Vucelja M. Mpemba index and anomalous relaxation. Phys Rev X (2019) 9:021060. doi:10.1103/PhysRevX.9.021060

CrossRef Full Text | Google Scholar

67. Chétrite R, Kumar A, Bechhoefer J. The metastable Mpemba effect corresponds to a non-monotonic temperature dependence of extractable work. Front Phys (2021) 9:654271. doi:10.3389/fphy.2021.654271

CrossRef Full Text | Google Scholar

68. Busiello DM, Gupta D, Maritan A. Inducing and optimizing Markovian Mpemba effect with stochastic reset. New J Phys (2021) 23:103012. doi:10.1088/1367-2630/ac2922

CrossRef Full Text | Google Scholar

69. Lin J, Li K, He J, Ren J, Wang J. Power statistics of Otto heat engines with the Mpemba effect. Phys Rev E (2022) 105:014104. doi:10.1103/PhysRevE.105.014104

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Holtzman R, Raz O. Landau theory for the Mpemba effect through phase transitions. arXiv [Preprint] (2022). doi:10.48550/arXiv.2204.03995

CrossRef Full Text | Google Scholar

71. Kumar A, Bechhoefer J. Exponentially faster cooling in a colloidal system. Nature (Lond.) (2020) 584:64–8. doi:10.1038/s41586-020-2560-x

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Kumar A, Chétrite R, Bechhoefer J. Anomalous heating in a colloidal system. Proc Natl Acad Sci U S A (2022) 119:e2118484119. doi:10.1073/pnas.2118484119

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Kovacs AJ. Transition vitreuse dans les polymères amorphes. Etude phénoménologique. Fortschr Hochpolym.-Forsch (1963) 3:394–507. doi:10.1007/BFb0050366

CrossRef Full Text | Google Scholar

74. Kovacs AJ, Aklonis JJ, Hutchinson JM, Ramos AR. Isobaric volume and enthalpy recovery of glasses. II. A transparent multiparameter theory. J Polym Sci Polym Phys Ed (1979) 17:1097–162. doi:10.1002/pol.1979.180170701

CrossRef Full Text | Google Scholar

75. Vega Reyes F, Santos A. Steady state in a gas of inelastic rough spheres heated by a uniform stochastic force. Phys Fluids (2015) 27:113301. doi:10.1063/1.4934727

CrossRef Full Text | Google Scholar

76. Megías A, Santos A. Translational and angular velocity cumulants in granular gases of inelastic and rough hard disks or spheres. In preparation (2022).

Google Scholar

77. Garzó V. A kinetic theory approach to granular gaseous flows. Switzerland: Springer Nature (2019).Granular gaseous flows.

Google Scholar

78. Megías A, Santos A. Driven and undriven states of multicomponent granular gases of inelastic and rough hard disks or spheres. Granul Matter (2019) 21:49. doi:10.1007/s10035-019-0901-y

CrossRef Full Text | Google Scholar

79. Megías A, Santos A. Energy production rates of multicomponent granular gases of rough particles. A unified view of hard-disk and hard-sphere systems. AIP Conf Proc (2019) 2132:080003. doi:10.1063/1.5119584

CrossRef Full Text | Google Scholar

80. Santos A. Interplay between polydispersity, inelasticity, and roughness in the freely cooling regime of hard-disk granular gases. Phys Rev E (2018) 98:012904. doi:10.1103/PhysRevE.98.012904

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Bird GA. The DSMC method. Scotts Valley, CA: CreateSpace Independent Publishing Platform (2013).

Google Scholar

82. Montanero JM, Santos A. Computer simulation of uniformly heated granular fluids. Granul Matter (2000) 2:53–64. doi:10.1007/s100350050035

CrossRef Full Text | Google Scholar

83. Scala A. Event-driven Langevin simulations of hard spheres. Phys Rev E (2012) 86:026709. doi:10.1103/PhysRevE.86.026709

PubMed Abstract | CrossRef Full Text | Google Scholar