The metastable Mpemba effect corresponds to a non-monotonic temperature dependence of extractable work

The Mpemba effect refers to systems whose thermal relaxation time is a non-monotonic function of the initial temperature. Thus, a system that is initially hot cools to a bath temperature more quickly than the same system, initially warm. In the special case where the system dynamics can be described by a double-well potential with metastable and stable states, dynamics occurs in two stages: a fast relaxation to local equilibrium followed by a slow equilibration of populations in each coarse-grained state. We have recently observed the Mpemba effect experimentally in such a setting, for a colloidal particle immersed in water. Here, we show that this metastable Mpemba effect arises from a non-monotonic temperature dependence of the maximum amount of work that can be extracted from the local-equilibrium state at the end of Stage 1.


I. INTRODUCTION
A generic consequence of the second law of thermodynamics is that a system, once perturbed, will tend to relax back to thermal equilibrium. Such relaxation is typically exponential. To understand why, consider energy relaxation and recall that the heat equation, for the temperature field T at position r and time t and thermal diffusivity κ has solutions that can be written in the form 1 Here, T ∞ (r) is the static temperature-field solution of Equation (1); it must account for boundary conditions. For t → ∞, an arbitrary initial condition T (r, 0) will relax to this state. In Equation (2), the v m (r) are spatial eigenfunctions, with corresponding eigenvalues λ m and coefficients a m , which represent the projections of the field [T (r, 0) − T ∞ (r)] onto the corresponding eigenfunction. For long but finite times, all but the slowest eigenmode will have decayed, and the temperature is, approximately, which, indeed, shows a simple exponential decay to T ∞ (r) for a probe at a fixed position r. * raphael.chetrite@unice.fr 1 We begin the eigenfunction expansion at m = 2 to be consistent with the analogous expansion of the Fokker-Planck solution in Lu and Raz [1].
Although exponential decays are typical, anomalous, non-exponential relaxation is also encountered. Large objects, for example, may have an asymptotic time scale λ −1 2 that exceeds experimental times, so that it is not possible to wait "long enough." Similarly, glassy systems and other complex materials may have a spectrum of exponents for mechanical and dielectric relaxation that have not only very long time scales but also many closely spaced values that are not resolved as a sequence of exponentials. Rather, they can collectively combine to approximate a power-law or even logarithmic time decay, with specific details that depend on the history of preparation [2].
Another class of anomalous systems shows unexpectedly fast relaxation in certain circumstances. The bestknown of these is the observation that, occasionally, a sample of hot water may cool and begin to freeze more quickly than a sample of cool or warm water prepared under identical conditions. Based on the scenario of exponential relaxation sketched above, one's naive intuition is that a hotter system will have to "pass through" all intermediate temperatures and thus take longer to equilibrate. More succinctly, the observation is that, in some systems, the equilibration time is a non-monotonic function of the initial temperature: the time for a system initially in equilibrium at a given temperature takes to cool and reach equilibrium with the bath temperature does not always increase with initial temperature.
While observations of this phenomenon date back two millennia to the ancient Greeks [3,4], its modern study began with observations by Mpemba in the 1960s [5]. The effect has since been observed in systems such as manganites [6], clathrates [7], polymers [8,9] and predicted in simulations of other systems, including carbon nanotube resonators [10], granular fluids [11], and spin glasses [12]. In all these Mpemba effects, the relaxation time shows a surprisingly complicated dependence on the deviation of initial temperature from equilibrium: increasing and then decreasing, and in some cases increas-ing again with increasing deviation. The relaxation time thus does not increase monotonically with the deviation from equilibrium, as one might naively expect.
One challenge in studying Mpemba effects is that the systems where they have been observed or predicted have been rather complicated, with many possible explanations for the effect. The explanations tend to be complicated and specific to a particular system. Even water is not as simple as it might seem: proposed mechanisms include evaporation [13][14][15], convection [16], supercooling [17], dissolved gases [18], and effects arising from hydrogen bonds [19].
In an effort to understand the Mpemba effect more generically, Lu and Raz recently proposed an explanation that is linked to the structure of eigenfunction expansions such as that in Equation (2) [1]. Their work was formulated for mesoscopic systems that are in the classical regime yet are small enough that thermal fluctuations make an important contribution to their dynamics. Such systems may be described by master equations and Fokker-Planck equations, for finite and continuous state spaces, respectively [20][21][22][23][24]. For the latter, the Fokker-Planck equation describes the evolution of the probability density function p(x, t) for a system described by a state vector x(t). 2 Its structure is similar to that of Equation (1): its linearity implies that solutions are also described by an infinite-series, eigenfunction expansion similar to that in Equation (2). The essence of Lu and Raz's explanation is that the projection of the initial state p(x, 0) -a Gibbs-Boltzmann distribution corresponding to an initial temperature T -onto the slowest eigenfunction, a 2 can be non-monotonic in T , or, equivalently, in β −1 ≡ k B T , where k B ≡ 1 (in our units) is Boltzmann's constant. Such a consequence implies a Mpemba effect because the long-time limit for the probability density function has the same form as Equation (3): with g β b (x) the Gibbs-Boltzmann distribution for the system at a temperature T b corresponding to the surrounding thermal bath with which the system is in contact and can exchange energy. The coefficient a 2 is a function of both the initial temperature and bath temperature: where the initial state p(x, 0) is assumed to be in equilibrium at a higher temperature β −1 and where u 2,β b (x) is the left eigenfunction of the Fokker-Planck operator, which is the dual-basis element corresponding to the right eigenfunction v 2,β b (x) of the same Fokker-Planck operator. Both u 2 and v 2 are calculated for the Markovian 2 In a many-body system, the dimension of x can be very large.
Langevin dynamics associated with white noise whose covariance is set by the bath temperature, β −1 b . We need to distinguish between left and right eigenfunctions because the operator generating Fokker-Planck dynamics is not self-adjoint, in contrast to the operator generating the heat-diffusion dynamics discussed in Equation (1). The Mpemba effect then translates to the non-monotonicity of a 2 as a function of the initial temperature β −1 : If a high-temperature initial condition has a smaller coefficient a 2 , then, in the long-time limit, the system will be closer to equilibrium than a cool-temperature initial condition with larger a 2 . This non-monotonicity in a 2 is easier to establish than the non-monotonicity of equilibration times that defines the Mpemba effect. The latter requires either an experiment or, at the very least, repeated numerical solution of the full Fokker-Planck equation.
Inspired by the scenario proposed by Lu and Raz [1], we have explored the Mpemba effect in a simple, mesoscopic setting that -unlike previous work -lends itself to quantitative experiments that straightforwardly connect with theory [25]. In particular, we explored the motion of a single micron-scale colloidal particle immersed in water and moving in a tilted double-well potential. The one-dimensional (1D) state space consists of the position x(t) of the particle. By choosing carefully the tilt of the potential, along with the energy-barrier height and the offset (asymmetry) of the double-well potential within a box that confines the particle motion at high temperatures, we could demonstrate convincingly the existence of the Mpemba effect and measure the non-monotonic temperature dependence of the a 2 coefficient. We even found conditions where a 2 = 0. At such a point, the slowest relaxation dynamics is ∼ e −λ3t , implying an exponential speed-up over the generic relaxation dynamics, ∼ e −λ2t . This strong Mpemba effect had been predicted by Klich et al. [26].
Although our recent experimental work gives strong support to the basic scenario proposed by Lu and Raz, it does not offer good physical insight into the conditions needed to produce or observe the Mpemba effect. What physical picture corresponds to the anomalous temperature dependence of the a 2 coefficient? In this Brief Research Report, we offer a more physical interpretation of the Mpemba effect explored in our previous work.

