Anomalous heating in a colloidal system.

We report anomalous heating in a colloidal system, an experimental observation of the inverse Mpemba effect, where for two initial temperatures lower than the temperature of the thermal bath, the colder of the two systems heats up faster when coupled to the same thermal bath. For an overdamped, Brownian colloidal particle moving in a tilted double-well potential, we find a nonmonotonic dependence of the heating times on the initial temperature of the system. Entropic effects make the inverse Mpemba effect generically weaker-harder to observe-than the usual Mpemba effect (anomalous cooling). We also observe a strong version of anomalous heating, where a cold system heats up exponentially faster than systems prepared under slightly different conditions.

Can an initially cold system heat up faster than an initially cool system that is otherwise nominally identical? Naively, one would assume that a slowly heating object relaxes to the temperature of its surroundings exponentially, passing through all the intermediate temperatures. A system that is initially at a cold temperature should then take longer to heat than a system initially at a cool temperature. However, for rapid heating, a system may evolve toward equilibrium so that its intermediate states are not in thermal equilibrium and are not characterized by a unique temperature. In such cases, the possibility of anomalously fast heating has recently been predicted and confirmed in numerical studies of an Ising antiferromagnet (1). Further numerical studies suggest that these effects may be seen in a wide variety of systems, including fluids with inelastic (2–4) and elastic (5, 6) collisions and spin glasses (7).

Although anomalous heating is a recent prediction, an analogous anomaly for cooling and freezing has been noted in observations of water dating back to 350 BC (8). Its first systematic study was done in 1969 by Mpemba and Osborne, who concluded that hot water could begin to freeze in a time shorter than that required for cold water (9). This phenomenon has since been dubbed the Mpemba effect and was followed up with further experiments on water (10–16), accompanied by some controversy, tracing back to the difficulty of obtaining reproducible results (17, 18). Proposed mechanisms for the effect include evaporation (10, 19, 20), convection currents (21–23), dissolved gases and solutes (11, 14, 21), supercooling (12, 13), and hydrogen bonds (24, 25).

In an effort to understand the Mpemba effect in more generic terms, Lu and Raz introduced a theoretical picture that related the effect to the geometry of system dynamics in a state space whose elements are defined by the amplitudes of eigenmodes of linear dynamical systems (1). A fast quench can then lead a system to follow a nonequilibrium path through state space to equilibrium that is shorter than the path traced out by a slowly cooling system. In recent work, we showed that this kind of Mpemba effect is present in a system consisting of a colloidal particle immersed in water and subject to a carefully designed potential (26). From this point of view, the dynamics of cooling and heating obey similar principles, and anomalous heating represents an inverse Mpemba effect. Yet, despite a formal similarity between the cases of heating and cooling (1, 27), anomalous heating has not previously been seen experimentally, neither in systems that exhibited the anomalous cooling effect (10–16, 26, 28–30) nor in any other system.

Here we present experimental evidence for the inverse Mpemba effect. We also observe a strong version (31) of the effect, where, for a carefully chosen initial temperature, a system heats up exponentially faster than systems that were initially at different temperatures. Surprisingly, as we shall see, subtle differences between high- and low-temperature limits generically make the inverse effect more difficult to observe experimentally. Moreover, the mechanism for the inverse effect that we find to be relevant in our experiments does not depend on the presence of metastability (32), which played a crucial role in the forward cases explored in previous experiments.

Experimental Setup

In our experiment, a Brownian particle (silica bead, Ø1.5 m) is subjected to forces exerted by an external potential. These forces are generated by rapidly shifting the focus of an optical tweezer according to a feedback rule (Fig. 1). Such feedback traps based on optical tweezers have been used to create arbitrary, complex energy landscapes having high spatial and temporal resolution (33, 34). Here we create a one-dimensional tilted double-well potential U(x) in a finite domain, with a very low barrier (Fig. 2) (Materials and Methods). The geometric asymmetry in the potential is defined by the parameter . Our setup has steep walls at the domain boundaries corresponding to the maximum force pN /nm applied by the optical tweezers (Materials and Methods). The nearly vertical walls confine particle motion to a box with size nm and asymmetry .

