Non-monotonic Mpemba effect in binary molecular suspensions

The Mpemba effect is a phenomenon in which an initially hotter sample cools sooner. In this paper, we show the emergence of a non-monotonic Mpemba-like effect in a molecular binary mixture immersed in a viscous gas. Namely, a crossover in the temperature evolution when at least one of the samples presents non-monotonic relaxation. The influence of the bath on the dynamics of the particles is modeled via a viscous drag force plus a stochastic Langevin-like term. Each component of the mixture interchanges energy with the bath depending on the mechanical properties of its particles. This discrimination causes the coupling between the time evolution of temperature with that of the partial temperatures of each component. The non-monotonic Mpemba effect -- and its inverse and mixed counterparts -- stems from this coupling. In order to obtain analytical results, the velocity distribution functions of each component are approximated by considering multitemperature Maxwellian distributions. The theoretical results derived from the Enskog kinetic theory show an excellent agreement with direct simulation Monte Carlo (DMSC) data.


Introduction
Under certain conditions, a sample of water at a hotter temperature freezes faster. This counterintuitive phenomenon is known in literature as the Mpemba effect according to E. B. Mpemba [1], who was the first researcher to make rigorous observations. Although several mechanisms have been proposed to explain this effect in water [2], there are still concerns regarding the true nature of it [3].
Leaving aside the origin of the Mpemba effect in water, it is clear that to observe this effect, the time evolution of temperature must be coupled to that of microscopic-or kinetic-variables. For this reason, the Mpemba effect is considered to be a memory effect [4]. This kind of outof-equilibrium processes have also been observed in various systems such as clathrate hydrates [5], polymers [6], and colloids [7]. In these systems, a Mpemba-like effect emerges when two samples at different initial temperatures evolve in such a way that at a given time t c the relaxation curves for the temperatures cross each other. Thus, freezing is avoided and a straightforward explanation can be provided.
However, there remain a lot of variables that have influence on the temperature evolution. As a consequence, setting the initial conditions giving rise to a crossing of the cooling curves can result in an arduous task. For this reason, from an analytical point of view, simple models have been usually considered to gain some insight into the problem. This is particularly the case of kinetic-theory-based models employed to describe the emergence of Mpembalike effect (and its inverse counterpart, namely when the two initial temperatures are below the asymptotic one) in granular [8][9][10][11]-i.e. with inelastic collisions-and molecular gases [12,13]. In the former case, inelasticity of collisions couples the evolution of the (granular) temperature with that of non-hydrodynamic variables-such as the fourth cumulant or kurtosis and/or the rotational-totranslational temperature ratio. This coupling is the reason of the occurrence of the Mpemba effect in granular gases during the relaxation towards the steady state. On the other hand, no inelasticity couples the temperature with any kinetic variable in the case of molecular gases, therefore other mechanisms must be considered. For example, in a recent paper [12], the molecular gas is driven by means of a non-linear drag force plus a stochastic force. The nonlinearity in the velocity dependence of the drag force turns out in the coupling between the temperature T -seen as a measure of the amount of kinetic energy-with the fourth cumulant a 2 . Nevertheless, in concordance with the granular case [14], a 2 is assumed to be small. As a result, initial temperatures must be chosen to be close enough to each other for the emergence of the Mpemba effect.
A stronger Mpemba-like effect has been recently reported in the case of a mixture of molecular gases [13]. In this paper, the mixture is in contact with a thermal reservoir. As usual [15], in the context of the Enskog kinetic theory, the influence of the surrounding fluid on the solid particles is modeled via an external force composed by two terms: (i) a viscous drag force proportional to the (instantaneous) velocity of the solid particles and (ii) a stochastic Langevin-like term. While the first term tries to mimic the loss of energy due to frictional forces through the drag coefficients γ i (i = 1, 2), the second term simulates the energy transmission due to instantaneous and random collisions with the bath. In this work, in accordance with the results obtained in lattice-Boltzmann simulations [16], the drag coefficients γ i are chosen for each species. This differentiation in the interaction with the background gas causes the coupling between the evolution equations of the temperature T (t) and the partial temperatures T i (t) of each component. This coupling allows that, with an appropriate choice of the initial values, the Mpemba effect (and its inverse equivalent) arises in the relaxation towards equilibrium-defined in terms of a Maxwellian distribution at the background temperature T ex . Moreover, the use of the partial temperatures as the control parameter permits more flexibility in the selection of the initial conditions. Namely, temperature curves can cross even when the initial differences are of the same order than the temperature themselves (large Mpemba effect). Furthermore, a Maxwellian approximation was used to estimate the partial production rates ξ i -which measure rate of energy interchanged in collisions among particles of species i with j. Thus, no cumulants are needed and the effect is easily understood.
The main goal of the present work is to extend the results obtained in [13] to those situations where the temperature presents a non-monotonic relaxation through an appropiately choice of the initial conditions. The non-monotonic Mpemba-like effect appears when at least one of the two samples whose temperature curves cross presents non-monotonic relaxation. Namely, when the temperature difference |T (t) − T ex |-of at least one of the samples-increases at the early stages of the evolution but, at a given time, T (t) reaches the equilibrium value T ex . To fulfil this goal, we consider the Enskog kinetic theory in combination with the Fokker-Planck suspension model described above. Since the theoretical predictions derived in this work are based on the use of the Maxwellian approach, a comparison between theory and simulations is convenient to assess the reliability of the theory. Here, we compare our theoretical results against the numerical solution of the kinetic equations by means of the DSMC method conveniently adapted to the suspension model [17,18].