II. THERMALIZATION IN A DOUBLE-WELL POTENTIAL WITH METASTABILITY
A common feature of experiments showing Mpemba effects is that they involve a temperature quench: the system is cooled very rapidly. We model this situation by making the high-temperature initial state an initial condition for dynamics that take place entirely in contact with a bath of fixed temperature. In effect, the quench is infinitely fast. The thermalization dynamics are then given by the Langevin equatioṅ with γ a friction coefficient and η(t) Gaussian white noise modeling thermal fluctuations from the bath, with enforces the fluctuation-dissipation relation [21,23]. The potential U (x) is a double-well potential with barrier height E 0 ≫ β b −1 and two coarsegrained states, denoted L and R in Figure 1A. The range of particle motions is also constrained to a finite range; the potential is implicitly infinite at the extremities. By tilting the potential, one state has a higher energy than the other (difference is ∆E) and becomes a toy model for the water-ice phase transition. However, the energy barrier E 0 , while high enough that the two states are well defined, is also low enough that many transitions over the barrier are observed during a typical experiment. Figure 1 illustrates the case studied by Kumar and Bechhoefer [25], with (a) showing the potential and (b) the dynamics of a quench from a high temperature. With a moderately high barrier, both wells have significant probability for the equilibrium state g β b (x) ( Figure 1B, right). For U (x), the barrier E 0 = 2.0, the energy difference between states is ∆E = 1.3, and the hot temperature β −1 h = 1000; all quantities are multiplied by β b and are, hence, dimensionless. The separation between wells was 80 nm. At a temperature corresponding to β −1 , the equilibrium free energy of the system is and the corresponding equilibrium Gibbs density is The metastability of U means that the system evolves on two very different time scales: Stage 1 is a fast relaxation to local equilibration. The initial, high-temperature Gibbs density rapidly evolves to a state that is at local equilibrium with respect to the bath temperature. A local equilibrium is a density that is similar locally to g β b , but with altered fractions of systems in the left or right wells. Using Bayes' theorem, we can write such a local-equilibrium state as ρ leq β,β b (x) = P ( be in the left well at β) P ( x| be in the left well at β b ) + P ( be in the right well at β) P ( x| be in the right well at β b ) .
More precisely, the local β, β b equilibrium is the density x > 0 (right well) , with 0 ≤ a L ≤ 1. Choosing a L + a R = 1 ensures normalization of the probability density.
In a fast quench, we assume that the fraction of initial systems at equilibrium at the higher temperature β −1 is unchanged when local equilibrium is established. Essentially, we ignore the diffusion of trajectories that start on one side of the barrier and end up on the other at the end of the transient. In this approximation, the fraction that ends up in each well corresponds to that of the initial state, g β . Thus, and As shown in Figure 1B, center, the local-equilibrium distribution ρ leq β,β b (x) is discontinuous at x = 0; higher barriers will reduce the discontinuity, of order e −β b E0 ≪ 1.
Stage 2 is a final relaxation to global equilibrium on a slow time scale: the overall populations in each well (coarse-grained state) change, and the density converge to the Gibbs density g β b . Local equilibrium is maintained during the evolution, which is illustrated schematically in Fig 1b. In this metastable regime, the equilibration time was analyzed by Kramers long ago [21,[27][28][29]. It also corresponds to the limit of Equation (4); as a result, the final relaxation is exponential, with decay rate λ 2 .

III. METASTABLE MPEMBA EFFECT
Given this scenario of thermal relaxation as a two-stage process, we can readily understand how the Mpemba effect can occur. The idea is to follow the dynamics in the function space of all admissible probability density functions p(x, t). If we expand the solution in eigenfunctions analogously to Equation (2), we see that the infinitedimensional function space is spanned by the eigenfunctions. To visualize the motion, we project it onto the 2D subspace spanned by the eigenfunctions v 2 (x) and v 3 (x). The system state is then characterized as a parametric plot of the amplitudes a 2 (t) and a 3 (t). Animations from a 3D projection spanning a 2 -a 3 -a 4 are available in the supplementary material. A similar geometric plot was used to explore quenching in an anti-ferromagnetic Ising spin system in [26]. Figure 2 shows the geometry of trajectories. They are organized about two static, 1D curves, labeled G and G leq . The red curve (G) represents the set of all equilibrium Gibbs-Boltzmann densities, g β , for 0 ≤ β < ∞. It is sometimes known as the quasi-static locus. The green curve (G leq ) represents the set of all local-equilibrium densities of the form of Equation (9) Probability-density dynamics in the plane defined by the a2 and a3 coefficients. The red curve G denotes the set of equilibrium densities, the green curve G leq the set of localequilibrium densities. Arrows indicate the slow relaxation along G leq to global equilibrium, at the intersection with G (denoted by the large hollow marker with a dot at its center at T = 1). Gray lines denote the rapid relaxation from an initial condition (temperature relative to the bath indicated by a marker along G). The time progression of p(x, t), projected onto the a2-a3 plane, is from dark to light. Curves are calculated from the double-well potential described in Kumar and Bechhoefer, with domain asymmetry α = 3 (see [25] for definitions).
is not shown in the figure.) The two curves intersect at a 2 = a 3 = 0, which describes the global equilibrium g β b with respect to the bath (large hollow marker with dot). The apparent crossing near a 2 ≈ 0.4 is spurious, as the 3D projections in the supplement show. The dynamical trajectories are represented by the variously shaded gray curves. At time t = 0, the systems are in equilibrium along the red curve at a variety of temperatures {1, 1.2, 1.5, 3, 50, 100, 1000} × T b , which are indicated by black markers. The curves then move rapidly towards the green curve (local equilibrium). The time course is suggested by the dark-to-light gradient. Once they reach the vicinity of G leq , they closely follow this green curve back to the global-equilibrium state.
Within this representation, we note the "arrival point" of each trajectory when it "hits" G leq . For small temperatures (1, 1.2, 1.5, 3), the distance between this arrival point and the global-equilibrium state increases monotonically with β. For larger temperatures (50, 100, 1000), however, the distance decreases until, at T = 1000T b , it nearly vanishes (denoting the strong Mpemba effect). Along G leq , the system is in the limit described by Equation (4) and relaxes exponentially to global equilibrium. Relaxation along G leq therefore must be monotonic with the distance away from global equilibrium. Trajectories that arrive along this curve that are farther from global equilibrium will take longer to relax. In the Appendix, we show that this notion of "distance" along G leq can be expressed as the Kullback-Leibler divergence D KL between the local equilibrium density ρ leq β,β b given in Equation (9) and the global equilibrium density g β b . In particular, is a monotonic function of a L (defined in Equation 10), which is the natural parameter for the manifold G leq . Now we can understand how the (metastable) Mpemba effect can arise. In the example shown in Figure 2, the distance along G leq initially increases with T and so does the total equilibration time. But then this distance decreases for higher temperatures, leading to the Mpemba effect. We note that in our approximation, the time to traverse the initial stage is much shorter than the time to relax along the green curve, so that variations in the length of the initial trajectory are irrelevant.
If the bath temperature were changed at a finite rate (rather than a hot system being quenched directly into the bath), then the dynamics would be different. For example, if the system is very slowly cooled from the initial temperature to final bath temperature, the trajectory would follow the quasi-static locus (red curve G) and no Mpemba effect would be possible. Having shown that no Mpemba effect is possible with an infinitely slow quench and that the effect can be observed in the limit of an infinitely rapid quench, we can conclude that the Mpemba effect requires a sufficiently fast temperature quench.