Fig. 1.

Feedback trap based on optical tweezers. (A) Schematic of the feedback trap setup. FI, Faraday isolator; M, mirror; SF, spatial filter; BS, beam splitter (nonpolarizing); L, lens; MO, microscope objective; SC, sample chamber; PBS, polarizing beam splitter; HW, half-wave plate; F, short-pass filter; QPD, quadrant photodiode; DM, dichroic mirror; CS, cover slip; Cam, camera; FPGA, field-programmable gate array. Planes conjugate to the back focal plane of the trapping objective are shown in red dashed lines. (Inset) An image of the cage system consisting of the trapping and detection objectives. (B) One cycle of a feedback trap consists of three stages: (i) measurement of the current position of the particle, (ii) calculation of the force based on the potential (black) to be imposed, and (iii) application of the force by shifting the trapping laser (green). Optical tweezers apply a linear restoring force, where k is the stiffness of the harmonic trap and is the expected trap displacement.

Fig. 2.

Schematic of the energy landscape U(x) used to explore the inverse Mpemba effect, set asymmetrically (α  =  2) within a box with potential walls with finite slope at the domain boundaries.

Experimental Protocol

The inverse Mpemba effect occurs if, for three different temperatures , the equilibration time of an initially cold system at temperature is smaller than the equilibration time of the initially cool system at temperature , when coupled to a common thermal bath at temperature . To study the effect, we impose an instantaneous heating quench in our experiments via a three-step process: 1) prepare the initial state of the system corresponding to the Boltzmann distribution at an initial temperature T0, 2) release a particle at a position sampled from the initial distribution , and 3) record the trajectories of the particle as it relaxes in a bath at temperature . The initial positions are sampled assuming U(x) to have infinite potential walls at the domain boundaries. Once the particle is released into the bath, it is always at the bath temperature. We repeat the quenching protocol times, with each cycle 60 ms long, to create a statistical ensemble of the state of the system at each time step s. The dynamics of the particle after the quench in the potential U(x) can be described by the overdamped Langevin equationwhere γ is the Stokes friction coefficient and η is Gaussian white noise, with and .

Although the initial and final states in our experiment obey Boltzmann distributions at temperatures T0 and , the intermediate states p(x,  t) are not in equilibrium: they are not a Boltzmann distribution for any temperature T. For this reason, instead of trying to define an intermediate effective temperature, we measure the distance between the intermediate state p(x,  t) and the equilibrium state (1, 26). From equivalent alternatives (1), we choose the L1 distance (35) for the analysis of particle trajectories in our experiments. This distance is defined as the integrated absolute value of the difference between the densities p and π,

Here is the frequency estimate of the probability for a measured position x at a time t in the interval , where , with and the number of bins. Similarly, π is the frequency estimate of the Boltzmann distribution at .

We note that other choices of distance functions may have a more physical interpretation. For example, the Kullback–Leibler (KL) divergence connects the nonequilibrium free energy of a system to the entropy produced during the relaxation process (36–41). However, for analyzing our experiments, the L1 distance is simpler to evaluate and thus more convenient than the KL divergence. Nevertheless, we confirm, as predicted in earlier theoretical work (1), that the detection of the inverse Mpemba effect is robust to the choice of distance metrics: whenever we observe it using the L1 distance, we see it also when using the KL divergence ().

Inverse Mpemba Effect in an Asymmetric Potential

To determine how the inverse Mpemba effect depends on the initial temperature of the system, we release the particle in a bath of fixed temperature (Fig. 2). After a particle is released in the bath at t = 0 at a low temperature T0, it moves stochastically in response to thermal fluctuations and potential gradient forces and finally equilibrates with the bath, which is at temperature .

