Mpemba Effect in Molecular Gases under Nonlinear Drag

Physics of Fluids ARTICLE scitation.org/journal/phf Mpemba effect in molecular gases under nonlinear drag Cite as: Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 Submitted: 2 June 2020 •Accepted: 4 July 2020 • Published Online: 22 July 2020 Andrés Santos1,a) and Antonio Prados2,b) AFFILIATIONS 1Departamento de Física and Instituto de Computación Científica Avanzada (ICCAEx), Universidad de Extremadura, 06006 Badajoz, Spain 2Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, 41080 Sevilla, Spain Note: This paper is part of the Special Topic, Advances in Micro/Nano Fluid Flows: In Memory of Prof. Jason Reese. a)Author to whom correspondence should be addressed: [email protected] b)Electronic mail: [email protected] ABSTRACT We look into the Mpemba effect—the initially hotter sample cools sooner—in a molecular gas with nonlinear viscous drag. Specifically, the gas particles interact among them via elastic collisions and with a background fluid at equilibrium. Thus, within the framework of kinetic theory, our gas is described by an Enskog–Fokker–Planck equation. The analysis is carried out using the first Sonine approximation, in which the evolution of temperature is coupled to that of excess kurtosis. This coupling leads to the emergence of the Mpemba effect, which is observed at an early stage of relaxation and when the initial temperatures of the two samples are close enough. This allows for the development of a simple theory, linearizing the temperature evolution around a reference temperature, namely, the initial temperature closer to the asymptotic equilibrium value. The linear theory provides a semiquantitative description of the effect, including expressions for crossover time and maximum temperature difference. We also discuss the limitations of our linearized theory. Published under license by AIP Publishing. https://doi.org/10.1063/5.0016243 ., s I. INTRODUCTION One of the signatures of nonequilibrium systems is the presence of memory effects.1A system displays memory when its time evolution from a given initial state is not uniquely determined by the initial values of its macroscopic—or hydrodynamic—variables. In other words, the system evolution depends on how it has been previously aged; memory effects are thus intimately related to aging,2which has been typically associated with glassy behavior.3–7 Notwithstanding this, in addition to being investigated in models for glasses,8–11 it has been found in many different physical systems: granular fluids,12,13 dense granular matter,14,15 ferroelectrics,16 disordered mechanical systems,17 and frictional interfaces,18 to name just a few. The Mpemba effect19–24 is a counterintuitive memory phenomenon: given two samples of fluid, the one that is initially hotter may cool more rapidly. Therefore, the curves describing the time evolution of the temperature of the two samples cross each other at a certain time tc, and the curve for the initially hotter sample stays below the other one for longer times, t>tc. It is important to characterize the range of values of the relevant physical quantities that allow for the emergence of the Mpemba effect; in general, the difference in initial temperatures must be small enough. Although first reported in the case of water,19,25 its existence for that liquid is still controversial.26,27 As a proof of concept, the feasibility of the Mpemba effect has recently been reported in granular gases.21,28–30 Therein, collisional inelasticity couples the evolution of (granular) temperature to other quantities—such as the kurtosis or the rotational-to-translational energy ratio—monitoring the nonequilibrium nature of the velocity distribution function (VDF), even in homogeneous and isotropic states.31,32 In this work, we show that the Mpemba effect is also present in homogeneous and isotropic states of molecular gases—i.e., with elastic collisions—driven by an external drag force with a velocitydependent friction coefficient. The particles of our system are supposed to be hard spheres, for the sake of simplicity, surrounded by a background fluid in equilibrium. Gas particles collide among Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-1 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf them, with these collisions being modeled by a Boltzmann–Enskog collision term in the evolution equation for the VDF. Gas particles also interact with the background fluid. This interaction translates into two forces: (i) a macroscopic, deterministic, nonlinear drag force and (ii) a stochastic force. The intensity of the latter follows from the fluctuation–dissipation theorem, which ensures that the gas VDF tends to a Maxwellian with the temperature of the background fluid in the long-time limit. In the evolution equation, the interaction between the gas and the fluid is described by a Fokker–Planck term; therefore, the VDF obeys an Enskog–Fokker–Planck equation with nonlinear drag. The framework of our work is thus nonlinear Brownian motion,33 but the Brownian particles are no longer independent since they interact through instantaneous hard collisions. If the particles of the background fluid, with mass mbf, are much lighter than the Brownian particles, with mass m, the drag force is usually assumed to be linear in the velocity of the particles, Fdrag =−mζ0v. In fact, this is the leading behavior found when an expansion in powers of mbf/mis performed, which leads to linear Brownian motion. Nevertheless, the drag force becomes nonlinear when higher order terms in the expansion are brought to bear. Specifically, the drag force can be written as Fdrag =−mζ(v)v, and there appears a velocity-dependent drag coefficient ζ(v), with ζ(v= 0) = ζ0. Therein, the first correction in ζ0introduces a quadratic dependence on v.34–36 In some situations, nonlinearities in the drag coefficient have quite strong physical implications.36–40 The main goal of this paper is to study the Mpemba effect in the kinetic theory framework we have just described, i.e., the Enskog– Fokker–Planck equation with nonlinear drag. To meet this end, we work using the first Sonine approximation, in which the time evolution of temperature is coupled to that of excess kurtosis. This coupling, which is absent in the case of ζ(v) = ζ0, is responsible for the emergence of the Mpemba effect. The value of excess kurtosis is assumed to be small in the Sonine approximation, which entails that the initial temperatures of the samples must be close to each other and the Mpemba crossover takes place at the early stage of relaxation. This allows us to linearize the problem and derive analytical expressions for the relevant physical quantities that characterize the Mpemba effect, like the crossing time in the temperature evolution, the maximum value of the initial temperature difference, and the magnitude of the effect. The plan of the paper is as follows. We put forward the model and the kinetic description in Sec. II, where the equations for the velocity moments are also derived. Section III is devoted to the Sonine approximation: the infinite hierarchy for the velocity moments is closed by expanding the VDF in Laguerre polynomials, retaining only the first nontrivial cumulant, namely, excess kurtosis. The Mpemba effect is analyzed in Sec. IV: we develop a linearized model, investigate the crossover time, construct the phase diagram in the space of parameters, quantify the magnitude of the effect, and finally study the accuracy of the linearized theory. Finally, Sec. V presents the main conclusions of our work. II. ENSKOG–FOKKER–PLANCK EQUATION: MOMENT EQUATIONS Let us consider a d-dimensional system of elastic hard spheres of mass mand diameter σin a uniform and isotropic fluidized state. The spheres are assumed to be suspended in a background fluid in equilibrium so that their one-body VDF f(v) satisfies the Enskog– Fokker–Planck equation, ∂tf(v)−∂ ∂v⋅[ζ(v)v+ξ2(v) 2 ∂ ∂v]f(v)=J[v∣f,f]. (2.1) On the one hand, the force exerted by the background fluid on the Brownian particles has two components: a nonlinear drag force Fdrag =−mζ(v)vand a white-noise stochastic force with nonlinear variance m2ξ2(v). On the other hand, collisions between Brownian particles are accounted for by the Boltzmann–Enskog collision operator,41,42 J[v1∣f,f]=σd−1g(σ)∫dv2∫dσΘ(v12 ⋅σ)v12 ⋅σ ×[f(v′ 1)f(v′ 2)−f(v1)f(v2)]. (2.2) Therein, g(σ)=limr→σ+g(r)is the contact value of the pair correlation function g(r), Θ(⋅) is the Heaviside step function, v12 ≡v1−v2 is the relative velocity, and v′ 1=v1−(v12 ⋅σ)σ,v′ 2=v2+(v12 ⋅σ)σ(2.3) are the postcollisional velocities. Note that in the spatially uniform states under consideration, the Enskog collision operator (2.2) is simply Boltzmann’s multiplied by the factor g(σ). The velocity-dependent coefficients ζ(v) and ξ(v) are related by the condition that Eq. (2.1) admits as a stationary solution the equilibrium VDF, fs(v)=n(m 2πkBTs)d/2 e−mv2/2kBTs, (2.4) where kBis the Boltzmann constant and Tsis the equilibrium temperature of the background fluid, which acts as a thermostat. This yields the fluctuation–dissipation relation, ξ2(v)=2kBTs mζ(v). (2.5) Equation (2.1) can then be rewritten as ∂tf(v)−∂ ∂v⋅ζ(v)(v+kBTs m ∂ ∂v)f(v)=J[v∣f,f], (2.6) which describes the Brownian motion of an ensemble of particles of mass mmoving in the background fluid. These Brownian particles are not independent, with their interaction being incorporated through the collision term. In this work, as the simplest nonlinear model, we consider the quadratic dependence of the drag coefficient on the velocity derived in Refs. 34–36. Thus, we restrict ourselves to ζ(v)=ζ0(1+γmv2 kBTs), (2.7) where γ>0 is a dimensionless parameter measuring the degree of nonlinearity of the drag force. When both Brownian particles and background fluid particles are three-dimensional hard spheres, it has been shown that γ=mbf/10m.34–36 See also Appendix A. In this work, we restrict ourselves to mbf ≤2m, i.e., γ≤0.2. Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-2 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf By taking velocity moments in Eq. (2.6), the evolution equation for temperature T=m kBd⟨v2⟩=m kBnd ∫dvv2f(v), (2.8) where n=∫dvf(v) is the number density, is obtained as ˙ T ζ0=2(Ts−T)[1+γ(d+ 2)T Ts]−2γ(d+ 2)T2 Ts a2, (2.9) in which we have introduced the excess kurtosis a2=d d+ 2 ⟨v4⟩ ⟨v2⟩2−1. (2.10) In the particular case of a linear drag, γ= 0, the solution to Eq. (2.9) is simply T(t)=Ts+[T(0)−Ts]e−2ζ0t. However, in the case of nonlinear drag, γ>0, the evolution of temperature is coupled to that of excess kurtosis. Imagine that T(0) >Ts; the larger the value of a2(0), the larger the initial cooling rate is and the sooner the temperature is expected to reach the thermostat value Ts. This property can give rise to an Mpemba phenomenon, as reported in the case of granular fluids.21 Similarly, the inverse Mpemba effect, in which the cooler system heats sooner,20,21,29 may also be expected for T(0) <Ts. Since the evolution equation (2.9) involves the excess kurtosis a2(t), we need to consider its evolution equation. This in turn involves sixth-degree moments, and so on, giving rise to an infinite hierarchy of moment equations. To derive this hierarchy, let us introduce the dimensionless VDF ϕ(c) as f(v)=n vd T(t)ϕ(c),c≡v vT(t), (2.11) where vT(t)≡√2kBT(t)/mis the thermal velocity. Then, the kinetic equation (2.6) becomes ∂tϕ(c)−∂ ∂c⋅[˙ T 2Tc+ζ0(1+γ2T Ts c2) ×(c+Ts 2T ∂ ∂c)]ϕ(c)=νs√T Ts I[c∣ϕ,ϕ], (2.12) in which we have defined νs≡g(σ)nσd−1√2kBTs/m, which is basically the collision frequency at the steady state, and the dimensionless collision operator I[c1∣ϕ,ϕ]=∫dc2∫dσΘ(c12 ⋅σ)c12 ⋅σ ×[ϕ(c′ 1)ϕ(c′ 2)−ϕ(c1)ϕ(c2)]. (2.13) In addition, we have employed the property vd T n∂tf(v)=∂tϕ(c)−˙ T 2T ∂ ∂c⋅[cϕ(c)]. (2.14) Multiplying both sides of Eq. (2.12) by cℓ, integrating over c, and making use of Eq. (2.9), we obtain the hierarchy of equations for the moments Mℓ≡⟨cℓ⟩, ˙ Mℓ=ℓζ0{[γ(ℓ−2)+γ(d+ 2)T Ts(1+a2)−Ts T]Mℓ −γ2T Ts Mℓ+2 +ℓ+d−2 2 Ts TMℓ−2}−νs√T Ts μℓ. (2.15) Here, we have introduced the collisional moments μℓas μℓ≡−∫dccℓI[c∣ϕ,ϕ]. (2.16) Since, by definition, M0= 1, M2=d/2, and M4=d(d+ 2)(1 + a2)/4 [see Eqs. (2.8) and (2.10)], it is easy to check that, as it should be, ˙ M0 =˙ M2= 0 (note that μ0=μ2= 0). Next, setting ℓ= 4 in Eq. (2.15), we get ˙ a2=ζ0γ4T Ts[2Ts T(1 + a2)+(d+ 2)(1 + a2)2−(d+ 4) ×(1+3a2−a3)]−ζ04Ts Ta2−4νs d(d+ 2)√T Ts μ4, (2.17) where we have introduced the sixth-degree cumulant a3by M6=d(d+ 2)(d+ 4) 8(1+3a2−a3). (2.18) Some comments are in order. First, note that two time scales compete in Eqs. (2.15) and (2.17). The inverse of the drag coefficient for low velocities, ζ−1 0, dictates the time scale over which particles feel the action of the background fluid. Meanwhile, the characteristic time for particle–particle collisions is the inverse of the stationary collision frequency, ν−1 s. Second, Eqs. (2.9),(2.15), and (2.17) are formally exact within the Enskog–Fokker–Planck description, but they do not make a closed finite set. Not only does ˙ Mℓexplicitly involve a higher-degree moment Mℓ+2 but also the collisional moment μℓis a nonlinear functional of the full VDF ϕ(c). An approximate closure is needed to deal with a finite set of equations. III. SONINE APPROXIMATION For isotropic states, the reduced VDF ϕ(c) can be expanded in a complete set of orthogonal polynomials as ϕ(c)=e−c2 πd/2[1 + ∞ ∑ ℓ=2 aℓL(d−2 2) ℓ(c2)], (3.1) where L(α) ℓ(x)is a generalized Laguerre (or Sonine) polynomial.43 Of course, the coefficients with ℓ= 2 and ℓ=3inEq.(3.1) are the same as the cumulants a2and a3, respectively, introduced before [see Eqs. (2.10) and (2.18)]. In the first Sonine approximation, all terms beyond ℓ= 2 in Eq. (3.1) are dropped, i.e., ϕ(c)≈e−c2 πd/2{1 + a2[c4 2−d+ 2 2c2+d(d+ 2) 8]}. (3.2) Inserting Eq. (3.2) into Eq. (2.16) with ℓ= 2 and neglecting terms quadratic in a2, one obtains Γ(d/2)μ4≈√2(d−1)πd−1 2a2.44–46 Therefore, Eq. (2.17) becomes ˙ a2=8ζ0γ(1−T Ts)−{ζ0[4Ts T−8γ+ 4γ(d+ 8)T Ts] +τ−1 s 8(d−1) d(d+ 2)√T Ts⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭a2, (3.3) Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-3 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf where we have introduced the mean free time at the steady state47,48 τs=√2Γ(d/2)π1−d 2ν−1 s, (3.4) and, for consistency, the terms a2 2and a3have been neglected. Equations (2.9) and (3.3) make a closed set to investigate the existence of the Mpemba effect. First, we define dimensionless temperature and time by θ≡T Ts ,t∗≡t τs . (3.5) The latter approximately measures the accumulated number of collisions per particle up to time t. With these variables, Eqs. (2.9) and (3.3) can be rewritten as ˙ θ ζ∗ 0=2(1−θ)[1 + γ(d+ 2)θ]−2γ(d+ 2)θ2a2, (3.6a) ˙ a2 ζ∗ 0=8γ(1−θ)−[4 θ−8γ+ 4γ(d+ 8)θ+8(d−1) d(d+ 2)√θ ζ∗ 0]a2, (3.6b) where we have introduced a dimensionless low-velocity drag coefficient as ζ∗ 0≡ζ0τs. (3.7) Now, the dot over θand a2denotes a derivative with respect to t∗. Equations (3.6) are linear in the excess kurtosis a2but nonlinear in the temperature ratio θ. They constitute our starting point for the analysis of the Mpemba effect, to be carried out in Sec. IV. In the dimensionless variables we are using, there are only two relevant parameters: (i) γ, which measures the strength of the nonlinearity in the drag, and (ii) ζ∗ 0, which compares the characteristic times for collisions, τs, and for the viscous drag, ζ−1 0. Note that the regime ζ∗ 0≪1 (ζ∗ 0≫1) means that the viscous drag acts over a much longer (shorter) time scale than collisions do. IV. MPEMBA EFFECT A. Linearized model Let us imagine two initial states A and B with {θ(0),a2(0)} ={θ0 A,a0 2A}and {θ0 B,a0 2B}, respectively. The corresponding solutions to Eq. (3.6) are denoted by {θA(t∗), a2A(t∗)} and {θB(t∗), a2B(t∗)}. Without loss of generality, we assume that θ0 A>θ0 B.49 Now, we show that both the Mpemba effect and its inverse version are expected to emerge when the initially hotter sample has a larger value of the excess kurtosis. First, we analyze the case in which both initial temperatures are higher than the stationary one, θ0 A>θ0 B>1, and the system cools down to reach the steady state. The Mpemba effect is present when θA(t∗) relaxes more rapidly than θB(t∗), which calls for the existence of a crossover time t∗ csuch that θA(t∗ c)=θB(t∗ c). Since the cooling rate increases with a2, the condition a0 2A>a0 2Bseems to be necessary for the Mpemba effect to emerge. Second, we look into the case in which both temperatures are lower than the stationary value, 1 >θ0 A>θ0 B, and the system heats up. The inverse Mpemba effect appears if θA(t∗) relaxes more slowly than θB(t∗), which again needs a0 2A>a0 2B. In general, the nonlinear dependence on θof the set of Eq. (3.6) impedes a fully analytical treatment. However, the excess kurtosis is supposed to be small in the Sonine approximation, which has allowed us to neglect nonlinear terms in a2. Since a2is the quantity controlling the appearance of the Mpemba effect, its smallness implies that both initial temperatures, θ0 Aand θ0 B, cannot be very far from each other for the Mpemba effect to emerge. In addition, the crossing of the curves θA(t∗) and θB(t∗) must take place at the early stage of evolution.50 Following the discussion above, we write θ(t∗) = θr+Ψ(t∗), where θr≈θ0 A≈θ0 Bis a certain reference temperature, and linearize Eq. (3.6) with respect to Ψ(t∗) and a2(t∗). The detailed solution of this linearization procedure is carried out in Appendix B; here, we present the results relevant for the analysis of the Mpemba effect. The time evolution is controlled by the matrix with elements Λ11 =2ζ∗ 0[1 + γ(d+ 2)(2θr−1)], (4.1a) Λ12 =2ζ∗ 0γ(d+ 2)θ2 r, (4.1b) Λ21 =8ζ∗ 0γ, (4.1c) Λ22 =ζ∗ 0[4 θr−8γ+ 4γ(d+ 8)θr]+8(d−1) d(d+ 2)√θr, (4.1d) and eigenvalues λ±=Λ11 +Λ22 ±δλ 2,δλ≡√(Λ11 −Λ22)2+ 4Λ12Λ21. (4.2) Let us consider the differences Δθ(t∗)≡θA(t∗)−θB(t∗) and Δa2(t∗)≡a2A(t∗)−a2B(t∗) between the time evolutions corresponding to the two different initial states A and B. Within the linearized theory, these differences are given by (see Appendix B) Δθ(t∗)=(λ+−Λ11 δλ Δθ0−Λ12 δλ Δa0 2)e−λ−t∗ −(λ−−Λ11 δλ Δθ0−Λ12 δλ Δa0 2)e−λ+t∗ , (4.3a) Δa2(t∗)=(λ+−Λ22 δλ Δa0 2−Λ21 δλ Δθ0)e−λ−t∗ −(λ−−Λ22 δλ Δa0 2−Λ21 δλ Δθ0)e−λ+t∗ . (4.3b) Note that both Δθand Δa2vanish in the long-time limit. B. Mpemba crossover The accuracy of the linearized theory developed above for describing the Mpemba effect—and the inverse Mpemba effect—is illustrated in Fig. 1. The linear theory remains valid even when the system is initially far from the steady state; the analytical expressions of the linearized theory, given by Eq. (B4), predict the crossover of the curves correctly but start to deviate from the “exact” numerical integration as time grows. For all these plots, an optimal choice for θris the initial temperature of the sample that is closer to the Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-4 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf FIG. 1. First stage in the evolution of θA(t∗) and θB(t∗) [panels (a), (c), (e), and (g)] and a2A(t∗) and a2B(t∗) [panels (b), (d), (f), and (h)] for d= 3, ζ∗ 0=1, and γ= 0.1. The initial states are a0 2A=0.5, a0 2B=−0.35, and [(a) and (b)] θ0 A=10, θ0 B=9, [(c) and (d)] θ0 A=2,θ0 B=1.8, [(e) and (f)] θ0 A=1.1,θ0 B=1.05, and [(g) and (h)] θ0 A=0.9,θ0 B=0.85. Circles correspond to the numerical solutions of Eq. (3.6), whereas solid lines correspond to the linearized model, given by Eq. (B4). The Mpemba effect is neatly observed in panels (a), (c), and (e). The linearized theory with θr=θ0 Bgives a correct account thereof—although it deviates from the numerical solution of Eq. (3.6) as time increases in panels (a) and (c), for which their initial temperatures are not close to the steady value. The inverse Mpemba effect is depicted in panel (g), and the linear theory also describes it correctly, but now we have chosen θr=θ0 A. steady state. Figure 1 also shows that the excess kurtosis relaxes to equilibrium more rapidly than temperature. Let us now restrict ourselves to a situation in which the Mpemba effect is present. Thus, there exists a crossover time such that Δθ(t∗ c)=0. According to Eq. (4.3a), it is given by t∗ c=1 δλ ln Λ12 −(λ−−Λ11)R0 Λ12 −(λ+−Λ11)R0,R0≡Δθ0 Δa0 2 . (4.4) Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-5 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf The crossover time t∗ cdepends on the initial preparation only through the reference temperature θr≈θ0 A≈θ0 Band the ratio R0 in this simplified description, for given values of ζ∗ 0and γ. Note that we have chosen Δθ0>0 and, for the Mpemba effect to exist, we need that Δa0 2>0, i.e., we have that R0>0. Figure 2 displays t∗ cas a function of R0for some illustrative cases. Different panels correspond to different values of the reference temperature. In all of them, the crossover time t∗ cvanishes in the limit as R0→0 and grows with R0.Figure 2 also includes the values of the crossover time obtained from the numerical solution of Eq. (3.6) for a0 2A=0.5 and a0 2B=−0.35 with θ0 B=θr in panels (a)–(c) and θ0 A=θrin panel (d). It is observed that the agreement with Eq. (4.4) improves as γincreases and R0decreases. In addition, Eq. (4.4) underestimates the crossover time for the direct Mpemba effect with initial temperatures far from that of the thermostat, while it tends to overestimate t∗ cfor the inverse Mpemba effect or when the initial temperatures are close to the thermostat one. Equation (4.4) shows that t∗ cdiverges in the limit as R0→R0 th, where R0 th isathreshold value for the ratio, given by R0 th =Λ12 λ+−Λ11 . (4.5) Thus, the Mpemba effect disappears if R0≥R0 th: in this region, t∗ c, as defined by Eq. (4.4), ceases to be a real number. In fact, if we define β≡λ−−Λ11 λ+−Λ11 =1−δλ λ+−Λ11 , (4.6) we can rewrite t∗ cas t∗ c=1 δλ ln 1−βR0/R0 th 1−R0/R0 th . (4.7) The emergence of the Mpemba effect is basically controlled by the strength of the drag nonlinearity γ. As expected on a physical basis, the Mpemba crossover takes place earlier as γincreases. Throughout Fig. 2, the curves for γ= 0.1 lie above those for γ= 0.2. Thus, the smaller γis, the smaller the threshold value R0 th we find. Recall that the drag becomes linear in the limit as γ→0+, for which the temperature obeys a closed first-order differential equation, independent of the value of the excess kurtosis, and the Mpemba effect is no longer present. Note also that terms beyond the quadratic one in the expansion of ζ(v) in powers of vmight be necessary as γ(or, equivalently, mbf/m) increases. C. Phase diagram A phase diagram in the (γ,R0) plane can be constructed, as illustrated in Fig. 3. The line R0=R0 th separates the regions in which the Mpemba effect is present (R0<R0 th) and absent (R0>R0 th), where R0 th is defined in Eq. (4.5). The range 0 <R0<R0 th for which the Mpemba effect emerges increases with γ,ζ∗ 0, and θr. It must be remarked that θr<1 corresponds to the inverse Mpemba effect, in which the system relaxes to equilibrium from below the steady temperature Ts. The threshold values R0 th obtained from the numerical solution of Eq. (3.6) for a0 2A=0.5 and a0 2B=−0.35 are also shown in Fig. 3. FIG. 2. Crossover time t∗ cas a function of R0≡Δθ0/Δa0 2[see Eq. (4.4)]. All the panels correspond to d= 3 and ζ∗ 0=1, but to different values of the reference temperature: θr= 10, 2, 1.05, and 0.5, from (a) to (d). In each panel, two values of the nonlinearity parameter are considered: γ= 0.1 (upper curves) and 0.2 (lower curves). The crossover time decreases as the nonlinearity coefficient γincreases, as expected. The symbols are the values obtained from the numerical solution of Eq. (3.6) with a0 2A=0.5 and a0 2B=−0.35. Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-6 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf FIG. 3. Phase diagram in the plane R0≡Δθ0/Δa0 2vs γfor d= 3. The values of ζ∗ 0are (a) ζ∗ 0=0.01 (solid curves and closed circles) and 0.1 (dashed curves and open triangles) and (b) ζ∗ 0=1(solid curves and closed circles) and 2 (dashed curves and open triangles). From the bottom to the top in each panel, θr= 0.5, 1.05, 2, and 10. In each case, the Mpemba effect is present (absent) below (above) the corresponding curve. Symbols are obtained by numerically solving Eq. (3.6) with a0 2A=0.5 and a0 2B=−0.35, whereas lines correspond to the analytical prediction given by Eq. (4.5). The agreement with the simplified model is quite good, especially in panel (a), for which ζ∗ 0<1. In panel (b), for which ζ∗ 0≥1, the linearized theory still gives a semi-quantitative picture and, notably, successfully captures the weak influence of both γand ζ∗ 0on R0 th if ζ∗ 0≥1. In any case, the linear model overestimates (underestimates) R0 th for θr= 10 and 2 (θr= 1.05 and 0.5), as anticipated from Fig. 2. Interestingly, the maximum ratio relative to the reference temperature, θ−1 rR0 th, keeps increasing with increasing θr. At fixed θr, the upper bound of R0 th corresponds to the limit γ→∞, which is independent of ζ∗ 0, namely, lim γ→∞R0 th =2(d+ 2)θ2 r/(d−2 + 12θr) 1 + √1 + 16(d+ 2)θ2 r/(d−2 + 12θr)2. (4.8) Notwithstanding, in our modeling, we are only retaining the first correction, quadratic in the velocities, in the drag coefficient ζ(v). Therefore, from a physical point of view, γis expected not to be very large; otherwise, higher order terms in the velocity should be incorporated into the drag coefficient. D. Magnitude of the Mpemba effect When the Mpemba effect is present, the temperature difference Δθvanishes at the crossover time t∗ c. Since Δθalso vanishes in the long time limit as t∗→∞, there must exist a certain time t∗ m>t∗ c where |Δθ(t∗)| reaches a local maximum. Therefore, one has that ∣Δθ(t∗)∣≤∣Δθ(t∗ m)∣for any time t∗>t∗ c. The above discussion can be used to introduce a quantitative measure of the magnitude of the Mpemba effect. Let us define28 Mp ≡∣Δθ(t∗ m)∣ (4.9) as a quantitative measure of its magnitude. From Eq. (4.3), one finds t∗ m=t∗ c+1 δλ ln λ+ λ− , (4.10a) Mp Δa0 2=Λ12(1−R0/R0 th λ+)λ+ δλ(λ− 1−βR0/R0 th )λ− δλ. (4.10b) Thus, the Mpemba magnitude Mp depends on the initial differences Δθ0and Δa0 2by a simple scaling law: the ratio Mp/Δa0 2is a function of the ratio R0≡Δθ0/Δa0 2.Figure 4 shows that the larger the ratio R0is, the smaller Mp/Δa0 2becomes. Of course, Mp vanishes as R0approaches its threshold value R0 th, as readily seen in Eq. (4.10b). Comparison with the values of Mp obtained from the numerical solution of Eq. (3.6) shows that the simple linearized model is qualitatively correct in capturing the dependence of the order of magnitude of Mp on the parameters of the problem. In agreement with what was observed in Figs. 2 and 3, the linearized model tends to overestimate (underestimate) Mp for the direct (inverse) Mpemba effect. Figure 4(c) shows that the prediction of Mp is especially accurate if the initial temperatures are close to the equilibrium one; in that case, Mp is slightly overestimated for small R0, while it is slightly underestimated as R0approaches its threshold value R0 th. E. Reliability of the linear theory The linear theory we have presented does not apply for all times, unless the reference temperature θr=θs= 1 and both initial temperatures θ0 Aand θ0 Bare close to the steady state. Nevertheless, as already stated before, our linear theory is not the standard linearization around the steady state but an approximate scheme to obtain a good approximation to the actual time evolution of the system at the early stage of its evolution, where the Mpemba effect is expected to come about. This means that our linear approach does have some limitations, as observed in Figs. 2–4. According to them, the linearized model becomes more accurate as |θr−1| and ζ∗ 0 decrease. While a complete account of the range of validity of the linear approximation is outside the scope of our paper, quite simple arguments can be presented in the limit of weak nonlinearity, γ→0+. The behavior of the threshold R0 th—below which the Mpemba effect is found—depends on that of its numerator, Λ12, and its denominator, λ+−Λ11. On the one hand, Eq. (4.1b) tells us that the former is linear in γ,Λ12 =O(γ). On the other hand, the behavior of λ+−Λ11 for small γdepends on the sign of the function φ(θr,ζ∗ 0)≡4(d−1) d(d+ 2)ζ∗ 0 θ3/2 r−θr+ 2. (4.11) Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-7 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf FIG. 4. Magnitude of the Mpemba effect, Mp/Δa0 2, as a function of R0≡ Δθ0/Δa0 2. The parameter values are the same as that given in Fig. 2: (i) in all the panels, we have that d= 3 and ζ∗ 0=1and (ii) from panel (a) to (d), different values of the reference temperature are considered, namely, θr= 10, 2, 1.05, and 0.5, with two values of the nonlinearity parameter being considered in each panel, γ= 0.1 (lower curves) and 0.2 (upper curves). Again, symbols are obtained by numerically solving Eq. (3.6) with a0 2A=0.5 and a0 2B=−0.35, whereas lines correspond to the analytical prediction given by Eq. (4.10b). Specifically, it can be readily shown that λ+−Λ11 ≈⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 2ζ∗ 0φ(θr,ζ∗ 0) θr ,φ(θr,ζ∗ 0)>0 −8ζ∗ 0(d+ 2)θ3 r φ(θr,ζ∗ 0)γ2,φ(θr,ζ∗ 0)<0, (4.12) and thus, R0 th ∼⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩(d+ 2)θ3 r γ φ(θr,ζ∗ 0),φ(θr,ζ∗ 0)>0 −1 4θr φ(θr,ζ∗ 0) γ,φ(θr,ζ∗ 0)<0. (4.13) In the limit as γ→0+, the drag becomes linear, the temperature obeys a closed equation, and no Mpemba effect can be present in the system. This is consistent with the behavior found for the threshold R0 th when φ(θr,ζ∗ 0)>0; therein, R0 th →0. However, R0 th →∞when φ(θr,ζ∗ 0)<0; this is an unphysical result that makes us conclude that the simplified linear model (4.3) ceases to be reliable if φ(θr,ζ∗ 0)<0 and γ≪1. The locus φ(θr,ζ∗ 0)=0 is plotted in Fig. 5. The above discussion should not be employed to disregard the linear model directly when φ(θr,ζ∗ 0)<0; the linearization can be useful unless γis very small. We can estimate the value of γfor which the linear theory is no longer accurate by asking the estimate for R0 th in (4.13) to be large. This leads to the condition γ≪γℓ≡∣φ(θr,ζ∗ 0)∣ 4θr . (4.14) We illustrate the above result in Fig. 6 for the three-dimensional case, specifically for ζ∗ 0=5 and θr= 9, in which case we have φ(θr,ζ∗ 0)≃−4.12 and γℓ≃0.114. While the (weak) Mpemba effect predicted by the linearized model with γ= 0.001 is actually absent, FIG. 5. Locus φ(θr,ζ∗ 0)=0in the plane ζ∗ 0vs θrfor d= 3. The function φ(θr,ζ∗ 0), given by Eq. (4.11), is negative above (positive below) the curve. The circle represents the minimum θr= 6, ζ∗ 0=4 5√6≃1.96, whereas the dashed line at θr= 2 is the vertical asymptote of the left branch of the locus. Thus, ζ∗ 0 must be large enough, specifically larger than its value at the minimum, to allow for negative values of φ. Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-8 Published under license by AIP Publishing 23 May 2025 14:04:34 Physics of Fluids ARTICLE scitation.org/journal/phf FIG. 6. Evolution of Δθ(t∗) = θA(t∗)−θB(t∗) for d= 3, ζ∗ 0=5. The two panels correspond to small values of γ, namely, (a) γ= 0.001 and (b) γ= 0.1. The initial states are {θ0 A,a0 2A}={10, 0.5}and {θ0 B,a0 2B}={9, −0.35}. Circles correspond to the numerical solutions of Eq. (3.6), whereas solid lines correspond to the linearized model (4.3) with θr=θ0 B. The value of γin panel (a) verifies condition (4.14) that controls the failure of the linearized theory. the model succeeds in locating the crossover time if γ= 0.1, which is quite close to γℓ. Note that γ= 0.1 corresponds to the case in which the mass of the Brownian particles and that of the surrounding fluid are identical, as shown in Appendix A. V. CONCLUSIONS We have neatly observed the Mpemba effect in a molecular gas with nonlinear drag. For the Mpemba effect—and for the inverse Mpemba effect, in which the initially cooler system heats sooner— to emerge, the initially hotter sample must have a sufficiently larger value of excess kurtosis a0 2; the larger a0 2is, the larger the cooling rate becomes. This behavior is completely analogous to that found in a granular gas of smooth hard spheres.21 The above analysis entails that the Mpemba effect is absent if both samples, A and B, are initially at equilibrium. In that case, a0 2A=a0 2B=0, and the parameter R0defined in Eq. (4.4) diverges to infinity. For the Mpemba effect to emerge, we need to prepare the samples in nonequilibrium states before coupling them to the common reservoir at temperature Ts. This can be achieved, for instance, by temporarily coupling the samples to their respective reservoirs at temperatures different from Ts. Analytical predictions have been obtained, in a wide range of values of the system parameters, within a linearized model. The linearization is carried around a reference temperature—specifically the initial temperature of the sample that is closer to the equilibrium value, not around the steady temperature. Therefore, our analytical framework is not limited to near-equilibrium situations. Within this scheme, we have found semi-quantitatively accurate expressions for (i) the crossover time, (ii) the maximum value of the initial temperature difference, and (iii) the magnitude of the Mpemba effect. In addition, we have looked into the limitations of the linearized model, especially for small values of the parameter γcharacterizing the nonlinearity. This work also opens avenues for further research. It is interesting to consider in more detail some specific limits of the present model, which are physically relevant: (i) small nonlinearity γ≪1, which appears naturally as the first correction to the usual linear drag and (ii) time scale separation between viscous drag and collisions, i.e., either ζ∗ 0≪1 or ζ∗ 0≫1. In both cases, a systematic—mainly perturbative—analytical approach seems to be feasible. In addition, it is important to deepen our understanding of the aging phenomena in this molecular gas. Specifically, looking into the Kovacs effect,51,52 which has attracted a lot of attention lately,17,18,30,53–63 is compelling. Finally, we plan to check the theoretical predictions of this work against computer simulations obtained from the Langevin equation with an interpretation of the multiplicative noise33,64–66 consistent with the Enskog–Fokker–Planck equation (2.6). ACKNOWLEDGMENTS Professor Jason Reese was a brilliant scientist in the field of kinetic theory, with very important contributions at both the theoretical and the computational level. With this paper, we respectfully pay tribute to his memory. A.S. and A.P. acknowledge financial support from the Spanish Agencia Estatal de Investigación through Grant Nos. FIS2016-76359-P and PGC2018-093998-B-I00, respectively. A.S. is also grateful to the Junta de Extremadura (Spain) for Grant No. GR18079. All the grants are partially financed by the European Regional Development Fund. APPENDIX A: HARD INTERACTION BETWEEN BROWNIAN PARTICLES AND BACKGROUND FLUID In Refs. 34–36, the authors consider the emergence of nonlinear Brownian motion when an ensemble of “heavy” Brownian particles (mass mand number density n) moves in a bath modeled as a dilute gas of “light” particles (mass mbf and number density nbf), which is at equilibrium at temperature Ts. A velocity-dependent drag coefficient ζ(v) is obtained as an expansion in powers of the mass ratio mbf/m, which is formally assumed to be small. The coefficients of the expansion of ζ(v) are given in terms of integrals that involve the differential cross section for the interaction between the Brownian particles and the bath particles. Explicit expressions for ζ(v) can be derived when simple potentials are employed for this interaction. For example, all particles are considered to be three-dimensional hard spheres in Ref. 36, both the Phys. Fluids 32, 072010 (2020); doi: 10.1063/5.0016243 32, 072010-9 Published under license by AIP Publishing 23 May 2025 14:04:34