IV. METASTABLE MPEMBA EFFECT IN TERMS OF EXTRACTABLE WORK
Our main goal is to express the criterion for the Mpemba effect in more physical terms. For the metastable setting described above, we will find such a criterion in terms of a thermodynamic work. We recall that the second law of thermodynamics for a system in contact with a single thermal bath of temperature β −1 b can be expressed in terms of work and free energy rather than entropy: where W is the work received by the system and △F neq denotes the difference in nonequilibrium free energies (final − initial values). See, for example, Gavrilov et al. [30], Equation 5 and associated references. We recall also that the nonequilibrium free energy generalizes the familiar notion of free energy to systems out of equilibrium. Thus, in analogy to Equation 7, we define where the average energy E(ρ) and Gibbs-Shannon en-tropy S(ρ) are given by These expressions reduce to their usual definitions for ρ = g β b but can be evaluated, as well, over nonequilibrium densities.
In the formulation of the second law of Equation (11), the initial and final states are arbitrary. In our case, the initial state is the (approximate) local equilibrium reached at the end of Stage 1. In the final state, the system is in equilibrium with the bath.
Physically −△F neq represents the maximum amount of work that may be extracted from the nonequilibrium isothermal protocol [31]. We will refer to this quantity as the extractable work.
In the Appendix, we show that the difference in nonequilibrium free energies △F neq may be expressed as a Kullback-Leibler divergence. Explicitly, In our set-up, the extractable work between the "intermediate" time (end of Stage 1) where F neq,β b = F neq,β b ρ leq β,β b , and the final time of the slow evolution (where F neq,β b = F eq,β b ), is given by Equation (15): In Sec. III and Figure 2, we saw that D KL (ρ leq β,β b , g β b ) can be non-monotonic as a function of β. We thus conclude that there can be a non-monotonic dependence on β of the function This is our main result: If the metastable Mpemba effect occurs, then the extractable work from the localequilibrium state at the end of Stage 1 is non-monotonic in the initial temperature β −1 . Figure 3 shows an example, again calculated for the potential considered by Kumar and Bechhoefer [25].
In addition to having a clear physical interpretation, W ex (β, β b ) is easily calculated as a simple numerical integral of equilibrium Gibbs-Boltzmann distributions for two temperatures. By contrast, to establish the nonmonotonicity of a 2 , the criterion of Lu and Raz [1], one must first find the left eigenfunction u 2 by solving the boundary-value problem associated with the adjoint Fokker-Planck operator.  Figure 1A.