Fig. 3 shows example time traces of evolution in the potential U(x). Fig. 4 shows the measured times to reach equilibrium for systems that start at different initial temperatures. The equilibration is defined to be the time when the distance curve reaches the noise-limited final distance (). As the initial temperature of the system decreases from to , the equilibration time increases monotonically and follows normal heating (). However, for the lower initial temperature range , the equilibration time decreases as the initial state of the system gets colder. Such a behavior corresponds to anomalous heating where a cold system takes less time to heat up than a cool system, i.e., . For lower temperatures (), the equilibration time increases again, exhibiting normal heating. Thus, we observe a sequence of normal, anomalous, and normal regimes for relaxation to thermal equilibrium (Figs. 4 and 5).

Fig. 3.

Dynamic trajectories relaxing to equilibrium at a hot temperature. Ten trajectories of a particle released from the equilibrium distribution at temperatures (black), (blue), and (red) into the hot bath, with the evolving probability density p(x,  t) shown for three times (estimated based on 5,000 trajectories) on a logarithmic time scale. The shaded gray region corresponds to a box size L = 240 nm.

Fig. 4.

Equilibration times for systems at different initial temperatures. Red markers denote initial temperature points whose dynamics are displayed in Fig. 5.

Fig. 5.

L1 distances for systems that heat up in a bath at temperature , starting at the temperatures T0 indicated above each graph. The thin vertical lines indicate the times when systems first reach thermal equilibrium (within noise levels).

Analysis in the High–Bath Temperature Limit

In Fig. 2, the variations in throughout the domain are , implying that dynamics at the bath temperature approximate ordinary diffusion, i.e., Brownian motion or Wiener process. To simplify the analysis of the relaxation trajectories in U(x) at a finite temperature , we model the bath as being at an effectively infinite temperature with no energy barrier. Further, we approximate the walls as being infinitely steep. The particle then freely diffuses in a domain with walls at and .

Approximating the bath as being at a very high temperature and the walls as infinitely steep simplifies the analysis in three ways: 1) the equilibrium state is a uniform distribution; 2) the Fokker-Planck operator is self-adjoint, so that left and right eigenfunctions are identical; and 3) the eigenfunctions have simple analytic expressions.

The free-diffusion picture also allows us to understand more intuitively why the inverse Mpemba occurs in this system. At very low initial temperatures, the tilt in the potential leads to a probability density of initial particle positions that is concentrated in the middle of the domain (i.e., within the lower well depicted in Fig. 2). At the high-temperature equilibrium, it is approximately uniformly distributed throughout the domain. Because the dynamics of the system in the bath are dominated by diffusion, most particles diffuse a distance that is typically . At intermediate temperatures, the initial probability density is split between the two well positions. A significant number of particles in the left well will need to diffuse a distance , implying an equilibration time times longer. Indeed, Fig. 4 shows that the equilibration time at a cold temperature can be roughly four times faster than that at intermediate temperatures (40 ms/10 ms = 4).

In the high-temperature approximation, the Fokker–Planck equation describing the probability density p(x,  t) of particle positions reduces to the heat equation,subject to no-flux boundary conditions at , where is the Fokker–Planck operator for a freely diffusing particle (). Note that this high-noise limit is complementary to but less familiar than the low-noise limit, which leads to metastability phenomena (42). Indeed, the high-noise limit of diffusion processes has recently stimulated wide interest because of its relation to the strong-measurement limit of quantum measurements (43, 44).

In the heat equation limit, the probability density p(x,  t) can be written as an infinite sum of eigenfunctions v of with associated eigenvaluesordered so that . Since the contribution of the eigenfunctions decreases exponentially with k, the term k = 2 dominates at large but finite times. Thus, in that regime, the probability density can be approximated by

where the coefficient a2 depends on the initial temperature T0, as well as on the bath temperature , and represents the dynamics of the second mode amplitude during thermalization. Generally, a2 is a measure of the overlap between the second left eigenfunction of the operator written within squared-brackets in Eq. and the initial state of the system (26),

In the high-temperature limit, the spatial eigenfunctions of the diffusion equation are (45)where is the normalization constant, defined such that with .

