# Mpemba effect in driven granular gases: role of distance measures

Apurba Biswas<EMAIL_ADDRESS>The Institute of Mathematical Sciences,
C.I.T. Campus, Taramani, Chennai 600113, India  Homi Bhabha National
Institute, Training School Complex, Anushakti Nagar, Mumbai 400094, India  V .
V. Prasad<EMAIL_ADDRESS>Department of Physics, Cochin University of
Science and Technology, Kochi 682022, India  R. Rajesh<EMAIL_ADDRESS>The Institute of Mathematical Sciences, C.I.T. Campus, Taramani, Chennai
600113, India  Homi Bhabha National Institute, Training School Complex,
Anushakti Nagar, Mumbai 400094, India

###### Abstract

Mpemba effect refers to the counterintuitive effect where a system which is
initially further from the final steady state equilibrates faster than an
identical system that is initially closer. The closeness to the final state is
defined in terms of a distance measure. For driven granular systems, the
Mpemba effect has been illustrated in terms of an ad-hoc measure of mean
kinetic energy as the distance function. In this paper, by studying four
different distance measures based on the mean kinetic energies as well as
velocity distribution, we show that the Mpemba effect depends on the
definition of the measures.

##  I Introduction

Usual experience tells us that a hotter object takes a longer time to cool
than a similar object initially at a lower temperature. Such phenomena are
well described using Newtons’s cooling law where rate of cooling is directly
proportional to the difference in temperature between the object and its
surrounding. Mpemba effect refers to a counterintuitive phenomenon where a
system which is initially at a higher temperature equilibrates faster than a
similar system which is initially at a lower intermediate temperature, when
both the systems are quenched to the same low temperature. This effect was
first observed in the case of water where freezing time was found to be lower
for water that is initially at higher temperature [1, 2, 3, 4, 5, 6]. However,
the effect is not limited to water and has been experimentally shown to exist
in many other physical systems such as magnetic alloys [7], clathrate hydrates
[8], polylactides [9], carbon nanotube resonators [10] and colloidal systems
[11, 12, 13], thus suggesting that the Mpemba effect is a much more general
anomalous relaxation phenomenon that can be studied in numerous physical
systems. Theoretical studies on the Mpemba effect have focused on spin systems
[14, 15, 16, 17, 18, 19], systems undergoing phase transitions [17, 20, 21],
Markovian systems with only a few states [22, 23], systems of single particle
diffusing in a potential [24, 25, 26, 27], active systems [28], spin glasses
[29], molecular gases in contact with a thermal reservoir [30, 31, 32, 33],
quantum systems [34] and granular systems [35, 36, 37, 38, 39, 40, 41].

The protocol followed to test for the existence of the Mpemba effect considers
two identical systems which are prepared in different initial states. Both the
systems are then quenched to a common final state, and the Mpemba effect is
said to be present if the system that is initially further away from the final
state equilibrates faster. To quantify which initial state is further away as
well as crossing of trajectories, a distance measure for points in the phase
space has to be defined.

For systems relaxing to thermal equilibrium, the distance to the final
equilibrium state is measured in terms of the probabilities of the different
states [22]. It was argued that any distance measure satisfying the following
properties should result in a unique definition of the Mpemba effect: (1) As
the system relaxes toward thermal equilibrium, the distance function should
decrease with time, (2) for three temperatures $T_{h}>T_{c}>T_{b}$, the
distance from $T_{b}$ is larger for $T_{h}$ compared to $T_{c}$, i.e., the
distance measure should be a monotonically increasing function of temperature
and (3) the distance function should be a continuous, convex function of
probability distribution [22]. Let $\pi$ denote the equilibrium Boltzmann
distribution, and $p_{i}(t)$ the probability of state $i$ at time $t$. Then,
examples of such distances are “entropic distance”
$D_{e}=\sum_{i}(p_{i}(t)-\pi_{i})E_{i}/T_{b}+p_{i}(t)\ln
p_{i}(t)-\pi_{i}\ln\pi_{i}$ [22], total variation distance
$L_{1}(t)=\sum_{i}|p_{i}(t)-\pi_{i}|$ [11] and Kullback-Leibler (KL)
divergence defined as $D_{KL}(t)=\sum_{i}p_{i}(t)\ln(p_{i}(t)/\pi_{i})$ [42].

The calculation of the above definitions of distance rely on knowing the
probability distribution $p(t)$ at all times. Analytical calculation of $p(t)$
is possible only for exactly solvable problems which, in the context of the
Mpemba effect, are restricted to simple single particle systems [11, 12, 13]
and to systems with only a few states [22]. For interacting many-particle
systems and more so in the context of out of equilibrium systems, the distance
measures defined in the probability space are, however, inaccessible through
direct measurements in experiments or computationally expensive to measure in
simulations. The more natural choices for the distance measures that have been
used in experiments are moments of the distribution which are directly
observable like mean energy per spin [29, 17], magnetization [7, 17], granular
temperature [35, 36, 37, 38, 39, 40, 41], etc. In addition, the generalization
to far from equilibrium systems, like driven granular systems is not clear.
There are no conclusive studies about the correspondence or similarity between
the directly observable measures and the distance measures defined in
probability space. Thus, it is not evident whether the Mpemba effect that has
been established using these proxy measures for distances is unique or its
existence depends on the distance measure used.

To address the above issue, we study the dependence of the Mpemba effect on
the choice of different measures in driven granular systems, a prototypical
interacting, many-particle far from equilibrium system. The advantage of
granular system is that it allows exact analysis for two-point correlation
functions, at least in the linearized regime. Being an athermal system,
granular temperature or the mean kinetic energy of the system has been used to
track the evolution of the system with higher granular temperature being
considered further away from the final steady state [35, 37, 36, 38, 39, 40].
However, it is not known a priori whether mean kinetic energy correctly
predicts the distance of the initial states from the final steady states. In
the following, we summarize the several results that demonstrates the
existence of the Mpemba effect in granular systems. For a system of smooth
monodispersed particles [35, 37, 38], the Mpemba effect is due to the coupling
of the translational granular temperature with the excess kurtosis of the
velocity distribution function or correlations between the velocities of the
different particles. In the case of rough granular gas [36], the Mpemba effect
is due to the coupling of granular temperatures defined for translational and
rotational degrees of freedom. In these cases, for the Mpemba effect to exist,
however the initial states need to be different from steady states. Existence
of the Mpemba effect for evolution of systems from initial steady states could
be realised in the context of binary inelastic gases or when inelastic gases
are driven anisotropically. For the case of binary gas, Mpemba effect was
traced to the energy-exchange between the energies of the smaller and bigger
particles [38], while in anisotropically driven gas [39, 40], the Mpemba
effect is due to the coupling between the granular temperatures along $x$ and
$y$ directions. However, in all the above mentioned analysis, the time
evolution of the system is projected onto only one of the variables which is
the total granular temperature defined by the second moment of velocity
distribution of the system.

In this paper, in addition to total energy, we introduce other measures such
as Manhattan measure ($\mathcal{L}_{1}$), Euclidean measure
($\mathcal{L}_{2}$) and KL divergence which can describe the evolution of the
system in the phase space of all the relevant variables. We perform the
analysis in the setup of anisotropically driven granular gas as well as driven
inelastic Maxwell gas [39, 40]. We derive the criteria for the existence of
the Mpemba effect with the various measures and determine the region of phase
space which shows the Mpemba effect. We show that these phase diagrams are
non-universal in the sense that they depend on the measure used.

The remainder of the paper is organized as follows. We first define the model
of anisotropically driven granular and Maxwell gases in Sec. II. In Sec. III,
we derive the time evolution of the two-point correlation functions, exactly
for the Maxwell gas, and in the linearized regime for the granular gas. It
allows us to exactly calculate their steady state values which characterize
the steady state. In Sec. IV, we define the Mpemba effect and the various
measures that we use to characterize the distance between steady states. In
Sec. V, we present the results for the Mpemba effect using the various
measures. We compare the various measures for the Mpemba effect in terms of
their phase diagrams. Section VI contains the conclusion and its implications.

##  II Models

In this paper, we analyze two models for the driven granular gas in two
dimensions: inelastic Maxwell model and the hard disc granular gas model.
Consider a collection of identical, inelastic particles. The velocities of
these particles change in time due to inelastic binary collisions that are
momentum conserving. The new velocities $\boldsymbol{v}^{\prime}_{i}$ and
$\boldsymbol{v}^{\prime}_{j}$ of the two colliding particles $i$ (velocity
$\boldsymbol{v}_{i}$) and $j$ (velocity $\boldsymbol{v}_{j}$) after the
collision is given by

$\displaystyle\begin{split}\boldsymbol{v}^{\prime}_{i}=\boldsymbol{v}_{i}-\alpha[(\boldsymbol{v}_{i}-\boldsymbol{v}_{j}).\boldsymbol{\hat{e}}]\boldsymbol{\hat{e}},\\\
\boldsymbol{v}^{\prime}_{j}=\boldsymbol{v}_{j}+\alpha[(\boldsymbol{v}_{i}-\boldsymbol{v}_{j}).\boldsymbol{\hat{e}}]\boldsymbol{\hat{e}},\end{split}$
(1)

where

$\displaystyle\alpha=\frac{1+r}{2},$ (2)

and $r$ is the coefficient of restitution and $\boldsymbol{\hat{e}}$ is the
unit vector along the line joining the centers of the particles at contact. In
addition, the particles are driven anisotropically at a constant rate,
$\lambda_{d}$. We describe the driving for each model separately.

Inelastic Maxwell gas: In the inelastic Maxwell gas model, the collision rates
of the particles are independent of the relative velocities of the colliding
particles. Thus, the spatial degrees of freedom are irrelevant. These
simplifications lead to the model being analytically tractable, yet keeping
essential characteristics of more complicated realistic models.

In addition to the collision rules described by Eq. (1), each particle is
driven with rate $\lambda_{d}$. When particle $i$ is driven, its velocity
$\boldsymbol{v}_{i}$ is changed to $\boldsymbol{v}^{\prime}_{i}$, given by

$\displaystyle v^{\prime}_{ix,y}=-r_{wx,y}v_{ix,y}+\eta_{ix,y},$ (3)

where $r_{wx,y}$ are scalar parameters associated with the driving along the
$x$ and $y$ directions, and $\boldsymbol{\eta}$ is a noise taken from a fixed
distribution $\phi({\boldsymbol{\eta}})$ with a finite second moment denoted
by

$\displaystyle\sigma^{2}_{k}=\int^{\infty}_{-\infty}d\eta_{k}\eta^{2}_{k}\phi(\boldsymbol{\eta}),\quad
k=x,y.$ (4)

When $\sigma^{2}_{x}\neq\sigma^{2}_{y}$, the driving is anisotropic and leads
to particles having an anisotropic velocity distribution. The physical
motivations for the form of driving in Eq. (3) may be found in Refs. [43, 44].

Let $P(\boldsymbol{v},t)$ denote the probability that a randomly chosen
particle has velocity $\boldsymbol{v}$ at time $t$. Its time evolution is
given by