Enskog kinetic theory
Let us consider an ensemble of hard spheres of masses m 1 and m 2 and diameters σ 1 and σ 2 immersed in a viscous gas at temperature T ex . At a kinetic-theory level, all the relevant information of the system is enclosed in the one-particle velocity distribution functions f i (r, v, t) of species i. For moderate densities and homogeneous states, the functions f i obey the set of coupled Enskog kinetic equations [19]: where is the Enskog collision operator and the operator F i represents the influence of the bath on the dynamics of the particles. For low-Reynolds numbers, this force term is usually represented by a drag force plus a Fokker-Planck collision operator [20]. In this way, the Enskog equation reads [21] ∂ (2) where γ i are the drag coefficients and k B is the Boltzmann constant. Here, according to the simulation results in gas-solid flows reported by Yin and Sundaresan [16], the coefficients γ i can be written as γ i = γ 0 R i , where γ 0 = 36η g / ρ(σ 1 + σ 2 ) , ρ = i m i n i is the total mass density, n i is the number density of the component i, and η g is the viscosity of the solvent. Moreover, dimensionless functions R i depend on the mole fraction x i = n i /(n 1 + n 2 ), mass m 1 /m 2 and size σ 1 /σ 2 ratios, and the partial volume fractions φ i = (π/6)n i σ 3 i , related with the total one by φ = φ 1 + φ 2 . More details of this kind of Langevin-like models can be found in Refs. [21,22].
For homogeneous situations, all the relevant information of the binary mixture is provided by the partial temperatures T i (t) of the component i-or, equivalently, by the temperature ratio θ(t) = T 1 (t)/T 2 (t) and the (total) temperature T (t) = x 1 T 1 (t) + x 2 T 2 (t) of the mixture. The partial temperatures are defined as The evolution equations for the temperature ratio θ and the (reduced) temperature T * = T/T ex can be obtained by multiplying both sides of the Enskog equation (2) by m i v 2 and integrating over velocity. They are given by [13] ∂ ∂t * T * = 2 where we have introduced the scaled time t * = 4(n 1 + n 2 )/(σ 1 + σ 2 ) 2 √ 4k B T ex (m 1 + m 2 )t. Here, µ i j = m i /(m i + m j ), χ 12 is the pair correlation function, In order to write Eq 5, the partial production rates have been calculated by approaching the distribution functions f i to their Maxwellian forms: . (8) More details can be found in Refs. [13,23]. If we freely let the mixture evolve, it will asymptotically achieve the equilibrium state where equipartition holds, i.e. T eq 1 = T eq 2 = T eq = T ex . However, during the time evolution towards this final state, partial temperatures are not equal (T 1 (t) T 2 (t)). This symmetry breaking and the fact that γ 1 γ 2 (because R 1 R 2 ) causes the emergence of the Mpemba effect in the mixture. Note that in the case that γ 1 = γ 2 , total temperature evolution decouples from that of the partial temperatures and the Mpemba effect dissapears.