Since anomalous heating (inverse Mpemba) is associated with a2 coefficients where for cool and cold temperatures, respectively, a nonmonotonic temperature dependence of the a2 coefficients leads to anomalous heating effects. However, these coefficients are not directly accessible in experiments. Instead, we extract them from decay curves of the distance function defined in Eq. . Measuring differences between the initial distance and noise-limited final distance gives a quantity proportional to (26). We measure this difference as a function of initial temperature.

Fig. 6 shows the nonmonotonic temperature dependence of . The values correlate with the measured equilibration times. To see the agreement of the measured values of for the potential at finite temperature with theoretical predictions based on the potential at a high temperature, we explicitly calculate a2 coefficients using Eqs. and . We fit the data to a single parameter, an overall proportionality constant. The fit leads to 1.48 ±0.03, which agrees to with the calculated value, . Thus, we have very good agreement between the experimental a2 coefficient and its theoretical prediction. Note that although is nonmonotonic with the initial temperature of the system, the L1 distance at t = 0 is a monotonically decreasing function of T0 (). As a consequence, even if the system is initially close to thermal equilibrium, it may not heat up more quickly than a system that starts farther away.

Fig. 6.

Measurements of decay amplitude for different initial temperatures T0. Markers denote experimental measurements, and the solid red line is based on the values calculated from the Fokker-Planck equation in the high-temperature limit. The arrow indicates the temperature at which the strong inverse Mpemba effect occurs. Error bars represent 1 SD and are calculated from the fits.

At initial temperature , where (Fig. 6, blue arrow), the decay is dominated by λ3 and represents an exponential speed-up of the heating process compared to decays at temperatures where . Such a situation corresponds to the strong inverse Mpemba effect (46). The transient decay at the time scale set by the eigenvalue ms disappears, and the system decays instead at a rate ms. In summary, for , the initially cool system lags the initially cold system, and the inverse Mpemba effect is observed. For even lower temperatures , the coefficient changes sign, leading to a low-temperature increase in .

Discussion

Our results give clear experimental evidence for the inverse Mpemba effect in a simple setup. The nonmonotonic dependence of the equilibration time on the initial temperature of the system can be understood through the nonmonotonicity of a2 coefficients. We observed the inverse Mpemba effect for a quench for the case of a heat bath whose average energy greatly exceeded the range of potential variation. We used this feature to model system dynamics in a high–bath temperature limit where the relaxation dynamics are governed by a simple heat diffusion equation. Using numerical integrals of analytic expressions for the eigenfunctions, we obtained the a2 coefficients as a function of initial temperature. We found evidence for the strong inverse Mpemba effect, special temperatures where the systems heat up exponentially faster than those at other initial temperatures.

One interesting feature of our system is that it shows that the inverse Mpemba effect can be observed in a system where metastability does not play an important role. Indeed, although our potential has two wells, we can accurately model dynamics in the bath state using an infinite-temperature approximation, where the potential is flat and the motion is purely diffusive. In this case, the form of the potential is important only for defining low-temperature probability densities for the system state. After the quench, it plays no role. By contrast, most scenarios explored previously for continuous-state systems for the forward Mpemba effect depend on the slow time scales corresponding to barrier hopping at low temperature (32). Note, however, that a recent study predicts the forward Mpemba effect in potentials with a single well (47).

Previous experiments on the forward Mpemba effect showed a clear separation of the time scales determined by the eigenvalues λ2 and λ3 (26). We can offer some insight as to why it is easier to observe the forward Mpemba effect than the inverse effect: when a system relaxes to a bath at temperature , the time scale separation between the decay curves corresponding to λ3 and λ2 depends on the ratio . In particular, in order to measure the values, one fits the part of the decay curve that corresponds to λ2. Thus, the greater the value Λ, the easier the accurate measurement of the values.

For the forward Mpemba effect studied in ref. 26, the system cools from a hot temperature to a cold temperature in a double-well potential, and the ratio Λ of eigenvalues λ3 to λ2 (i.e., ) is . However, for the inverse Mpemba effect studied here, the ratio of eigenvalues is . Thus, is about four times smaller than in the case of heating, implying that the forward effect will be easier to observe than the inverse effect. Indeed, our observations of the inverse effect required an ensemble of 5,000 trajectories to obtain results that are statistically similar to results for the forward case obtained with only 1,000 trajectories.