$\displaystyle\frac{dP(\boldsymbol{v},t)}{dt}$
$\displaystyle=\lambda_{c}\int\int\int
d\boldsymbol{\hat{e}}d\boldsymbol{v}_{1}d\boldsymbol{v}_{2}P(\boldsymbol{v}_{1},t)P(\boldsymbol{v}_{2},t)\delta(\boldsymbol{v}_{1}-\alpha[(\boldsymbol{v}_{1}-\boldsymbol{v}_{2}).\boldsymbol{\hat{e}}]\boldsymbol{\hat{e}}-\boldsymbol{v})-\lambda_{c}P(\boldsymbol{v},t)$
$\displaystyle+\lambda_{d}\int\int
d\boldsymbol{\eta}d\boldsymbol{v}_{1}\Phi(\boldsymbol{\eta})P(\boldsymbol{v}_{1},t)\delta[-r_{wx}v_{1x}+\eta_{x}-v_{x}]\delta[-r_{wy}v_{1y}+\eta_{y}-v_{y}]-\lambda_{d}P(\boldsymbol{v},t),$
(5)

where the first and second terms on the right hand side describe the gain and
loss terms due to collisions while the third and fourth terms describe the
gain and loss terms due to driving.

Hard disc granular gas: In the hard disc granular gas, the collision rates are
proportional to the relative velocities of the colliding particles, as
expected for ballistic motion. The time evolution of the velocity distribution
function $f(\boldsymbol{v},t)$, defined as the number density of particles
having velocity $\boldsymbol{v}$ at time $t$ is described using the Enskog-
Boltzmann equation [45]

$\displaystyle\frac{\partial}{\partial t}f(\boldsymbol{{v}},t)=\chi
I(f,f)+\Big{(}\frac{\sigma^{2}_{x}}{2}\frac{\partial^{2}}{\partial
v^{2}_{x}}+\frac{\sigma^{2}_{y}}{2}\frac{\partial^{2}}{\partial
v^{2}_{y}}\Big{)}f(\boldsymbol{{v}},t),$ (6)

where $\chi$ is the pair correlation function [46], $I(f,f)$ (see Appendix B,
Eq. (60) for details) is the collision integral which accounts for the gain
and loss terms in Eq. (5) due to collisions, and $\sigma^{2}_{x}$ and
$\sigma^{2}_{y}$ are the variances or strengths of the white noise along the
$x$\- and $y$-directions respectively.

The driving terms in Eq. (6), corresponding to the diffusion-like terms, are
related to the driving for the Maxwell gas given in Eq. (3). By Taylor
expanding the terms related to driving in Eq. (5), and truncating for small
noise, it is easily shown that the driving terms in Eq. (6) corresponds to the
special case $r_{wx}=r_{wy}=1$ [43, 44, 39].

Note that we have considered a spatially homogeneous system to describe the
system such that spatial degrees of freedom are ignored. Moreover, we have
also assumed molecular chaos hypothesis to use product measure for the joint
velocity distribution function in the collision integral.

##  III Characterizing the steady states

In this section, we define the relevant two point correlation functions for
both the models. The time evolutions of these two point correlations were
already derived in Refs. [39, 40]. However, we summarize the derivations and
compute their steady state values as they form an integral part of the present
analysis. Sections III.1 and III.2 contain the derivations for the inelastic
Maxwell model and the hard disc granular gas model respectively. The text
follows closely the text in Refs. [39, 40].

###  III.1 Inelastic Maxwell model

For the inelastic Maxwell model, we are interested in the following two-point
correlation functions:

$\displaystyle E_{x}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum^{N}_{i=1}\langle v^{2}_{ix}(t)\rangle,$
$\displaystyle E_{y}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum^{N}_{i=1}\langle v^{2}_{iy}(t)\rangle,$
$\displaystyle E_{xy}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum^{N}_{i=1}\langle v_{ix}(t)v_{iy}(t)\rangle,$
$\displaystyle C_{x}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\sum^{N}_{i=1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\langle v_{ix}(t)v_{jx}(t)\rangle,$ (7)
$\displaystyle C_{y}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\sum^{N}_{i=1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\langle v_{iy}(t)v_{jy}(t)\rangle,$ $\displaystyle
C_{xy}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\sum^{N}_{i=1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\langle v_{ix}(t)v_{jy}(t)\rangle,$

where $E_{x}(t)$ and $E_{y}(t)$ denote the mean kinetic energies of the
particles along $x$ and $y$ directions respectively. $E_{xy}(t)$ denote the
correlations between $v_{x}$ and $v_{y}$ of the same particle whereas
$C_{x}(t)$, $C_{y}(t)$ and $C_{xy}(t)$ denote the velocity-velocity
correlations between pairs of particles. The time evolution of these
correlation functions is derived using Eq. (5) and is compactly written in the
form

$\displaystyle\frac{d\boldsymbol{\tilde{\Sigma}}(t)}{dt}=\boldsymbol{\tilde{R}}\boldsymbol{\tilde{\Sigma}}(t)+\boldsymbol{\tilde{D}},$
(8)

where the column vectors $\boldsymbol{\tilde{\Sigma}}(t)$ and
$\boldsymbol{\tilde{D}}$ are given by:

$\displaystyle\begin{split}&\boldsymbol{\tilde{\Sigma}}(t)=[E_{x}(t),E_{y}(t),E_{xy}(t),C_{x}(t),C_{y}(t),C_{xy}(t)]^{T},\\\
&\boldsymbol{\tilde{D}}=[\lambda_{d}\sigma^{2}_{x},\lambda_{d}\sigma^{2}_{y},0,0,0,0]^{T}.\end{split}$
(9)

The components of the matrix $\boldsymbol{\tilde{R}}$ are given in Appendix A.
In the steady state, in the thermodynamic limit, the velocity-velocity
correlations and the correlation between $v_{x}$ and $v_{y}$ of the same
particle, i.e., $E_{xy}$ vanish as shown in Ref. [40]. Moreover, if these
correlations are zero for the initial state, then it remains zero for all
times. In that case, we only write for the time evolution of the non-zero mean
kinetic energies, i.e., $E_{x}$ and $E_{y}$ in a compact form as

$\frac{d\boldsymbol{\Sigma}(t)}{dt}=~{}\boldsymbol{R}\boldsymbol{\Sigma}(t)+\boldsymbol{S},$
(10)

where

$\displaystyle\boldsymbol{\Sigma}(t)=\begin{bmatrix}E_{x}(t),E_{y}(t)\end{bmatrix}^{T},$
(11)
$\displaystyle\boldsymbol{S}=\begin{bmatrix}\lambda_{d}\sigma^{2}_{x},\lambda_{d}\sigma^{2}_{y}\end{bmatrix}^{T},$
(12)

and $\boldsymbol{R}$ is a $2\times 2$ matrix, whose entries are given by

$\displaystyle\begin{split}R_{11}=&\frac{3}{4}\lambda_{c}\alpha^{2}-\lambda_{c}\alpha-\lambda_{d}(1-r^{2}_{wx}),\quad
R_{12}=\frac{\lambda_{c}}{4}\alpha^{2},\\\
R_{22}=&\frac{3}{4}\lambda_{c}\alpha^{2}-\lambda_{c}\alpha-\lambda_{d}(1-r^{2}_{wy}),\quad
R_{21}=\frac{\lambda_{c}}{4}\alpha^{2}.\\\ \end{split}$ (13)

We will be using a different set of variables than $E_{x}$ and $E_{y}$ for the
analysis. We define the new set of variables, namely the total kinetic energy,
$E_{tot}$, and the difference of energies, $E_{dif}$, as:

$\displaystyle E_{tot}=E_{x}+E_{y},$ (14) $\displaystyle E_{dif}=E_{x}-E_{y}.$
(15)

The time evolution equations for $E_{tot}$ and $E_{dif}$ can be expressed,
starting from Eq. (10), as

$\frac{d\boldsymbol{E}(t)}{dt}=-~{}\boldsymbol{\chi}\boldsymbol{E}(t)+\boldsymbol{D},$
(16)

where

$\displaystyle\boldsymbol{E}(t)=\begin{bmatrix}E_{tot}(t),E_{dif}(t)\end{bmatrix}^{T},$
(17)
$\displaystyle\boldsymbol{D}=\begin{bmatrix}\lambda_{d}(\sigma^{2}_{x}+\sigma^{2}_{y}),\lambda_{d}(\sigma^{2}_{x}-\sigma^{2}_{y})\end{bmatrix}^{T},$
(18)

and $\boldsymbol{\chi}$ is a $2\times 2$ matrix with the components of the
matrix given by

$\displaystyle\begin{split}\chi_{11}&=\frac{2\lambda_{c}\alpha(1-\alpha)+\lambda_{d}(2-r^{2}_{wx}-r^{2}_{wy})}{2},\\\
\chi_{22}&=\frac{\lambda_{c}\alpha(2-\alpha)+\lambda_{d}(2-r^{2}_{wx}-r^{2}_{wy})}{2},\\\
\chi_{12}&=\frac{\lambda_{d}(r^{2}_{wy}-r^{2}_{wx})}{2},\quad\chi_{21}=\frac{\lambda_{d}(r^{2}_{wy}-r^{2}_{wx})}{2}.\end{split}$
(19)

Equation (16) can be solved exactly by linear decomposition using the
eigenvalues $\lambda_{\pm}$ of $\boldsymbol{\chi}$:

$\displaystyle\lambda_{\pm}=$
$\displaystyle\frac{1}{4}\Big{[}2\lambda_{d}(2-r^{2}_{wx}-r^{2}_{wy})+\alpha\lambda_{c}(4-3\alpha)$
(20)
$\displaystyle\pm\sqrt{4\lambda^{2}_{d}(r^{2}_{wy}-r^{2}_{wx})^{2}+\alpha^{4}\lambda^{2}_{c}}\Big{]}.$

It is straightforward to show that $\lambda_{\pm}>0$ with
$\lambda_{+}>\lambda_{-}$. The solution for $E_{tot}(t)$ and $E_{dif}(t)$ is

$\displaystyle\begin{split}E_{tot}(t)-E^{st}_{tot}&=K_{+}e^{-\lambda_{+}t}+K_{-}e^{-\lambda_{-}t},\\\
E_{dif}(t)-E^{st}_{dif}&=L_{+}e^{-\lambda_{+}t}+L_{-}e^{-\lambda_{-}t},\end{split}$
(21)

where $E^{st}_{tot}$ and $E^{st}_{dif}$ are steady state values of
$E_{tot}(t)$ and $E_{dif}(t)$ respectively. The coefficients
$K_{+},K_{-},L_{+}$ and $L_{-}$ along with $E^{st}_{tot}$ and $E^{st}_{dif}$
are given in Appendix A.

###  III.2 Hard disc granular gas model

For the model of driven hard disc granular gas whose evolution is described
using Eq. (6), the mean kinetic energies along the $x$ and $y$ directions are
defined as

$\displaystyle E_{i}(t)=\frac{2}{n}\int
d\boldsymbol{v}\frac{1}{2}mv^{2}_{i}f(\boldsymbol{{v}},t),\quad i=x,y,$ (22)

where $n=\int d\boldsymbol{v}f(\boldsymbol{{v}},t)$ is the number density, $m$
is the mass of the particles. In order to derive the time evolution of the
mean kinetic energies using Eq. (6), we assume an anisotropic form of the
Gaussian for the velocity distribution function (which is valid when the
system is weakly inelastic) as described in Ref. [39],

$f(\boldsymbol{v},t)=\frac{mn}{2\pi\sqrt{E_{x}(t)E_{y}(t)}}\exp\left[-\frac{mv^{2}_{x}}{2E_{x}(t)}-\frac{mv^{2}_{y}}{2E_{y}(t)}\right].$
(23)

With this approximation for $f(\boldsymbol{v},t)$, the time evolution of
$E_{tot}$ and $E_{dif}$ [defined as in Eqs. (14) and (15)] form a coupled set
of non-linear differential equations and is given by

$\displaystyle\begin{split}\frac{\partial}{\partial
t}E_{tot}(t)=\mathcal{F}(E_{tot},E_{dif}),\\\ \frac{\partial}{\partial
t}E_{dif}(t)=\mathcal{G}(E_{tot},E_{dif}).\end{split}$ (24)

The functional forms for $\mathcal{F}(E_{tot},E_{dif})$ and
$\mathcal{G}(E_{tot},E_{dif})$ are given in Appendix B. In order to perform an
exact analysis, we linearize the non-linear equations by considering only the
initial states that are close to the final steady state. We define $\delta
E_{tot}(t)=E_{tot}(t)-E^{st}_{tot}$ and $\delta
E_{dif}(t)=E_{dif}(t)-E^{st}_{dif}$ as the time-dependent deviation of the
energies from the stationary state values. The linearized form for the non-
linear differential equations [see Eq. (24)] about the stationary state values
obey

$\frac{d}{dt}\begin{bmatrix}\delta E_{tot}(t)\\\ \delta
E_{dif}(t)\end{bmatrix}=-~{}\boldsymbol{\chi}\begin{bmatrix}\delta
E_{tot}(t)\\\ \delta E_{dif}(t)\end{bmatrix},$ (25)

where $\boldsymbol{\chi}$ is a $2\times 2$ matrix whose coefficients are given
in the Appendix B. The solutions for $\delta E_{tot}(t)$ and $\delta
E_{dif}(t)$ are then given by

$\displaystyle\begin{split}\delta
E_{tot}(t)=K_{+}e^{-\lambda_{+}t}+K_{-}e^{-\lambda_{-}t},\\\ \delta
E_{dif}(t)=L_{+}e^{-\lambda_{+}t}+L_{-}e^{-\lambda_{-}t},\end{split}$ (26)

where the coefficients $K_{+},K_{-},L_{+}$ and $L_{-}$ are given by

$\displaystyle\begin{split}K_{+}&=\frac{1}{\gamma}\Big{[}\chi_{12}\delta
E_{dif}(0)-(\lambda_{-}-\chi_{11})\delta E_{tot}(0)\Big{]},\\\
K_{-}&=\frac{1}{\gamma}\Big{[}-\chi_{12}\delta
E_{dif}(0)+(\lambda_{+}-\chi_{11})\delta E_{tot}(0)\Big{]},\\\
L_{+}&=\frac{1}{\gamma}\Big{[}(\lambda_{+}-\chi_{11})\delta E_{dif}(0)\\\
&-\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}\delta
E_{tot}(0)\Big{]},\\\
L_{-}&=\frac{1}{\gamma}\Big{[}-(\lambda_{-}-\chi_{11})\delta E_{dif}(0)\\\
&+\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}\delta
E_{tot}(0)\Big{]}.\end{split}$ (27)

Here, $\lambda_{\pm}$ are the eigenvalues of $\boldsymbol{\chi}$ and
$\gamma=\lambda_{+}-\lambda_{-}$.

##  IV  Mpemba effect and distance measures in phase space

We first define the protocol that we will follow for illustrating the Mpemba
effect. Consider two systems $P$ and $Q$ which have identical parameters
except for the driving strengths. Both the systems, in their respective steady
states, are then quenched to a common steady state. This is achieved by
instantaneously changing the driving strengths of $P$ and $Q$ to the common
driving strength of the final steady state, keeping all the other parameters
of both the systems fixed. The two initial steady states $P$ and $Q$ differ in
their initial distance (to be appropriately defined) from the final steady
state. The evolution of the systems $P$ and $Q$ correspond to two different
trajectories in phase space. Then the Mpemba effect is said to exist if the
trajectory that was initially at a larger distance from the final steady state
approaches the final state faster than the trajectory that was initially at a
shorter distance. For granular systems, several distance measures have been
used to explore the Mpemba effect. We derive the criterion for the existence
of the Mpemba effect for the different measures and for each measure, we
illustrate the Mpemba effect. We also ask how much the Mpemba effect depends
on the distance measure being used.

The steady state of the two models for driven granular systems considered in
this paper is completely specified by the velocity distribution
$P({\boldsymbol{v}})$. Distance between two probability distributions may be
defined in terms of an information theoretic quantity known as Kullback-
Leibler (KL) divergence [47, 48, 49]. However, the velocity distribution
cannot be solved exactly and therefore the KL divergence becomes unwieldy for
studying Mpemba effect, though we will study this numerically. Instead, the
steady state of the system has been parametrised by the moments of the
velocity. Note that the equations for the two point correlation functions
close among themselves for both the Maxwell gas as well as the linearized
granular gas. This motivates using the second moments of velocity distribution
to characterize the steady states. For the system of anisotropically driven
granular gas, the mean kinetic energies, $E_{x}$ and $E_{y}$ along $x$\- and
$y$\- directions respectively, are different. Thus, the total energy,
$E_{tot}=E_{x}+E_{y}$ and the difference of energies, $E_{dif}=E_{x}-E_{y}$
serve as the appropriate variable to describe a state of the system. Thus, a
steady state is characterized by $(E_{tot},E_{dif})$.

In this section, we introduce different measures that have been used to define
the distance between two steady states of the system. When the steady state is
defined through $(E_{tot},E_{dif})$, the distance measures can be defined in
terms of the difference of the total energy of the two states (Sec. IV.1), in
terms of two dimensional Euclidean distance (Sec. IV.2) or two-dimensional
Manhattan distance (Sec. IV.3) of the phase space variables of the two states.
In terms of the velocity distribution, the convergence to the steady state can
be characterized through KL divergence (Sec. IV.4). The different measures
track the temporal evolution of the system as it evolves from an initial state
to a final state, and we derive the conditions for the Mpemba effect to exist.

###  IV.1 Total energy

The most common variable that is used in literature of driven granular gases
to describe its state is the mean kinetic energy. The existence of the Mpemba
effect in driven granular systems has been shown in previous studies [38, 39,
40] in terms of these variables. In this section, we briefly discuss the
condition for the existence of the Mpemba effect in an anisotropically driven
granular gas. We consider two systems $P$ and $Q$ whose initial steady states
are denoted by $[E^{P}_{tot}(0),E^{P}_{dif}(0)]$ and
$[E^{Q}_{tot}(0),E^{Q}_{dif}(0)]$ respectively. Here, the distance of the
initial states compared to the final state is measured in terms of the total
energy. The initial steady states for the two systems are prepared such that
$E^{P}_{tot}(0)>E^{Q}_{tot}(0)$. Both the systems are then quenched to a
common steady state having the total energy lower than the initial total
energies of both the systems.

The Mpemba effect is present if the two trajectories $E^{P}_{tot}(t)$ and
$E^{Q}_{tot}(t)$ cross each other at some finite time $t=\tau$ at which

$E^{P}_{tot}(\tau)=E^{Q}_{tot}(\tau).$ (28)

To obtain the value of $\tau$, we equate the total energies of $P$ and $Q$
using either Eq. (21) or Eq. (26) depending on the inelastic Maxwell gas or
hard disc granular gas respectively to obtain

$K^{P}_{+}e^{-\lambda_{+}\tau}+K^{P}_{-}e^{-\lambda_{-}\tau}=K^{Q}_{+}e^{-\lambda_{+}\tau}+K^{Q}_{-}e^{-\lambda_{-}\tau},$
(29)

whose solution is

$\tau=\frac{1}{\lambda_{+}-\lambda_{-}}\ln\Big{[}\frac{K^{P}_{+}-K^{Q}_{+}}{K^{Q}_{-}-K^{P}_{-}}\Big{]}.$
(30)

In terms of the parameters of the initial steady states, $\tau$ reduces to

$\tau=\frac{1}{\lambda_{+}-\lambda_{-}}\ln\Big{[}\frac{\chi_{12}\Delta
E_{dif}-(\lambda_{-}-\chi_{11})\Delta E_{tot}}{\chi_{12}\Delta
E_{dif}-(\lambda_{+}-\chi_{11})\Delta E_{tot}}\Big{]},$ (31)

where

$\displaystyle\begin{split}&\Delta E_{tot}=E^{P}_{tot}(0)-E^{Q}_{tot}(0),\\\
&\Delta E_{dif}=E^{P}_{dif}(0)-E^{Q}_{dif}(0).\end{split}$ (32)

For the Mpemba effect to be present, we require that $\tau>0$. Since
$\lambda_{+}>\lambda_{-}$, the argument of logarithm in Eq. (31) should be
greater than one.

Figure 1 illustrates the existence of the Mpemba effect where the trajectories
of the initial states leading to final steady state are defined in terms of
the total energy and their crossing time is given by Eq. (31).

Figure 1: The time evolution of anisotropically driven (a) inelastic Maxwell
and (b) hard disc granular gas is illustrated in terms of the mean kinetic
energy, $E_{tot}$ for two systems $P$ and $Q$. The initial conditions for the
inelastic Maxwell gas in (a) are $E^{P}_{tot}(0)=1.148$,
$E^{Q}_{tot}(0)=0.92$, $E^{P}_{dif}(0)=-0.595$ and $E^{Q}_{dif}(0)=0.455$. The
choice of the other parameters defining the system are $r=0.3$, $r_{wx}=0.88$,
$r_{wy}=0.39$, $\sigma_{x}=0.1$ and $\sigma_{y}=0.05$. The initial conditions
for the hard disc granular gas in (b) are $E^{P}_{tot}(0)=10.05$,
$E^{Q}_{tot}(0)=10.04$, $E^{P}_{dif}(0)=2.797$ and $E^{Q}_{dif}(0)=1.979$. The
choice of the other parameters defining the system are $r=0.65$, $m=1$,
$n=0.02$, $E^{st}_{tot}=10.0$ and $E^{st}_{dif}=2.82$. $P$ relaxes to the
steady state faster than $Q$, though it is initially at a larger distance
compared to the final steady state.

###  IV.2 Euclidean Distance

Mean kinetic energy as a measure of distance has the possible issue that the
trajectories may appear to cross in this one-dimensional projection, though in
the two dimensional phase they do not cross. More natural definitions are to
use Euclidean or Manhattan distances. As the initial and the final states are
defined by the pair of variables $\\{E_{tot}(t),E_{dif}(t)\\}$, we can define
an Euclidean measure for the trajectory connecting the two states as

$\displaystyle\mathcal{L}_{2}(t)=\sqrt{(E_{tot}(t)-E^{st}_{tot})^{2}+(E_{dif}(t)-E^{st}_{dif})^{2}}.$
(33)

Here, the initially “hotter” system or equivalently the system which is
initially farther from the final steady state has an initially larger
$\mathcal{L}_{2}$ compared to the “colder” system. For this measure, we define
the Mpemba effect as follows. Let us consider two systems $P$ and $Q$ such
that $P$ is initially at a larger distance from the final steady state
compared to $Q$, i.e., $\mathcal{L}^{P}_{2}(0)>\mathcal{L}^{Q}_{2}(0)$. Here,
the systems $P$ and $Q$ are identical in all respect except for the pair of
driving strengths $(\sigma_{x},\sigma_{y})$ that is required to prepare the
systems in their initial steady states. Then both the systems are quenched to
a common steady state by applying a same pair of driving strengths. In this
case, the Mpemba effect is said to exist if the two trajectories for the
systems $P$ and $Q$ quantified in terms of $\mathcal{L}^{P}_{2}(t)$ and
$\mathcal{L}^{Q}_{2}(t)$ cross each other at some finite time $t=\tau$ at
which

$\mathcal{L}^{P}_{2}(\tau)=\mathcal{L}^{Q}_{2}(\tau).$ (34)

To obtain the value of $\tau$, we equate the Euclidean measures for the
systems $P$ and $Q$ at $t=\tau$ as

$\displaystyle\sqrt{(E^{P}_{tot}(\tau)-E^{st}_{tot})^{2}+(E^{P}_{dif}(\tau)-E^{st}_{dif})^{2}}$
$\displaystyle=\sqrt{(E^{Q}_{tot}(\tau)-E^{st}_{tot})^{2}+(E^{Q}_{dif}(\tau)-E^{st}_{dif})^{2}}.$
(35)

For the inelastic Maxwell gas, using Eq. (21) we obtain the crossing times as

$\displaystyle\tau^{\pm}=\frac{1}{(\lambda_{+}-\lambda_{-})}\ln\left[\frac{2a}{d\pm\sqrt{d^{2}-4ab}}\right],$
(36)

where,

$\displaystyle\begin{split}a&=\big{[}(K^{P}_{+})^{2}+(L^{P}_{+})^{2}-(K^{Q}_{+})^{2}-(L^{Q}_{+})^{2}\big{]},\\\
b&=\big{[}(K^{P}_{-})^{2}+(L^{P}_{-})^{2}-(K^{Q}_{-})^{2}-(L^{Q}_{-})^{2}\big{]},\\\
d&=2\big{[}K^{Q}_{+}K^{Q}_{-}+L^{Q}_{+}L^{Q}_{-}-K^{P}_{+}K^{P}_{-}-L^{P}_{+}L^{P}_{-}\big{]}.\end{split}$
(37)

A similar expression for the crossing time [see Eq. (36)] is also obtained for
the case of hard disc granular gas for which the constants $\lambda_{\pm}$,
$a$, $b$, and $d$ are computed using Eq. (27).

Note that there are two possibilities for the crossing time in Eq. (36).
Figure 2 illustrates such a scenario where both the crossings are present.
However, the presence of two crossings will eventually lead to no Mpemba
effect. Thus, we are interested in only those cases or the initial conditions
where there exists only one crossing of the trajectories of the initial states
leading to the final steady state. For the Mpemba effect to be present, we
require that $\tau^{+}>0$ or $\tau^{-}>0$, but not both positive. Since
$\lambda_{+}>\lambda_{-}$, the argument of logarithm in Eq. (36) should be
greater than one for $+(-)$ and less than one for $-(+)$. Figure 3 illustrates
the existence of the Mpemba effect where the trajectories of the initial
states leading to final steady state are defined in terms of the Euclidean
measure and their crossing time is given by Eq. (36).

Figure 2: The time evolution of the anisotropically driven inelastic Maxwell
gas in terms of the Euclidean measure, $\mathcal{L}_{2}(t)$ for two systems
$P$ and $Q$ with initial conditions $\mathcal{L}^{P}_{2}(0)=6.34$ and
$\mathcal{L}^{Q}_{2}(0)=6.15$, shows two crossings as illustrated in (a) and
(b) for the different times. The multiple crossing times are obtained using
Eq. (36). The choice of the other parameters defining the system are $r=0.1$,
$r_{wx}=0.95$, $r_{wy}=0.39$, $\sigma_{x}=1.6$ and $\sigma_{y}=1.1$.  Figure
3: The time evolution of anisotropically driven (a) inelastic Maxwell and (b)
hard disc granular gas is illustrated in terms of Euclidean measure,
$\mathcal{L}_{2}(t)$ for two systems $P$ and $Q$. The initial conditions for
the inelastic Maxwell gas in (a) are $\mathcal{L}^{P}_{2}(0)=3.13$ and
$\mathcal{L}^{Q}_{2}(0)=2.79$. The choice of the other parameters defining the
system are $r=0.3$, $r_{wx}=0.95$, $r_{wy}=0.39$, $\sigma_{x}=1.6$ and
$\sigma_{y}=1.1$. The initial conditions for the hard disc granular gas in (b)
are $\mathcal{L}^{P}_{2}(0)=2.02$ and $\mathcal{L}^{Q}_{2}(0)=1.186$. The
choice of the other parameters defining the system are $r=0.3$, $m=1$ and
$n=0.02$. $P$ relaxes to the steady state faster than $Q$, though its initial
Euclidean distance from the final steady state is larger.

###  IV.3 Manhattan Distance

Similar to the Euclidean measure, we can define another measure for the time
evolution of the trajectory connecting the initial and the final states in
terms of Manhattan distance. In the plane of $\\{E_{tot}(t),E_{dif}(t)\\}$,
the Manhattan measure is defined as

$\displaystyle\mathcal{L}_{1}(t)=|E_{tot}(t)-E^{st}_{tot}|+|E_{dif}(t)-E^{st}_{dif}|.$
(38)

Here, the system having larger $\mathcal{L}_{1}(0)$ at time $t=0$ is termed
the “hotter” system while compared to the initially “colder” system having a
comparatively smaller $\mathcal{L}_{1}(0)$. For this measure, we define the
Mpemba effect as follows. Let us consider two systems $P$ and $Q$ such that
$P$ is initially at a larger distance from the final steady state compared to
$Q$ in terms of Manhattan measure, i.e.,
$\mathcal{L}^{P}_{1}(0)>\mathcal{L}^{Q}_{1}(0)$. Here, the systems $P$ and $Q$
are identical except for the pair of driving strengths ($\sigma_{x}$,
$\sigma_{y}$). Both the systems are then subjected to same pair of driving
strengths. As a result, both the systems are driven to a common steady state.
Then the Mpemba effect is said to exist if the two trajectories for the
systems $P$ and $Q$ quantified in terms of $\mathcal{L}^{P}_{1}(t)$ and
$\mathcal{L}^{Q}_{1}(t)$ cross each other at some finite time $t=\tau$ at
which

$\mathcal{L}^{P}_{1}(\tau)=\mathcal{L}^{Q}_{1}(\tau),$ (39)

or equivalently,

$\displaystyle|(E^{P}_{tot}(\tau)-E^{st}_{tot})|+|(E^{P}_{dif}(\tau)-E^{st}_{dif})|$
$\displaystyle=|(E^{Q}_{tot}(\tau)-E^{st}_{tot})|+|(E^{Q}_{dif}(\tau)-E^{st}_{dif})|.$
(40)

Now, using Eq. (21) for the case of inelastic Maxwell gas and solving Eq.
(40), we obtain eight different solutions for the crossing time as

$\displaystyle\begin{split}\tau_{c,1}&=\frac{1}{\gamma}\ln\Big{(}\frac{-K^{P}_{+}+K^{Q}_{+}+L^{P}_{+}+L^{Q}_{+}}{K^{P}_{-}-K^{Q}_{-}-L^{P}_{-}-L^{Q}_{-}}\Big{)},\\\
\tau_{c,2}&=\frac{1}{\gamma}\ln\Big{(}\frac{-K^{P}_{+}-K^{Q}_{+}-L^{P}_{+}+L^{Q}_{+}}{K^{P}_{-}+K^{Q}_{-}+L^{P}_{-}-L^{Q}_{-}}\Big{)},\\\
\tau_{c,3}&=\frac{1}{\gamma}\ln\Big{(}\frac{-K^{P}_{+}-K^{Q}_{+}+L^{P}_{+}-L^{Q}_{+}}{K^{P}_{-}+K^{Q}_{-}-L^{P}_{-}+L^{Q}_{-}}\Big{)},\\\
\tau_{c,4}&=\frac{1}{\gamma}\ln\Big{(}\frac{-K^{P}_{+}+K^{Q}_{+}-L^{P}_{+}-L^{Q}_{+}}{K^{P}_{-}-K^{Q}_{-}+L^{P}_{-}+L^{Q}_{-}}\Big{)},\\\
\tau_{c,5}&=\frac{1}{\gamma}\ln\Big{(}\frac{-K^{P}_{+}-K^{Q}_{+}+L^{P}_{+}+L^{Q}_{+}}{K^{P}_{-}+K^{Q}_{-}-L^{P}_{-}-L^{Q}_{-}}\Big{)},\\\
\tau_{c,6}&=\frac{1}{\gamma}\ln\Big{(}\frac{-K^{P}_{+}+K^{Q}_{+}-L^{P}_{+}+L^{Q}_{+}}{K^{P}_{-}-K^{Q}_{-}+L^{P}_{-}-L^{Q}_{-}}\Big{)},\\\
\tau_{c,7}&=\frac{1}{\gamma}\ln\Big{(}\frac{-K^{P}_{+}+K^{Q}_{+}+L^{P}_{+}-L^{Q}_{+}}{K^{P}_{-}-K^{Q}_{-}-L^{P}_{-}+L^{Q}_{-}}\Big{)},\\\
\tau_{c,8}&=\frac{1}{\gamma}\ln\Big{(}\frac{-(K^{P}_{+}+K^{Q}_{+}+L^{P}_{+}+L^{Q}_{+})}{K^{P}_{-}+K^{Q}_{-}+L^{P}_{-}+L^{Q}_{-}}\Big{)}.\end{split}$
(41)

For the Mpemba effect to be present, there should be odd number of crossings.
Thus, we look for those initial conditions in the phase space that lead to odd
number of solutions to the Eq. (40). Figure 4 illustrates the existence of the
Mpemba effect where the trajectories of the initial states leading to final
steady state are defined in terms of the Manhattan measure.

Figure 4: The time evolution of anisotropically driven (a) inelastic Maxwell
and (b) hard disc granular gas is illustrated in terms of Manhattan measure,
$\mathcal{L}_{1}(t)$ for two systems $P$ and $Q$. The initial conditions for
the inelastic Maxwell gas in (a) are $\mathcal{L}^{P}_{1}(0)=7.54$ and
$\mathcal{L}^{Q}_{1}(0)=7.06$. The choice of the other parameters defining the
system are $r=0.9$, $r_{wx}=0.88$, $r_{wy}=0.3$, $\sigma_{x}=1.6$ and
$\sigma_{y}=1.0$. The initial conditions for the hard disc granular gas in (b)
are $\mathcal{L}^{P}_{1}(0)=2.12$ and $\mathcal{L}^{Q}_{1}(0)=1.62$. The
choice of the other parameters defining the system are $r=0.3$, $m=1$ and
$n=0.02$. $P$ relaxes to the steady state faster than $Q$, though its initial
Manhattan distance from the final steady state is larger.

###  IV.4 KL Divergence

In this section, we define the Mpemba effect in terms of an information
theoretic quantity known as Kullback-Leibler (KL) divergence. The measure is
defined as

$\displaystyle D_{KL}(t)=\int
d\boldsymbol{v}P(\boldsymbol{v},t)\ln\Big{(}\frac{P(\boldsymbol{v},t)}{P^{st}(\boldsymbol{v})}\Big{)},$
(42)

where $P(\boldsymbol{v},t)$ is the instantaneous velocity distribution of the
particles at any time $t$ and $P^{st}(\boldsymbol{v})$ is the final steady
state velocity distribution function. The above quantity is not a true measure
of geometric distance as it is not symmetric between two given distribution
functions. But it is a good candidate for the study of relaxation dynamics of
an arbitrary initial state to a given reference steady state for the following
reasons: (a) it is a monotonically non-increasing function of time as has been
shown in earlier studies [47, 48, 49], (b) it provides information regarding
the measure of deviation and its temporal evolution of any two initial states
from the final reference state.

Although being a good candidate for the study of relaxation dynamics, the
above measure is not much fruitful in the context of granular gases. It is
because Eq. (42) requires to have information regarding the time evolution of
the velocity distribution function at any time $t$ which is not analytically
feasible in the case of granular gases. The previous theoretical studies on
velocity distribution of driven granular gases attempts to find the velocity
distribution of the steady state for different contexts [50, 51, 43, 44, 52].
But in all the cases, it was not possible to derive the form of instantaneous
velocity distribution at an instant of time $t$ even for the simplest model of
inelastic Maxwell gas.

As a result, we use numerical methods to compute the KL divergence. We briefly
discuss the procedure for the numerical computation of KL divergence for the
case of inelastic Maxwell gas. Starting from a random configuration of
velocities for $N$ particles, the system is evolved to an initial steady state
corresponding to an initial pair of driving strengths
$(\sigma_{x},\sigma_{y})$. The pair of driving strengths acting on the
particles is then changed to those of the desired final steady state. The
system evolves in time as follows. At each time step, either two particles
collide or a particle is driven, based on the corresponding rates. At each
time step we measure the time evolution of the velocity distribution,
$P(\boldsymbol{v},t)$ of the $N$ particles. The distributions are averaged
over $10^{4}$ realizations. The final steady state velocity distribution is
measured separately and is averaged over $10^{5}$ different realizations.
Thus, having information about $P(\boldsymbol{v},t)$ and
$P^{st}(\boldsymbol{v})$, we can then numerically determine the time evolution
of KL divergence using Eq. (42).

We now discuss the notion of a “hotter” and a “colder” system using the
measure of KL divergence when two systems are quenched to a final steady
state. The system having larger divergence at $t=0$ is referred to as the
“hotter” system as it is farther off from the final steady state whereas the
system having comparatively smaller divergence is referred to as the “colder”
system. Similar to the previous definitions of the Mpemba effect for the other
measures, if the hotter system equilibrates faster than the colder system to
the final steady state, the Mpemba effect is said to exist.

Figure 5: The time evolution of anisotropically driven (a) inelastic Maxwell
and (b) hard disc granular gas is illustrated in terms of KL divergence
measure, $D_{KL}(t)$ for two systems $P$ and $Q$. The initial conditions for
the inelastic Maxwell gas in (a) are $D^{P}_{KL}(0)=0.82$ and
$D^{Q}_{KL}(0)=0.64$. The choice of the other parameters defining the system
are $r=0.3$, $r_{wx}=0.88$, $r_{wy}=0.39$, $\sigma_{x}=1.5$ and
$\sigma_{y}=0.9$. The initial conditions for the hard disc granular gas in (b)
are $D^{P}_{KL}(0)=0.01$ and $D^{Q}_{KL}(0)=0.0078$. The choice of the other
parameters defining the system are $r=0.65$, $m=1$, $n=0.02$,
$\sigma_{x}=6.62$ and $\sigma_{y}=0.026$. $P$ relaxes to the steady state
faster than $Q$, though its initial KL divergence with respect to the final
steady state is larger.

Figure 5 illustrates the existence of the Mpemba effect where the trajectories
of the initial states leading to final steady state are defined in terms of
the KL divergence.

##  V Comparison between various measures

We now compare the effect of different distance measures on the Mpemba effect.
For a given initial condition in terms of intrinsic system parameters and the
pair of driving strengths, we look for the existence of the Mpemba effect
using the different distance measures.

We first show that for the same initial condition, different measures can give
different results, as illustrated in Fig. 6 for the inelastic Maxwell gas. In
the example shown, the Mpemba effect is absent for distance measure
corresponding to the total energy distance [see Fig. 6(a)] which it is present
for the other three measures [see Fig. 6(b)–(d)]. In addition to the non-
uniqueness of the Mpemba effect, even the notion of “hot” and “cold” system in
terms of distance from the final steady state for a pair of initial conditions
is not unique among the various definitions. What is initially hotter (shown
in red) in Fig. 6(a)–(c) is initially colder (shown in blue) in Fig. 6(d)

Figure 6: The time evolution of the two systems $P$ and $Q$ for an
anisotropically driven inelastic Maxwell gas is illustrated for the following
measures: (a) total energy, (b) Euclidean measure, (c) Manhattan measure and
for (d) KL divergence. The initial conditions for the various measures are
identical and are given in terms of the driving strengths for systems $P$ and
$Q$ as $(\sigma^{P}_{x}=1.9,\sigma^{P}_{y}=1.2)$ and
$(\sigma^{Q}_{x}=1.55,\sigma^{Q}_{y}=2.0)$. The choice of the other parameters
defining the system are $r=0.3$, $r_{wx}=0.95$, $r_{wy}=0.39$,
$\sigma_{x}=1.6$ and $\sigma_{y}=1.1$. The existence of the Mpemba effect and
the notion of “hot” and “cold” system in terms of distance from the final
steady state for a given pair of initial conditions is not unique among the
various measures.

Given that there is non-uniqueness among the definitions, we now check whether
all definitions show the Mpemba effect as well as whether there is any overlap
in the phase space regions that corresponds to the Mpemba effect for the
different measures. The phase space is labeled by the free parameters: final
steady states in terms of ($E^{st}_{tot},~{}E^{st}_{dif}$), initial steady
states of $P$ and $Q$, i.e., ($E^{P}_{tot},~{}E^{P}_{dif}$) and
($E^{Q}_{tot},~{}E^{Q}_{dif}$) for a fixed value of other intrinsic parameters
such as $r$, $r_{wx}$ and $r_{wy}$. However, we notice that the equations for
the correlation functions [see Eqs. (21) and (25)] are linear equations in the
differences in energies $\delta E_{tot}=E_{tot}-E^{st}_{tot}$, $\delta
E_{dif}=E_{dif}-E^{st}_{dif}$ for both $P$ and $Q$, giving four variables. We
also note that the equations for the correlation function, and hence the
crossing times, are invariant if all the $\delta E$s are scaled by the same
factor, thus reducing the number of factors by one. To make it a two
dimensional phase diagram, we fix the initial values of $\delta E^{P}_{tot}$
and $\delta E^{P}_{dif}$ for $P$, and determine the phase diagram in terms of
($\delta E^{Q}_{dif}/\delta E^{P}_{tot}$) and ($\delta E^{Q}_{tot}/\delta
E^{P}_{tot})$.

Figure 7: The phase diagram in the ($\delta E^{Q}_{dif}/\delta
E^{P}_{tot}$)-($\delta E^{Q}_{tot}/\delta E^{P}_{tot})$ plane shows the
existence of the Mpemba effect in the driven inelastic Maxwell gas for the use
of different distance measures: (a) total energy, (b) Manhattan, (c)
Euclidean, and (d) KL-divergence. The red (green) region corresponds to
absence (presence) of Mpemba effect, while the white regions are not
accessible. The choice of the parameters defining the system are $r=0.4$,
$r_{wx}=0.44$, $r_{wy}=0.95$, $\delta E^{P}_{tot}=1.00$ and $\delta
E^{P}_{dif}/\delta E^{P}_{tot}=0.53$, which are kept constant across all the
various measures.  Figure 8:  The phase diagram in the ($\delta
E^{Q}_{dif}/\delta E^{P}_{tot}$)-($\delta E^{Q}_{tot}/\delta E^{P}_{tot})$
plane shows the existence of the Mpemba effect in the driven hard disc
granular gas for the use of different distance measures: (a) total energy, (b)
Manhattan, and (c) Euclidean. The red (green) region corresponds to absence
(presence) of Mpemba effect, while the white regions are not accessible. The
choice of the parameters defining the system are $n=0.02$, $\sigma=1$, $m=1$,
$r=0.1$, $\delta E^{P}_{tot}=1.00$ and $\delta E^{P}_{dif}/\delta
E^{P}_{tot}=0.11$, which are kept constant across all the various measures.

The phase diagrams for the measures: total energy, Manhattan and Euclidean
distance, are obtained using the exact criteria derived for the existence of
the Mpemba effect, both for the inelastic Maxwell gas and hard disc granular
gas (linearized analysis). However, for the case of KL divergence, the phase
diagram is obtained using discrete set of points sampled from the phase space
and checked individually for the existence of the Mpemba effect through the
analysis of evolution of the trajectories of systems $P$ and $Q$. The KL
divergence is used only for the Maxwell gas as in the case of granular gas, it
is not possible to make a direct comparison with the linearized theories.
Figures 7 and 8 show the phase diagram of the Mpemba effect in terms of
various measures for the case of inelastic Maxwell gas and hard disc granular
gas respectively. The colored region corresponds to the allowed parameter
space, while the green (red) regions correspond to the Mpemba effect being
present (absent). It is clear that all the different distance measures show
the Mpemba effect. However, the phase boundary of the initial conditions that
show the Mpemba effect is very different for the various measures. We conclude
that the usual measures used for driven granular gases lead to a non-universal
definition of the Mpemba effect.

##  VI Summary and conclusion

To summarize, we studied the Mpemba effect in the driven inelastic Maxwell and
the driven granular gases using different definitions of distances between
points in phase space. These definitions were motivated by the choices that
have been used earlier in the literature for establishing the presence of the
Mpemba effect in granular systems. We studied the Mpemba effect in terms of
the following measures: total kinetic energy (as has been used earlier),
Euclidean measure ($\mathcal{L}_{2}$), Manhattan measure ($\mathcal{L}_{1}$)
and KL divergence. While the first three are based on average kinetic energy
in different directions, KL divergence is based on the probability
distribution of velocity. The analysis was performed for the anisotropically
driven gases [39, 40].

We derived the criteria for the existence of the Mpemba effect with the
various measures and showed the existence of the Mpemba effect for all the
different choices of distance measures. However, the phase diagrams for the
different measures are not the same with the presence or absence of the Mpemba
effect for a particular phase space point depending on the measure, resulting
in a non-unique definition of the Mpemba effect. Moreover, the notion of ‘hot’
and ‘cold’ initial states is different for the different choices of distance
measures from the final steady state.

In a Markov process, the relaxation to the steady state at large times is
dominated by the first “excited” state of the transition matrix, i.e.,
$P(C,t)=P_{1}(C)+a_{2}\exp(-t/\tau)P_{2}(C)+\ldots$, where $P(C,t)$ is the
probability distribution at time $t$, $P_{1}$ is the steady state and $P_{2}$
is the first excited state. The quantities $\tau$, $P_{1}$ and $P_{2}$ depend
only on the transition matrix, and independent of the initial states, while
the pre-factor $a_{2}$ is determined by the overlap of the initial steady
state with the first excited state. The system with smaller value of $|a_{2}|$
will equilibrate faster [22]. The same is true for the time evolution of an
arbitrary observable $f(C,t)$ as it also follows a similar structure as
$P(C,t)$. Therefore, if system $P$ approaches equilibrium faster according to
KL-divergence, then we will find similar relaxation behaviour with total
energy as the distance measure. However, what is cold and what is hot may not
be the same with both measures (see Fig. 6(a) and 6(d) for an example), thus
giving contradictory results for the Mpemba effect. Since it is difficult to
come up with an observable $f(C,t)$ that has the same characteristics as KL-
divergence, we generalised the total variation distance
$\left|P(C,t)-P_{1}(C)\right|$ to the distances $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$ used in this paper, based on the mean energies.

The ambiguity of characterising the Mpemba effect in granular systems is in
contrast to the unique identification of the Mpemba effect for the case of
single particle Markov systems, irrespective of the choice of distance
measure. This is due to the existence of a common criteria across all the
measures which is given in terms of the coefficient $a_{2}$ [22] associated
with the second eigenfunction in the eigenspectrum analysis of probability
distribution function. However, it is very challenging to calculate $a_{2}$
for an interacting many particle system like the granular system considered in
this paper. The natural choices are observables like kinetic energy, rather
than probability distributions, which are easier to track both in experiments
and calculations. However, these choices for distance measures may not
necessarily decrease monotonically with time, a requirement that was put
forward in Ref. [22]. Based on the results in this paper, the characterization
of the Mpemba effect with such ad-hoc distance measures should be done with
caution.

##  Appendix A  Time evolution of two point correlations for inelastic Maxwell
gas

In this appendix, we solve for the time evolution of the various two-point
correlations for the inelastic Maxwell gas. The two-point correlations that we
are interested in are discussed in Eq. (7) whose time evolution is given by
$\frac{d\boldsymbol{\tilde{\Sigma}}(t)}{dt}=\boldsymbol{\tilde{R}}\boldsymbol{\tilde{\Sigma}}(t)+\boldsymbol{\tilde{D}}$
where the matrices $\boldsymbol{\tilde{\Sigma}}$ and $\boldsymbol{\tilde{D}}$
are given in Eq. (9). The form of the matrix $\boldsymbol{\tilde{R}}$ is given
by

$\displaystyle\tilde{\boldsymbol{R}}\\!=$
$\displaystyle\\!\left[\begin{array}[]{cccccc}\frac{A_{1}(N-1)}{N}+A^{xx}_{4}&\frac{A_{2}(N-1)}{N}&0&\frac{-A_{1}(N-1)}{N}&\frac{-A_{2}(N-1)}{N}&0\\\
\frac{A_{2}(N-1)}{N}&\frac{A_{1}(N-1)}{N}+A^{yy}_{4}&0&\frac{-A_{2}(N-1)}{N}&\frac{-A_{1}(N-1)}{N}&0\\\
0&0&\frac{-A_{3}(N-1)}{N}+A^{xy}_{4}&0&0&\frac{A_{3}(N-1)}{N}\\\
-\frac{A_{1}}{N}&-\frac{A_{2}}{N}&0&\frac{A_{1}}{N}+2A^{x}_{5}&\frac{A_{2}}{N}&0\\\
-\frac{A_{2}}{N}&-\frac{A_{1}}{N}&0&\frac{A_{2}}{N}&\frac{A_{1}}{N}+2A^{y}_{5}&0\\\
0&0&\frac{A_{3}}{N}&0&0&-\frac{A_{3}}{N}+A^{x}_{5}+A^{y}_{5}\end{array}\right],$
(49)

where the constants $A_{1},A_{2},A_{3},A^{ij}_{4},A^{i}_{5}$, where
$i,j\in(x,y)$ are given by

$\displaystyle\begin{split}&A_{1}=\frac{3}{4}\lambda_{c}\alpha^{2}-\lambda_{c}\alpha,\quad
A_{2}=\frac{\lambda_{c}\alpha^{2}}{4},\\\
&A_{3}=\lambda_{c}\alpha(1-\frac{\alpha}{2}),\quad
A^{ij}_{4}=-\lambda_{d}(1-r_{wi}r_{wj}),\\\
&A^{i}_{5}=-\lambda_{d}(1+r_{wi}),\quad\textnormal{where}~{}i,j\in(x,y).\end{split}$
(50)

In the thermodynamic limit $N\rightarrow\infty$, the matrix
$\tilde{\boldsymbol{R}}$ simplifies to

$\displaystyle\boldsymbol{\tilde{R}}\\!=$
$\displaystyle\\!\left[\begin{array}[]{cccccc}A_{1}+A^{xx}_{4}&A_{2}&0&-A_{1}&-A_{2}&0\\\
A_{2}&A_{1}+A^{yy}_{4}&0&-A_{2}&-A_{1}&0\\\ 0&0&-A_{3}+A^{xy}_{4}&0&0&A_{3}\\\
0&0&0&2A^{x}_{5}&0&0\\\ 0&0&0&0&2A^{y}_{5}&0\\\
0&0&0&0&0&A^{x}_{5}+A^{y}_{5}\end{array}\right].$ (57)

In that case, the non-zero two-point correlations that exist are $E_{x}$ and
$E_{y}$ whose time evolution is given by Eq. (10). In the new set of variables
$E_{tot}=E_{x}+E_{y}$ and $E_{dif}=E_{x}-E_{y}$, the time evolutions reduces
to Eq. (16) whose solution is given by Eq. (21) where the coefficients
$K_{+},K_{-},L_{+}$ and $L_{-}$ along with $E^{st}_{tot}$ and $E^{st}_{dif}$
are given by

$\displaystyle
K_{+}=\frac{1}{\gamma}\Big{[}-(\lambda_{-}-\chi_{11})E_{tot}(0)+\chi_{12}E_{dif}(0)-\frac{\lambda_{d}}{\lambda_{+}}\big{[}\big{(}\chi_{12}-(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{y}\big{]}\Big{]},$
$\displaystyle
K_{-}=\frac{1}{\gamma}\Big{[}(\lambda_{+}-\chi_{11})E_{tot}(0)-\chi_{12}E_{dif}(0)+\frac{\lambda_{d}}{\lambda_{-}}\big{[}\big{(}\chi_{12}-(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{y}\big{]}\Big{]},$
$\displaystyle
E^{st}_{tot}=\frac{\lambda_{d}}{\gamma}\Big{[}\frac{\big{(}\chi_{12}-(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{y}}{\lambda_{+}}-\frac{\big{(}\chi_{12}-(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{y}}{\lambda_{-}}\Big{]},$
$\displaystyle
L_{+}=\frac{1}{\gamma}\Big{[}-\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}E_{tot}(0)+(\lambda_{+}-\chi_{11})E_{dif}(0)-\frac{\lambda_{d}}{\lambda_{+}\chi_{12}}\big{[}(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})(\sigma^{2}_{x}-\sigma^{2}_{y})\big{]}\Big{]},$
$\displaystyle
L_{-}=\frac{1}{\gamma}\Big{[}\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}E_{tot}(0)-(\lambda_{+}-\chi_{11})E_{dif}(0)+\frac{\lambda_{d}}{\lambda_{-}\chi_{12}}\big{[}(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})(\sigma^{2}_{x}-\sigma^{2}_{y})\big{]}\Big{]},$
$\displaystyle
E^{st}_{dif}=\frac{\lambda_{d}}{\chi_{12}\gamma}\Big{[}(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})(\sigma^{2}_{x}-\sigma^{2}_{y})\big{(}\frac{1}{\lambda_{+}}-\frac{1}{\lambda_{-}}\big{)}\Big{]},$
$\displaystyle\gamma=\lambda_{+}-\lambda_{-}.$ (58)

##  Appendix B  Hard disc granular gas

###  B.1  Form of the collision integral

The time evolution of the velocity distribution function $f(\boldsymbol{v},t)$
for an anisotropically driven hard disc granular is described using the
Enskog-Boltzmann equation [45]

$\displaystyle\frac{\partial}{\partial t}f(\boldsymbol{{v}},t)=\chi
I(f,f)+\Big{(}\frac{\sigma^{2}_{x}}{2}\frac{\partial^{2}}{\partial
v^{2}_{x}}+\frac{\sigma^{2}_{y}}{2}\frac{\partial^{2}}{\partial
v^{2}_{y}}\Big{)}f(\boldsymbol{{v}},t),$ (59)

where $\chi$ [46] is the pair correlation function, $I(f,f)$ is the collision
integral which is given by

$\displaystyle I(f,f)=\sigma\int d\boldsymbol{{v}}_{2}\int
d\boldsymbol{{e}}\Theta(-\boldsymbol{{v}}_{12}.\boldsymbol{{e}})$
$\displaystyle\left|\boldsymbol{{v}}_{12}.\boldsymbol{{e}}\right|\Big{[}\frac{1}{r^{2}}f(\boldsymbol{v}^{{}^{\prime\prime}}_{1},t)f(\boldsymbol{v}^{{}^{\prime\prime}}_{2},t)-f(\boldsymbol{{v}}_{1},t)f(\boldsymbol{{v}}_{2},t)\Big{]},$
(60)

where
$(\boldsymbol{v}^{{}^{\prime\prime}}_{1},\boldsymbol{v}^{{}^{\prime\prime}}_{2})$
are the pre-collision velocities that give
$(\boldsymbol{v}_{1},\boldsymbol{v}_{2})$ upon collision.

###  B.2  Time evolution of $E_{tot}$ and $E_{dif}$

In this appendix, we write the explicit expressions for the time evolution of
$E_{tot}$ and $E_{dif}$ for the hard disc granular gas. The time evolutions of
$E_{tot}(t)$ and $E_{dif}(t)$ are written in a compact as given in Eq. (24)
where the functions $\mathcal{F}(E_{tot},E_{dif})$ and
$\mathcal{G}(E_{tot},E_{dif})$ are given by

$\displaystyle\mathcal{F}(E^{st}_{tot},E^{st}_{dif})={\frac{m(\xi^{2}_{0x}+\xi^{2}_{0y})}{\nu_{0}}}+\frac{n\chi\sigma(1+r)\sqrt{E^{st}_{tot}-E^{st}_{dif}}}{15\nu_{0}\sqrt{2\pi
m}E^{st}_{dif}}\times\Bigg{[}\Big{(}(3r-7)T^{2}_{tot}+4(7r-3)E^{st}_{tot}E^{st}_{dif}+3(3r-7)T^{2}_{dif}\Big{)}$
$\displaystyle{\textit{E}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}-(3r-7)(E^{st}_{tot}-3E^{st}_{dif})\sqrt{T^{2}_{tot}-T^{2}_{dif}}\textit{{E}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}-(E^{st}_{tot}+E^{st}_{dif})\Big{(}(3r-7)E^{st}_{tot}+(7r-3)E^{st}_{dif}\Big{)}$
$\displaystyle\textit{{K}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}+(3r-7)\sqrt{E^{st}_{tot}+E^{st}_{dif}}(E^{st}_{tot}-E^{st}_{dif})^{3/2}\textit{{K}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\Bigg{]},$
(61)
$\displaystyle\mathcal{G}(E^{st}_{tot},E^{st}_{dif})={\frac{m(\xi^{2}_{0x}-\xi^{2}_{0y})}{\nu_{0}}}+\frac{n\chi\sigma(1+r)(3r-7)\sqrt{E^{st}_{tot}-E^{st}_{dif}}}{15\nu_{0}\sqrt{2\pi
m}E^{st}_{dif}}\times\Bigg{[}(E^{st}_{tot}+3E^{st}_{dif})(E^{st}_{tot}+E^{st}_{dif})\textit{{E}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}$
$\displaystyle+(E^{st}_{tot}-3E^{st}_{dif})\sqrt{T^{2}_{tot}-T^{2}_{dif}}\textit{{E}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}-(E^{st}_{tot}+E^{st}_{dif})^{2}\textit{{K}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}$
$\displaystyle-\sqrt{E^{st}_{tot}+E^{st}_{dif}}(E^{st}_{tot}-E^{st}_{dif})^{3/2}\textit{{K}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\Bigg{]},$
(62)

where $K(x)$ and $E(x)$ are elliptic integrals of first and second kind
respectively. Here, time $t$ is measured in terms of $1/\nu_{0}$ such that it
is dimensionless where $\nu_{0}$ is the frequency of interparticle collisions
given by

$\displaystyle\nu_{0}=\chi\sigma
n\sqrt{\frac{2(E^{st}_{tot}+E^{st}_{dif})}{\pi
m}}E\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}.$ (63)

However, in order to make the analysis analytically tractable, we linearize
the non-linear time evolution equations in (24) by considering initial states
that are close to final stationary state denoted by $E^{st}_{tot}$ and
$E^{st}_{dif}$. In that case, the linearized time evolutions are written in a
compact form as given in Eq. (25) where the elements of the matrix
$\boldsymbol{\chi}$ are given by

$\displaystyle\begin{split}\chi_{11}=&-\frac{n\chi\sigma(1+r)}{30\nu_{0}\sqrt{2\pi
m}E^{st}_{dif}\sqrt{E^{st}_{tot}-E^{st}_{dif}}}\times\Big{[}((3-47r)(E^{st}_{dif})^{2}+8(9r-1)E^{st}_{dif}E^{st}_{tot}+5(3r-7)(E^{st}_{tot})^{2})\textit{{E}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&+\frac{7-3r}{\sqrt{(E^{st}_{tot})^{2}-(E^{st}_{dif})^{2}}}\Big{(}5(E^{st}_{dif})^{3}+3(E^{st}_{dif})^{2}E^{st}_{tot}-13E^{st}_{dif}(E^{st}_{tot})^{2}+5(E^{st}_{tot})^{3}\Big{)}\textit{{E}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&+\frac{(7-3r)(E^{st}_{tot}-E^{st}_{dif})^{3/2}}{\sqrt{E^{st}_{tot}+E^{st}_{dif}}}\Big{(}E^{st}_{dif}(-4\sqrt{E^{st}_{tot}+E^{st}_{dif}}+1)-E^{st}_{tot}(2\sqrt{E^{st}_{tot}+E^{st}_{dif}}+3)\Big{)}\textit{{K}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&+\Big{(}(13r-17)(E^{st}_{dif})^{2}+2(1-9r)E^{st}_{dif}E^{st}_{tot}+5(7-3r)(E^{st}_{tot})^{2}\Big{)}\textit{{K}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&+\frac{\pi
E^{st}_{dif}}{2(E^{st}_{tot}-E^{st}_{dif})}\Big{(}3(3r-7)(E^{st}_{dif})^{2}+4(7r-3)E^{st}_{dif}E^{st}_{tot}+(3r-7)(E^{st}_{tot})^{2}\Big{)}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi(3r-7)E^{st}_{dif}\sqrt{E^{st}_{tot}-E^{st}_{dif}}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}\Big{(}3(E^{st}_{dif})^{2}-4E^{st}_{dif}E^{st}_{tot}+(E^{st}_{tot})^{2}\Big{)}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi
E^{st}_{dif}(E^{st}_{tot}+E^{st}_{dif})}{2(E^{st}_{tot}-E^{st}_{dif})}\Big{(}(7r-3)E^{st}_{dif}+(3r-7)E^{st}_{tot}\Big{)}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi(3r-7)E^{st}_{dif}(E^{st}_{tot}-E^{st}_{dif})^{5/2}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\Big{]},\end{split}$
(64)
$\displaystyle\begin{split}\chi_{12}=&-\frac{n\chi\sigma(1+r)}{30\nu_{0}\sqrt{2\pi
m}(E^{st}_{dif})^{2}\sqrt{E^{st}_{tot}-E^{st}_{dif}}}\\\
&\times\Big{[}\Big{(}9(7-3r)(E^{st}_{dif})^{3}-10(r+3)(E^{st}_{dif})^{2}E^{st}_{tot}-(7-3r)E^{st}_{dif}(E^{st}_{tot})^{2}+2(7-3r)(E^{st}_{tot})^{3}\Big{)}\textit{{E}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-(7-3r)\frac{(E^{st}_{tot}-E^{st}_{dif})^{3/2}}{\sqrt{E^{st}_{tot}+E^{st}_{dif}}}\Big{(}15(E^{st}_{dif})^{2}+9E^{st}_{dif}E^{st}_{tot}-2(E^{st}_{tot})^{2}\Big{)}\textit{{E}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&-\frac{(7-3r)\sqrt{E^{st}_{tot}-E^{st}_{dif}}}{\sqrt{E^{st}_{tot}+E^{st}_{dif}}}\Big{(}5(E^{st}_{dif})^{3}-4(E^{st}_{dif})^{2}E^{st}_{tot}-3(E^{st}_{tot})^{2}E^{st}_{dif}+2(E^{st}_{tot})^{3}\Big{)}\textit{{K}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&+\Big{(}-(3-7r)(E^{st}_{dif})^{3}+4(1+r)(E^{st}_{dif})^{2}E^{st}_{tot}+5(1-5r)(E^{st}_{tot})^{2}E^{st}_{dif}+2(7-3r)(E^{st}_{tot})^{3}\Big{)}\textit{{K}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi
E^{st}_{dif}E^{st}_{tot}}{2(E^{st}_{tot}-E^{st}_{dif})}\Big{(}3(7-3r)(E^{st}_{dif})^{2}+4(3-7r)E^{st}_{dif}E^{st}_{tot}+(7-3r)(E^{st}_{tot})^{2}\Big{)}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi(7-3r)E^{st}_{dif}E^{st}_{tot}\sqrt{E^{st}_{tot}-E^{st}_{dif}}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}\Big{(}3(E^{st}_{dif})^{2}-4E^{st}_{dif}E^{st}_{tot}+(E^{st}_{tot})^{2}\Big{)}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi
E^{st}_{dif}E^{st}_{tot}(E^{st}_{tot}+E^{st}_{dif})}{2(E^{st}_{tot}-E^{st}_{dif})}\Big{(}(3-7r)E^{st}_{dif}+(7-3r)E^{st}_{tot}\Big{)}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi(7-3r)E^{st}_{dif}E^{st}_{tot}(E^{st}_{tot}-E^{st}_{dif})^{5/2}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\Big{]},\end{split}$
(65)
$\displaystyle\begin{split}\chi_{21}=&\frac{n\chi\sigma(1+r)(7-3r)}{30\nu_{0}\sqrt{2\pi
m}E^{st}_{dif}\sqrt{E^{st}_{tot}-E^{st}_{dif}}}\times\Big{[}\Big{(}-5(E^{st}_{dif})^{2}+8E^{st}_{dif}E^{st}_{tot}+5(E^{st}_{tot})^{2}\Big{)}\textit{{E}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&+\Big{(}\frac{(E^{st}_{tot}+E^{st}_{dif})^{2}}{2}+(3E^{st}_{tot}-E^{st}_{dif})(E^{st}_{tot}-3E^{st}_{dif})\frac{\sqrt{E^{st}_{tot}-E^{st}_{dif}}}{\sqrt{E^{st}_{tot}+E^{st}_{dif}}}\Big{)}\textit{{E}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&-\Big{(}\frac{4}{\sqrt{(E^{st}_{tot})^{2}-(E^{st}_{dif})}}\big{(}4(E^{st}_{dif})^{2}-E^{st}_{dif}E^{st}_{tot}-(E^{st}_{tot})^{2}\big{)}\Big{)}\textit{{E}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&-\frac{(E^{st}_{tot}-E^{st}_{dif})^{3/2}}{\sqrt{E^{st}_{tot}+E^{st}_{dif}}}\Big{(}3E^{st}_{dif}+5E^{st}_{tot}\Big{)}\textit{{K}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}-\Big{(}(5E^{st}_{tot}-3E^{st}_{dif})(E^{st}_{tot}+E^{st}_{dif})\Big{)}\textit{{K}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi
E^{st}_{dif}(E^{st}_{tot}+3E^{st}_{dif})(E^{st}_{tot}+E^{st}_{dif})}{2(E^{st}_{tot}-E^{st}_{dif})}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}-\frac{\pi
E^{st}_{dif}(E^{st}_{tot}+E^{st}_{dif})^{2}}{2(E^{st}_{tot}-E^{st}_{dif})}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&+\frac{\pi
E^{st}_{dif}(E^{st}_{tot}-3E^{st}_{dif})(E^{st}_{tot}-E^{st}_{dif})^{3/2}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}+\frac{\pi
E^{st}_{dif}(E^{st}_{tot}-E^{st}_{dif})^{5/2}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\Big{]},\end{split}$
(66)

and

$\displaystyle\begin{split}\chi_{22}=&\frac{n\chi\sigma(1+r)(7-3r)}{30\nu_{0}\sqrt{2\pi
m}(E^{st}_{dif})^{2}\sqrt{E^{st}_{tot}-E^{st}_{dif}}}\times\Big{[}\Big{(}-9(E^{st}_{dif})^{3}+2(E^{st}_{dif})^{2}E^{st}_{tot}+E^{st}_{dif}(E^{st}_{tot})^{2}-2(E^{st}_{tot})^{3}\Big{)}\textit{{E}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&+\frac{\sqrt{E^{st}_{tot}-E^{st}_{dif}}}{\sqrt{E^{st}_{tot}+E^{st}_{dif}}}\Big{(}9(E^{st}_{dif})^{3}+2(E^{st}_{dif})^{2}E^{st}_{tot}-E^{st}_{dif}(E^{st}_{tot})^{2}-2(E^{st}_{tot})^{3}\Big{)}\textit{{E}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&+\frac{(E^{st}_{tot}-E^{st}_{dif})^{3/2}}{\sqrt{E^{st}_{tot}+E^{st}_{dif}}}\Big{(}3(E^{st}_{dif})^{2}+3E^{st}_{dif}E^{st}_{tot}+2(E^{st}_{tot})^{2}\Big{)}\textit{{K}}\Big{(}\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&+\Big{(}3(E^{st}_{dif})^{2}-3E^{st}_{dif}E^{st}_{tot}+2(E^{st}_{tot})^{2}\Big{)}\textit{{K}}\Big{(}\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&+\frac{\pi
E^{st}_{dif}E^{st}_{tot}(E^{st}_{tot}+E^{st}_{dif})(E^{st}_{tot}+3E^{st}_{dif})}{2(E^{st}_{tot}-E^{st}_{dif})}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi
E^{st}_{dif}E^{st}_{tot}(E^{st}_{tot}-3E^{st}_{dif})(E^{st}_{tot}-E^{st}_{dif})^{3/2}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}{}_{2}F_{1}\Big{(}\frac{1}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\\\
&+\frac{\pi
E^{st}_{dif}E^{st}_{tot}(E^{st}_{tot}+E^{st}_{dif})^{2}}{2(E^{st}_{tot}-E^{st}_{dif})}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{-2E^{st}_{dif}}{E^{st}_{tot}-E^{st}_{dif}}\Big{)}\\\
&-\frac{\pi
E^{st}_{dif}E^{st}_{tot}(E^{st}_{tot}-E^{st}_{dif})^{5/2}}{2(E^{st}_{tot}+E^{st}_{dif})^{3/2}}{}_{2}F_{1}\Big{(}\frac{3}{2},\frac{3}{2};2;\frac{2E^{st}_{dif}}{E^{st}_{tot}+E^{st}_{dif}}\Big{)}\Big{]},\end{split}$
(67)

where ${}_{2}F_{1}(a,b;c;z)^{\prime}s$ are the hypergeometric function.

## References

  * Mpemba and Osborne [1969] E. B. Mpemba and D. G. Osborne, Cool?, Phys. Educat. 4, 172 (1969).
  * Wojciechowski  _et al._ [1988] B. Wojciechowski, I. Owczarek, and G. Bednarz, Freezing of aqueous solutions containing gases, Crystal Research and Technology 23, 843 (1988).
  * Mirabedin and Farhadi [2017] S. M. Mirabedin and F. Farhadi, Numerical investigation of solidification of single droplets with and without evaporation mechanism, Int. J. Refrig. 73, 219 (2017).
  * Vynnycky and Kimura [2015] M. Vynnycky and S. Kimura, Can natural convection alone explain the Mpemba effect?, Int. J. Heat Mass Transf. 80, 243 (2015).
  * Katz [2009] J. I. Katz, When hot water freezes before cold, Am. J. Phys. 77, 27 (2009).
  * Auerbach [1995] D. Auerbach, Supercooling and the Mpemba effect: When hot water freezes quicker than cold, Am. J. Phys. 63, 882 (1995).
  * Chaddah  _et al._ [2010] P. Chaddah, S. Dash, K. Kumar, and A. Banerjee, Overtaking while approaching equilibrium, arXiv preprint arXiv:1011.3598 (2010).
  * Ahn  _et al._ [2016] Y.-H. Ahn, H. Kang, D.-Y. Koh, and H. Lee, Experimental verifications of Mpemba-like behaviors of clathrate hydrates, Korean J. Chem. Eng. 33, 1903 (2016).
  * Hu  _et al._ [2018] C. Hu, J. Li, S. Huang, H. Li, C. Luo, J. Chen, S. Jiang, and L. An, Conformation directed Mpemba effect on polylactide crystallization, Cryst. Growth Des. 18, 5757 (2018).
  * Greaney  _et al._ [2011] P. A. Greaney, G. Lani, G. Cicero, and J. C. Grossman, Mpemba-like behavior in carbon nanotube resonators, Metallurgical and Materials Transactions A 42, 3907 (2011).
  * Kumar and Bechhoefer [2020] A. Kumar and J. Bechhoefer, Exponentially faster cooling in a colloidal system, Nature 584, 64 (2020).
  * Kumar  _et al._ [2022] A. Kumar, R. Chétrite, and J. Bechhoefer, Anomalous heating in a colloidal system, Proceedings of the National Academy of Sciences 119, e2118484119 (2022).
  * Bechhoefer  _et al._ [2021] J. Bechhoefer, A. Kumar, and R. Chétrite, A fresh understanding of the Mpemba effect, Nature Reviews Physics 1 (2021).
  * Gal and Raz [2020] A. Gal and O. Raz, Precooling strategy allows exponentially faster heating, Phys. Rev. Lett. 124, 060602 (2020).
  * Klich  _et al._ [2019] I. Klich, O. Raz, O. Hirschberg, and M. Vucelja, Mpemba index and anomalous relaxation, Phys. Rev. X 9, 021060 (2019).
  * Klich and Vucelja [2018] I. Klich and M. Vucelja, Solution of the Metropolis dynamics on a complete graph with application to the Markov chain Mpemba effect, arXiv preprint arXiv:1812.11962 (2018).
  * Das and Vadakkayil [2021] S. K. Das and N. Vadakkayil, Should a hotter paramagnet transform quicker to a ferromagnet? Monte Carlo simulation results for Ising model, Phys. Chem. Chem. Phys. (2021).
  * González-Adalid Pemartín  _et al._ [2021] I. González-Adalid Pemartín, E. Mompó, A. Lasanta, V. Martín-Mayor, and J. Salas, Slow growth of magnetic domains helps fast evolution routes for out-of-equilibrium dynamics, Phys. Rev. E 104, 044114 (2021).
  * Teza  _et al._ [2021] G. Teza, R. Yaacoby, and O. Raz, Relaxation shortcuts through boundary coupling, arXiv preprint arXiv:2112.10187 (2021).
  * Holtzman and Raz [2022] R. Holtzman and O. Raz, Landau theory for the Mpemba effect through phase transitions, Communications Physics 5, 280 (2022).
  * Zhang and Hou [2022] S. Zhang and J.-X. Hou, Theoretical model for the Mpemba effect through the canonical first-order phase transition, Physical Review E 106, 034131 (2022).
  * Lu and Raz [2017] Z. Lu and O. Raz, Nonequilibrium thermodynamics of the Markovian Mpemba effect and its inverse, Proc. Natl. Acad. Sci. USA 114, 5083 (2017).
  * Van Vu and Hasegawa [2021] T. Van Vu and Y. Hasegawa, Toward relaxation asymmetry: Heating is faster than cooling, Phys. Rev. Research 3, 043160 (2021).
  * Walker and Vucelja [2021] M. R. Walker and M. Vucelja, Anomalous thermal relaxation of Langevin particles in a piecewise-constant potential, Journal of Statistical Mechanics: Theory and Experiment 2021, 113105 (2021).
  * Busiello  _et al._ [2021] D. M. Busiello, D. Gupta, and A. Maritan, Inducing and optimizing Markovian Mpemba effect with stochastic reset, New Journal of Physics 23, 103012 (2021).
  * Lapolla and Godec [2020] A. Lapolla and A. Godec, Faster uphill relaxation in thermodynamically equidistant temperature quenches, Physical Review Letters 125, 110602 (2020).
  * Degünther and Seifert [2022] J. Degünther and U. Seifert, Anomalous relaxation from a non-equilibrium steady state: An isothermal analog of the Mpemba effect, Europhysics Letters 139, 41002 (2022).
  * Schwarzendahl and Löwen [2022] F. J. Schwarzendahl and H. Löwen, Anomalous cooling and overcooling of active colloids, Phys. Rev. Lett. 129, 138002 (2022).
  * Baity-Jesi  _et al._ [2019] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvión, A. Gordillo-Guerrero, D. Iñiguez, A. Lasanta, A. Maiorano, E. Marinari, _et al._ , The Mpemba effect in spin glasses is a persistent memory effect, Proc. Natl. Acad. Sci. USA 116, 15350 (2019).
  * Santos and Prados [2020] A. Santos and A. Prados, Mpemba effect in molecular gases under nonlinear drag, Phys. Fluids 32, 072010 (2020).
  * Gómez González  _et al._ [2021] R. Gómez González, N. Khalil, and V. Garzó, Mpemba-like effect in driven binary mixtures, Physics of Fluids 33, 053301 (2021).
  * Patrón  _et al._ [2021] A. Patrón, B. Sánchez-Rey, and A. Prados, Strong nonexponential relaxation and memory effects in a fluid with nonlinear drag, Phys. Rev. E 104, 064127 (2021).
  * Megías  _et al._ [2022] A. Megías, A. Santos, and A. Prados, Thermal versus entropic Mpemba effect in molecular gases with nonlinear drag, Phys. Rev. E 105, 054140 (2022).
  * Carollo  _et al._ [2021] F. Carollo, A. Lasanta, and I. Lesanovsky, Exponentially accelerated approach to stationarity in Markovian open quantum systems through the Mpemba effect, Phys. Rev. Lett. 127, 060401 (2021).
  * Lasanta  _et al._ [2017] A. Lasanta, F. Vega Reyes, A. Prados, and A. Santos, When the hotter cools more quickly: Mpemba effect in granular fluids, Phys. Rev. Lett. 119, 148001 (2017).
  * Torrente  _et al._ [2019] A. Torrente, M. A. López-Castaño, A. Lasanta, F. V. Reyes, A. Prados, and A. Santos, Large Mpemba-like effect in a gas of inelastic rough hard spheres, Phys. Rev. E 99, 060901 (2019).
  * Mompó  _et al._ [2021] E. Mompó, M. López-Castaño, A. Lasanta, F. Vega Reyes, and A. Torrente, Memory effects in a gas of viscoelastic particles, Physics of Fluids 33, 062005 (2021).
  * Biswas  _et al._ [2020a] A. Biswas, V. Prasad, O. Raz, and R. Rajesh, Mpemba effect in driven granular Maxwell gases, Physical Review E 102, 012906 (2020a).
  * Biswas  _et al._ [2022a] A. Biswas, V. Prasad, and R. Rajesh, Mpemba effect in an anisotropically driven granular gas, Europhysics Letters 136, 46001 (2022a).
  * Biswas  _et al._ [2022b] A. Biswas, V. Prasad, and R. Rajesh, Mpemba effect in anisotropically driven inelastic Maxwell gases, Journal of Statistical Physics 186, 1 (2022b).
  * Megías and Santos [2022] A. Megías and A. Santos, Mpemba-like effect protocol for granular gases of inelastic and rough hard disks, Frontiers in Physics , 739 (2022).
  * Kullback and Leibler [1951] S. Kullback and R. A. Leibler, On information and sufficiency, The annals of mathematical statistics 22, 79 (1951).
  * Prasad and Rajesh [2019] V. V. Prasad and R. Rajesh, Asymptotic behavior of the velocity distribution of driven inelastic gas with scalar velocities: Analytical results, J. Stat. Phys. 176, 1409–1433 (2019).
  * Prasad  _et al._ [2019] V. V. Prasad, D. Das, S. Sabhapandit, and R. Rajesh, Velocity distribution of driven granular gases, J. Stat. Mech.: Theory Exp. 2019 (6), 063201.
  * van Noije and Ernst [1998] T. van Noije and M. Ernst, Velocity distributions in homogeneous granular fluids: the free and the heated case, Granul. Matter 1, 57 (1998).
  * Brilliantov and Pöschel [2010] N. V. Brilliantov and T. Pöschel,  _Kinetic theory of granular gases_ (Oxford University Press, 2010).
  * Marconi  _et al._ [2013] U. M. B. Marconi, A. Puglisi, and A. Vulpiani, About an $H$-theorem for systems with non-conservative interactions, Journal of Statistical Mechanics: Theory and Experiment 2013, P08003 (2013).
  * Plata and Prados [2017] C. A. Plata and A. Prados, Global stability and $H$-theorem in lattice models with non-conservative interactions, Phys. Rev. E 95, 052121 (2017).
  * Megías and Santos [2020] A. Megías and A. Santos, Kullback–Leibler divergence of a freely cooling granular gas, Entropy 22, 1308 (2020).
  * Prasad  _et al._ [2017a] V. V. Prasad, S. Sabhapandit, A. Dhar, and O. Narayan, Driven inelastic Maxwell gas in one dimension, Phys. Rev. E 95, 022115 (2017a).
  * Prasad  _et al._ [2017b] V. V. Prasad, D. Das, S. Sabhapandit, and R. Rajesh, Velocity distribution of a driven inelastic one-component Maxwell gas, Phys. Rev. E 95, 032909 (2017b).
  * Biswas  _et al._ [2020b] A. Biswas, V. Prasad, and R. Rajesh, Asymptotic velocity distribution of a driven one dimensional binary granular Maxwell gas, Journal of Statistical Mechanics: Theory and Experiment 2020, 013202 (2020b).