Results and discussion
In this section, we consider temperature crossings where a non-monotonic relaxation is present. For this reason, the non-linearity of the evolution equations 4-5 plays a crucial role. The linear counterpart has been previously analysed in Ref. [13] where a perfect agreement between exact theoretical expressions-for the (scaled) crossing time t * cand DSMC simulations has been found. Unfortunately, unlike the linear case, we are lack of simple analytical expressions in the far-from-equilibrium framework. Thus, more qualitative analysis is required to establish the necessary but no sufficient conditions for the ocurrence of the Mpemba effect.
Let us consider two homogeneous states A and B far away from equilibrium. The time evolution of these states is fully determined by their initial temperatures T * 0,A and T * 0,B and their temperature ratios θ 0,A and θ 0,B . Thus, we can establish a necessary condition-for a crossover in the temperature relaxation-based on the selection of the initial slopes Φ(T * 0,A , θ 0,A ) and Φ(T * 0,B , θ 0,B ) of each state through their dependence on the initial values T * 0 and θ 0 . Under this assumption, a necessary condition for the Mpemba effect is Φ(T * 0,B , θ 0,B ) > Φ(T * 0,A , θ 0,A ), where we have considered the state A to be initially hotter than B, namely T * 0,A > T * 0,B . Next step is to ensure that the function Φ(T * , θ) exhibit a monotonic dependence on θ, so that, with an appropriate selection of the temperature ratios θ 0,A and θ 0,B , we can ensure the temperatures T * A and T * B will get closer or away from each other at the early stages. Only the first option is considered here as a simple way to attain the Mpemba effect. This will be always the case for those systems that behave monotonically and for those with non-monotonic relaxation but in which temperature evolution presents only one relative extreme value. The aim of this paper is the analysis of the latter situation. However, the determination of some criterion on whether or not temperature evolution behaves monotonically with the temperature ratio falls outside the purposes of this work.
Hence, to test the emergence of the Mpemba effect we perform the derivative of Φ with respect to θ for fixed T * . The result is which is always a positive (negative) function if γ * 2 > γ * 1 (γ * 2 < γ * 1 ). Therefore, the initial values must satisfy the following conditions Same conditions were obtained in the linear analysis developed in Ref. [13]. However, Eq. 10 does not constrain the region the initial conditions must belong to. In other words, differences between Φ(T * 0,A , θ 0,A ) and Φ(T * 0,B , θ 0,B ) must be chosen to be large enough for the occurrence of the Mpemba effect.  Figure 1. Evolution of the (reduced) temperature T * over the time t * for m 1 /m 2 = 10, σ 1 /σ 2 = 1, φ = 0.1, T * ex = 1, and x 1 = 0.5. Solid lines represent theoretical values and symbols DSMC data. From top to bottom, panel (a) corresponds to the NMME (cooling) and to the NMIME (heating) and panel (b) to the MME. To illustrate the occurrence of the non-monotonic Mpemba effect, we consider the following approach for the pair correlation functions [24] where M ℓ = i x i σ ℓ i . The non-monotonic Mpemba effect for a cooling, a heating and a mixed transition are plotted in Fig. 1. More specifically, the mixture under consideration is composed of equally distributed (x 1 = 1 2 ) hard spheres of same diameters (σ 1 = σ 2 ) and different masses (10m 1 = m 2 ) . In addition, we consider a moderate dense system (φ = 0.1). Lines are the theoretical results derived from the Enskog equation and symbols represent the DSMC data. For this parameter space, γ * 1 = 0.445 and γ * 2 = 4.451 and, consequently, the ratio T * 0,A − T * 0,B / θ 0,A − θ 0,B is chosen to be less than 0 in concordance with Eq. 10. For the Mpemba effect to occur, slopes Φ(T * , θ) are independently selected for the cooling and the heating cases (see table 1 for more details). Fig. 1(a) shows the non-monotonic Mpemba effect (NMME)-and its inverse counterpart (NMIME)-to emerge even when the initial temperature differences are around 10% (large nonmonotonic Mpemba effect). On the other hand, in Fig. 1(b) the system is initiated at two different temperatures: one cooler and the other hotter than T ex . Because of the non-linearity of the evolution equation of the temperature, a crossover is still possible. This curious effect is the so-called mixed Mpemba effect (MME) [11]. Moreover, Fig. 1 highlights an excellent agreement between the Enskog theory and the DSMC simulations. This excellent agreement ensures the accuracy of the Maxwellian approximation (Eq. 8) to model the velocity distribution function in these out-of-equilibrium systems.