V. DISCUSSION
The anomalous relaxation process known as the Mpemba effect is defined by a non-monotonic dependence of relaxation time on initial temperature. Lu and Raz [1] showed that an equivalent criterion is the nonmonotonicity of the a 2 projection coefficient derived from an associated Fokker-Planck equation. In this Brief Research Report, we have shown that, for a 1D potential U (x) with a metastable and a stable minimum, the Mpemba effect can be viewed as a simple two-stage relaxation in the function space of all admissible probability densities. In the fast Stage 1, the system relaxes to a local equilibrium. In the slow Stage 2, the populations in the two coarse-grained states equilibrate. In such a situation, we have shown that the Mpemba effect is associated with a non-monotonic temperature dependence of the maximum extractable work of the local equilibrium stage reached at the end of Stage 1.Relative to the a 2 coefficient, extractable work is a much more physical quantity that is also much easier to calculate.
The physical picture offered here, for a double-well potential, meets our goal: We can relate the existence of the Mpemba effect to a non-monotonicity of the extractable work. However, we have not carefully characterized the range of validity of the approximations used in our analysis. For example, in writing Equation (9), we assume that the fraction of initial systems that start in either state (x < 0 or x > 0) is preserved after the initial fast transient. In fact, even during the brief transient, calculating the fraction of systems in each region is subtle, a point emphasized by van Kampen [32] in a careful study that would be the starting point for a more detailed theoretical investigation.
Although our arguments assume a 1D potential with two local states, they generalize easily to many dimensions and many local states. In such cases, the state vector has a large number of dimensions, and solving the Fokker-Planck equation or even calculating its eigenfunctions is difficult. But calculating the extractable work remains easy. Of course, our arguments do not imply that the Mpemba effect can occur only in potentials with metastable states and leave open the possibility for other scenarios. APPENDIX 1. Monotonicity of Kullback-Leibler divergence along G leq . The Kullback-Leibler divergence [33] can be written in terms of Equation (9) as In the second line, we omit the corresponding a R terms. In the third line, a * L ≡ 0 −∞ dx g β b (x) and a * R ≡ ∞ 0 dx g β b (x). In the fourth line, the vectors represent two-state probability distributions. Note that in the "short Stage 1" approximation of Equation (10) We then investigate the monotonicity of D KL a L a R , a * L a * R by differentiating: which is positive for a L > a * L and negative for a L < a * L . (Recall that a L + a R = a * L + a * R = 1.) Thus, D KL ρ leq , g is monotonic in a L on either side of equilibrium. 2. Proof of Equation (15). The relationship is well known [34] and holds for any distribution, including ones describing local equilibrium. Below, to simplify notation, we write ρ leq for ρ leq β,β b and g for g β b .
which is equivalent to Equation (15).