Are these general features of and , or are they special to our potential? Since the dynamics of correspond to hops over the barrier, we expect that the ratio of eigenvalues λ3 to λ2 depends on the barrier height as , an intuition confirmed by a rigorous analysis in general (48, 49) and by numerical solution of the Fokker–Planck equation for our potential in particular (Fig. 7, red curve). However, for the high-temperature limit, Eq. shows that always (Fig. 7, blue curve). Thus, the ratio of eigenvalues Λ can be much higher in the forward case than in the reverse case, and as a result, the forward effect is generically easier to observe experimentally than the inverse effect.

Fig. 7.

Ratio Λ of eigenvalues λ3 to λ2 of the Fokker–Planck operator as a function of barrier height at the bath temperature . The red curve is for the double-well potential used in ref. 26 for the forward Mpemba experiments, and the blue curve is for the approximately flat potential used in the inverse Mpemba experiments. The hollow red marker denotes the ratio used in the forward Mpemba experiments, which corresponds to . The hollow blue marker denotes the ratio for the flat potential () used in the inverse Mpemba experiments.

In this paper, we offer evidence for anomalous heating in a colloidal system, complementing the more familiar scenario for anomalous cooling. Other memory-dependent relaxation phenomena (50) are worth exploring further. For example, Gal and Raz show that an initial cooling quench followed by a heating quench can speed up heating times exponentially, even in systems that would not otherwise exhibit the inverse Mpemba effect (51). In the Kovacs effect, protocols that drive the system out of equilibrium can produce nonmonotonic relaxations that show crossings in a quantity that monitors distance from equilibrium (7, 52–55). It may also be interesting to apply the feedback trap techniques used here to test experimentally the recent prediction that warming can be generically faster than cooling, for equal positive and negative free-energy deviations from thermodynamic equilibrium (27, 56, 57). The concepts of anomalous relaxation that we explore here for the inverse Mpemba effect may lead to a unified picture where all these effects can be seen as representatives of a broad class of physical phenomena.

Materials and Methods

Our optical tweezers setup is built on a vibration isolation table supporting a home-built microscope. We trap a colloidal particle (silica bead, Ø1.5 m, Bangs Laboratories). A linearly polarized 532-nm laser (Nd:YAG, Coherent Genesis MX STM-series, 1 Watt) is used for trapping and detection (Fig. 1). A small portion (10%) of the laser is separated from the main laser as the detection beam before entering the acoustooptic deflectors (AODs). In our experiment, we use AODs to steer the trap position in the trapping plane, placed at a plane conjugate to the back focal plane of the trapping objective. The details of the experimental setup are described in ref. 33. A feedback scheme is used to create the virtual potentials used in the inverse Mpemba experiments (Fig. 1) (26, 33). In a feedback trap, one 1) observes the position of the particle, 2) calculates the force based on its position in the user-defined potential, and 3) applies that force in each loop at a deterministic time step of s (58, 59).

In our experiment, the force is applied by moving the trap center relative to the bead position. The force generated by the displacement of the trap center is approximated as , where is a proportional feedback constant where t is the relaxation time of the underlying physical potential, and where x is the particle position at time .

In our experiment, a freely diffusing bead is subjected to forces exerted by an external, one-dimensional double-well potential, created by a feedback trap based on optical tweezers (33). We place the potential asymmetrically in the domain aswhere is given by

with a very low barrier , with the Boltzmann constant and the bath temperature (Fig. 2). The position x is measured in units of nm. The geometric asymmetry in the potential is defined by the parameter .

Compared to our previous results based on the newly developed feedback traps (33, 60), we have improved the mechanical stability of the setup by installing the trapping and detection objectives within a cage system (Fig. 1). Thus, mechanical drifts due to temperature changes of the surroundings have reduced effects on the particle position. Nonetheless, we observe drifts in the particle position at rates that are on the order of 1 nm s-1. We thus limit the heating cycle to 0.1 s and correct for the small drift before the next cycle.