Accidental Origin of the Mpemba Effect in Freezing: Reversals in the Freezing Order are Stochastic. Experiment and Theory

1. Introduction

Nearly 60 years ago, the surprising claim that “the water that starts hotter may freeze first” was dubbed the Mpemba effect [1]; the classical Mpemba effect implies that, under certain conditions, any (or almost any) vial of hot water freezes before any (or almost any) vial of cool water. Similar counterintuitive phenomena, seemingly defying basic physics and chemistry, have since been observed in liquids, polymers, clathrate hydrates, colloids, and even quantum systems [2,3,4,5,6,7,8,9,10,11,12,13,14,15] and discussed in many theoretical studies [8,9,11,13,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. This reflects burgeoning research on this effect, which is of interest to chemists in various areas, highlights that the underlying physicochemical mechanisms remain elusive [31], and that even the seemingly simple process of water freezing is not fully understood.

Our water freezing experiments demonstrate that the Mpemba effect in freezing occurs stochastically, accidentally, and only when the freezer temperature is very close to that of ice nucleation. As the theory shows, all this is due to the stochastic uncertainty in the timing of ice nucleation, this first-order phase transition, and does not require any hypothesis about any special freezing mechanism or “water memory”.

2. Results

2.1. Experiment

To clarify the basis of the in-water Mpemba paradox, we conducted experiments replicating previously reported studies on water freezing [1,2,3,4,7,15] under carefully controlled conditions (see Materials and Methods and Supplementary Materials 1).

Consistent with seminal studies [1,3], we focused [34] on the onset of freezing, rather than on the complete freezing of ice or the cooling of water to 0 °C. The latter is suggested by careful experiments on water cooling to 0 ℃ without freezing [4], which show that hot water does not cool faster than cold water. This allowed us to neglect water evaporation [35,36] (especially since in our, as well as in some other [1,2,3,4] experiments, the water samples were kept in closed vials), and convection [37,38] as causes of the Mpemba effect, and focus on the phase transition that apparently drives this effect in various systems, including milk [1], water [1,2,15], clathrate hydrates [5], and polymers [6].

Figure 1 shows the dynamics of cooling and freezing of water samples (placed in closed vials) at different thermostat temperatures T_th.

In each experiment, the four “cool” water samples started at +20 °C, and the four “hot” ones started at +70 °C; all of these were positioned in the freezer simultaneously. When cooling, all samples of the initially hot water reached 0 °C approximately 500 seconds later than all samples of the initially cool water (which agrees with Ref. 4); with further cooling at T_th>-7 °C, some samples of the initially hot water froze earlier, and some later, than the samples of the initially cool water, but at T_th<-8 °C all the samples of the initially hot water froze later than all the samples of the initially cool water.

The main panel of Figure 1a illustrates the full temperature evolution in one experiment with 8 water samples at a thermostat temperature T_th above -7 °C, while the main panel of Figure 1b illustrates the full temperature evolution at a thermostat temperature T_th below -8 °C. The insets in panels 1(a) and 1(b) show the low-temperature dynamics of the behavior of water samples in two other experiments (among many others; see Table 1a and Table 1b) conducted at thermostat temperatures above -7 °C (a) and below -8 °C (b).

Table 1a and 1b present detailed data on the cooling and freezing of all initially hot and cool water samples at different thermostat temperatures T_th.

At the freezer temperature T_th≥-7.2 °C (panel 1(a) and Table 1a), one can see a large variation in freezing onset times indicated by sharp temperature increases [39]. For the “initially cool” (+20 °C) and “initially hot” (+70 °C) water samples, the ranges of the onset times do overlap, so that at this relatively high T_th (-7.2 °C and higher) the Mpemba effect is observed, but it is stochastic– it is observed for some pairs of samples, but not for the others.

After each experiment (for T_th≥-7.2 °C, we conducted 9 experiments with 4 “initially hot” and 4 “ initially cool” water samples; see Table 1a) we did 16 comparisons of these 4 “hot” and 4 “cool” samples; the comparisons were not entirely independent, as we were dealing with only 8 independent figures; however, the ratio of the sums of the resulting figures—shown in the last two columns—provides a simple, general overview of the relationship between the “hot freezes before cool” and the “cool freezes before hot” cases. In total, we had 4×4×9 = 144 comparisons. The Mpemba effect (the “hot” sample freezes first) occurred in 39% of the cases. For example, the main panel of Figure 1a shows that “hot” samples 1h – 4h begin to freeze earlier than “cool” samples 1c, 3c, but later than the “cool” sample 2c.

Conversely, at T_th<-8 °C (panel 1(b) and Table 1b), one can see a small variation in freezing onset times. For the +20 °C and +70 °C samples, the onset times practically do not overlap, so that initially cool water virtually always (in 170 cases out of 176) freezes first (and below -10 °C, initially cool water always freezes first), so that at the low T_th (below -8 °C) the Mpemba effect is practically not observed. This is proved by 4×4×11 = 176 comparisons of 4 “initially hot” and 4 “ initially cool” water samples (the caveat made earlier—that these comparisons are not entirely independent—does not alter this general conclusion).

Thus, our water freezing experiments highlight (i) the stochastic nature of the onset of freezing, especially at T_th>-7 °C, and (ii) the significant role of the “not very cold” (-7.2 °C and higher) freezer temperature T_th in inducing this stochastic Mpemba effect.

Previous experimental papers have sometimes hinted at the stochastic, accidental nature of the Mpemba effect, using cautious phrases such as “the water that starts hotter may freeze first” [1] or “hot water may well begin to freeze quicker than cold” [13], but no comprehensive investigation into this critical aspect of the Mpemba effect related to ice nucleation has been reported.

Here, we will outline and try to explain the conditions under which the stochasticity in ice nucleation time sometimes leads to the observation of the Mpemba effect.

It should be noted that the statement about the stochasticity of the Mpemba effect applies only to events involving first-order phase transitions (see below). Systems without these transitions, such as colloids [8,11], some quantum [13,14,29,30,32] and mechanical [40,41] systems, can exhibit a deterministic rather than stochastic Mpemba effect.

2.2. Theory

The freezer temperature T_th impacts both freezing and cooling rates.

The latter, according to Newton’s law of thermal conductivity, can be described as

dT_t/dt = −(T_t − T_th);

(1)

here,

is the thermal conductivity coefficient attributed to the surroundings of the water sample (having temperature T_t), and t is time. This equation implies that d[ln(T_t − T_th)] = −dt; so,

T_t − T_th = (T₁ − T_th)exp(−t),

(2)

where T₁ > T_th represents the temperature of the sample at the initial moment (t₁=0) of the cooling. Thus, the time required for the sample to cool from T₁ to T_t is

$$t_{\mathit{th}}\left(T_1,T_t\right)={\displaystyle \frac{1}{}}\mathit{ln}{\displaystyle \frac{T_1-T_{\mathit{th}}}{T_t-T_{\mathit{th}}}};$$

(3)

here, ~0.002 s⁻¹ corresponds to cooling by a freezer set at T_th of about −5 − 15 ℃ (as in our experiments and in seminal works [1,2], where cooling from T₁ ≈+80 ℃ to T_t ≈+20 ℃ took about 10 minutes).

This ~0.002 s⁻¹ estimate will be used throughout our study for both “initially hot” and “initially cool” water samples without fitting to each cooling trace, since the exact value of

is not crucial for understanding the nature of the Mpemba effect.

The water freezing can start only when T_t becomes sufficiently [2,39,42] below T₀≡273 ^oК, i.e., 0 ℃.

According to the conventional Gibbs transition state theory for first-order phase transitions [39,43,44,45] (like water freezing), the time of appearance of an ice nucleus (corresponding to the most unstable, i.e., transition state during ice formation) around one given H₂O molecule in pure, still water can be estimated as

$$t_T^{\left(1\right)}~\tau ·\mathit{exp}\left({\displaystyle \frac{{+G}_{\mathit{nucl}}^{\#}\left(T\right)}{k_{B}T}}\right).$$

(4)

Here, $G_{\mathit{nucl}}^{\#}\left(T\right)$ is the free energy of the ice nucleus at temperature T<T₀; k_B is the Boltzmann constant; τ is the typical time for one H₂O molecule to attach to ice. The value of τ is no less than τ₀~10^-12 s (the typical time [45] of thermal vibrations at 0 ℃); but, in fact, τ is estimated (see [45], chs. 3.2, 8.2) from the difference between the rates of H₂O molecule attachment to and detachment from ice. The final estimate of τ also includes the “Zeldovich factor”, i.e., the number of steps spent to cross the top of the activating barrier, and the number of possible parallel nucleus growth pathways [42,46]. With only the dominant terms considered [42,45,46], τ~10^-7 s at “normal” sub-zero temperatures ($-{1^{\mathrm{o}}\mathrm{to}-20}^{\mathrm{o}}С$). This value of τ depends on temperature, but much less than $\mathit{exp}\left({\displaystyle \frac{{+G}_{\mathit{nucl}}^{\#}\left(T\right)}{k_{B}T}}\right)$, which rapidly tends toward infinity as the temperature T approaches 0 ℃ (so that the water freezing temperature must be strictly below 0 ℃ for kinetic reasons [39,42]).

The free energy $G_{\mathit{nucl}}^{\#}$ is very large for a 3-dimensional (3D) ice nucleus in bulk water [42], where ice, for kinetic reasons, arises only below -35 ℃ [47].

However, the 2-dimensional (2D) ice nucleus on an ice-binding surface has a smaller $G_{\mathit{nucl}}^{\#}\left(T\right)$; thus, ice can arise [42] on an ice-binding border of water at negative temperatures closer to 0 ℃. Then $G_{\mathit{nucl}}^{\#}\left(T\right)$ can be estimated as $\approx 4B_2\left({\displaystyle \frac{B_2}{-∆\mu \left(T\right)}}\right)$ [42,45], where Δμ(T) is the chemical potential of H₂O in ice minus that in water (so that ∆μ(T₀)=0, and ∆μ(T)<0 at T<T₀, where ${\displaystyle \frac{∆\mu \left(T\right)}{k_{B}T}}\approx {\displaystyle \frac{T-T_0}{{100}^{\mathrm{o}}}}$, see [45], with data taken from [49]), and B₂ is the additional free energy of an ice molecule at the ice/water interface; B₃≈0.85 k_BT₀ for a 3D ice/water interface [49] (see also [45], with data from [48]); this estimate seems to remain more or less valid for the perimeter of a 2D ice layer (i.e., B₂≈k_BT₀) [50]. With these ∆μ(T) and B₂ estimates, and T<T₀ (but close to T₀), we have [42,46,50,51]

$${\displaystyle \frac{G_{\mathit{nucl}}^{\#}\left(T\right)}{k_{B}T}}={\displaystyle \frac{T_{*}}{T_0-T}}$$

(5)

where T_*≈400 ℃ at B₂≈k_BT₀. This estimate captures the main: an enormous growth of $G_{\mathit{nucl}}^{\#}\left(T\right)$ when T approaches the ice-and-water coexistence temperature, 0 ℃.

Since at a given temperature T<T₀, the time of appearance of an ice nucleus around one given on-border H₂O molecule is determined by Eqs. (4), (5), the time of appearance of an ice nucleus around one of the N on-border H₂O molecules existing on the surface of the vessel (where, in the absence of dust, ice can only appear at sub-zero temperatures) can be estimated as

$$t_T^{\left(N\right)}~{\displaystyle \frac{\tau }{N}}·\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T}}\right).$$

(6)

With the number of on-border waters N ~ 10¹⁵ in a ~1 mL test tube (since a single water molecule with a diameter of ≈3 Å covers ≈10^–15 cm² of the test tube wall, and different water molecules can be roughly considered as independent potential nucleation sites), and (see above) τ~10^-7 s approximately, one can see from this equation that the value of $t_T^{\left(N\right)}$ amounts—theoretically—to minutes at the water temperature T = T₀ − 7 °C, days at T = T₀ − 6 °C, and ages at T = T₀ − 5 °C.

Thus, freezing of pure water is only possible at water temperatures T below T_freeze≈−6 ℃. It is worth noting that changing the values of N and τ by an order of magnitude each leads to roughly the same estimates of $t_T^{\left(N\right)}$ (minutes - days - ages at ≈-7, ≈-6, ≈-5 °C) when the fitting parameter T_* is varied by only 10%.

The probability of appearance of an ice nucleus around one given on-border H₂O molecule during the above-specified time τ at temperature T<T₀ is $p_{T<T_0,\tau }^{\left(+1\right)}\equiv {\displaystyle \frac{\tau }{t_T^{\left(1\right)}}}=\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T}}\right)$, while $p_{T>T_0,\tau }^{\left(+1\right)}=0$ at T≥T₀. Thus,

$$p_{T<T_0,\mathit{dt}}^{\left(+1\right)}={\displaystyle \frac{\mathit{dt}}{\tau }}·\mathit{exp}\left({\displaystyle \frac{-T_{*}}{T_0-T}}\right)$$

(7)

is the probability of the nucleus appearing around a given on-border H₂O molecule during a short time dt, and $p_{T,\mathit{dt}}^{\left(-1\right)}=1-p_{T,\mathit{dt}}^{\left(+1\right)}$ is the probability that a nucleus does not appear during this time dt around this on-border water molecule. So, the value ${p_{T<T_0,\mathit{dt}}^{\left(-N\right)}=\left[1-p_{T<T_0,\mathit{dt}}^{\left(+1\right)}\right]}^N\approx \exp \left(-N{p}_{T<T_0,\mathit{dt}}^{\left(+1\right)}\right)=\mathit{exp}\left[-\mathit{dt}{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T}}\right)\right]$ is the probability that, at T<T₀, in a system of N on-border waters not a single nucleus arises during the time dt. Then ${\displaystyle \prod _{i=1}^m}{p}_{T_i<T_0,d{t}_i}^{\left(-N\right)}=$ $\mathit{exp}\left[-{\sum }_{i=1}^{m}d{t}_i{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_i}}\right)\right]$ is the probability that not a single nucleus arises in this system during the time interval dt₁+…+dt_m when all of temperatures T₁,…,T_m are below T₀.

A freezer set at T_th<T₀ needs (see Eq. (3)) the time ${t_0=t}_{\mathit{th}}(T_1,T_0)$ for cooling water from T₁ to T₀, only after which ice can arise. At time above t₀ (when T_t falls below T₀),

$$S_{t^{*}>t_0}^{-}=\mathit{exp}\left[-{\int }_{t_0}^{t^{*}}\mathit{dt}{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_t}}\right)\right]$$

(8)

is the probability that no nucleus arises among N of on-border waters during time t₀t^* (note that $S_{{t\le t}_0}^{-}\equiv 1$, and $S_{t^{*}{>t}_0}^{-}$ tends to 0 at t^* → ∞ when all of T_t<T₀ at t₀≤t≤t*). Then

$$S_{t^{*}>t_0}^{+}\equiv 1-S_{t^{*}>t_0}^{-}=1-\mathit{exp}\left[-{\int }_{t_0}^{t^{*}}\mathit{dt}{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_t}}\right)\right]$$

(9)

is the probability that the first nucleus arises in the on-border water during the time interval t₀?t* (and $S_{{t\le t}_0}^{+}=0$). So, the probability of ice nucleation during the short time interval t*?(t*+dt*) is $d{t}^{*}\left\{{\displaystyle \frac{d}{d{t}^{*}}}{S}_{t^{*}>t_0}^{+}\right\}$, which equals to $\left\{\mathit{exp}\left[-{\int }_{t_0}^{t^{*}}\mathit{dt}{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_t}}\right)\right]·\left[{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_{t^{*}}}}\right)\right]\right\}·d{t}^{*}$. Thus,

$$P_{t^{*}>t_0}^{+}{\displaystyle \frac{d}{d{t}^{*}}}{S}_{t^{*}>t_0}^{+}=\mathit{exp}\left[-{\int }_{t_0}^{t^{*}}\mathit{dt}{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_t}}\right)\right]·\left[{\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_t}}\right)\right]$$

(10)

is the momentary probability of ice nucleation at t^* > t₀ (note: $P_{t^{*}\le t_0}^{+}$=0, ${\int }_{t_0}^{\infty }{P}_{t^{*}>t_0}^{+}d{t}^{*}$=1, while Eq. (3) gives the time t₀ = t_th(T₁,T₀) needed to cool water to the temperature T₀, and the time $t_{\mathit{th}}\left(T_1,T_{t^{*}}\right)$ needed to cool an ice-free water to the temperature $T_{t^{*}}$).

Figure 2 shows the computed cooling dynamics of water and computed probabilities $P_{t^{*}}^{+}$ of ice nucleation as functions of the cooling time t^*. These computations use Eqs. (2, 10) with N=10¹⁵, τ=10^-7 s, T_* = 400 ℃, =0.002 s⁻¹ for cooling from the initial temperatures T₁=+20 ℃ (panels 2 (a, c)) and T₁=+70 ℃ (panels 2 (b, d)) at various thermostat temperatures T_th.

Thus, as illustrated by Figure 2 (in accordance with Figure 1 and the Graphical Abstract) for both initial temperatures, the expected ranges of ice nucleation times are quite wide (>>1000 sec, see caps on the tops of $P_{t^{*}}^{+}$ peaks) and overlap at freezer temperatures above −7 ^oC, though at lower (below −8 ^oC and especially −10 ^oC) freezer temperatures, they are narrow (less than ≈500-second delay in the cooling of hot water) and do not overlap.

This theoretical result agrees with experimental observations across water samples cooled at different freezer temperatures (Figure 1a,b), supporting the conclusion that the “Mpemba effect” manifests itself—but only in a fraction of experiments—as a stochastic event, and only at relatively high freezer temperatures T_th above ≈−7 ^oC.

3. Discussion

Figure 1 and Table 1a, 1b clearly demonstrate the following striking features:

(i) water (in a closed vial) freezes not at 0 ℃ but (for kinetic reasons [39,42,45,46,50,51]) only below ≈-6 ℃, i.e., in a super-cooled state [2,39], and even then, starts to freeze not immediately;

(ii) the typical time of the onset of water freezing at cooling is very large at T_th ≈ -6 ÷ -7℃ [2,39,45]; this time decreases significantly (approximately tenfold) with decreasing thermostat temperature T_th from -6.3 ÷ -6.9 to -8.6 ÷ -8.7 ℃, i.e., by only ≈ 2°; this accords to the above (Figure 2(c, d)) estimates;

(iii) at the freezer temperature T_th above -7.2 ℃, freezing occurs at water temperature T_fr that is only slightly, typically by ≲0.5 ℃, above this T_th, while at the lower temperature T_th, of -8 ÷ -12 °C, the freezing occurs at water temperature T_fr which is close to ≈-6.5 ℃, and thus noticeably, by about 2° ÷ 5°, above this low T_th;

(iv) at the thermostat temperature T_th ≈ -8.6 ℃ and below, the scattering of the freezing onset times is relatively small – for all samples of both initially cool and initially hot water this scattering is less than 500 seconds (that is, less than the typical delay in the cooling of initially hot water), but at the thermostat temperature T_th ≳ -7 ℃, there is a huge (~10000 seconds) scattering of freezing onset times both for initially cool and initially hot water (which is much beyond the typical ≈500 second delay in cooling of initially hot water), so that initially hot water sometimes freezes faster and sometimes slower than initially cool water; this also corresponds to what is seen in Figure 2.

To grasp why higher freezer temperatures T_th result in wider freezing time ranges, while lower T_th yield narrower freezing time scattering, we can find an approximate solution to Eq. (10) by simplifying the consideration of the cooling & freezing process by dividing it into two stages:

(a) cooling, when the temperature drops from T₁ to T_fr, which takes time t_c = t_th(T₁,T_fr), and

(b) freezing onset at the constant temperature T_fr < T₀.

It should be noted that the purpose of this simplification is to provide a general understanding of the cooling & freezing process, rather than a precise description of its details. The approximation used is justified by the fact that the change in water temperature on the timescales at which ice formation begins is quite small (see Figure 2), and it allows for the use of deliberately approximate, rounded numerical parameters.

During cooling with a freezer set at the temperature T_th, the most probable ice-nucleation temperature T_fr(T_th) corresponds to the minimum (by T_fr) of the sum of the cooling time, t_th(T₁,T_fr), see Eq. (3), and the following (after reaching T_fr) mean freezing onset time ${〈t_{\mathrm{f}}〉}_{T_{\mathit{fr}}}={\int }_0^{\infty }t·P_{t>0}^{+}\mathit{dt}$. Denoting the constant, at T_t = T_fr = const, value ${\displaystyle \frac{N}{\tau }}\mathit{exp}\left({\displaystyle \frac{{-T}_{*}}{T_0-T_{\mathit{fr}}}}\right)$ as $A_{T_{\mathit{fr}}}$, we have – from Eq. (10) – $P_{t>0}^{+}=A_{T_{\mathit{fr}}}\exp (-A_{T_{\mathit{fr}}}t)$, and

$${〈t_f〉}_{T_{\mathit{fr}}}={\int }_0^{\mathrm{\infty }}t·A_{T_{\mathit{fr}}}\mathit{exp}\left(-A_{T_{\mathit{fr}}}t\right)·\mathit{dt}={\displaystyle \frac{1}{A_{T_{\mathit{fr}}}}}\equiv {\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right),$$

(11)

which drastically grows when T_fr approaches T₀.

The minimum (by T_fr) of the sum $t_{\mathit{th}}\left(T_1,T_{\mathit{fr}}\right)+{〈t_f〉}_{T_{\mathit{fr}}}$= ${\displaystyle \frac{1}{}}\ln {\displaystyle \frac{T_1-T_{\mathit{th}}}{T_{\mathit{fr}}-T_{\mathit{th}}}}+{\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right)$ is achieved when

$${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}={\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right){\displaystyle \frac{1}{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$$}\right.}}.$$

(12)

The solution of Eq. (12) gives T_fr as a function of T_th (see Figure 3).

This solution has two regimes: (A) when ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}>1$, and (B) when ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}<1$, separated by a point where ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}=1={\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right)$. At this point (marked by a green closed circle in Figure 3), $T_{\mathit{fr}}-{T_0=-T}_{*}/\mathit{ln}\left({\displaystyle \frac{N}{\tau }}\right)$, and $T_{\mathit{th}}=T_{\mathit{fr}}{-T}_{*}/{\left[\mathit{ln}\left({\displaystyle \frac{N}{\tau }}\right)\right]}^2$.

The solution of Eq. (12) has two regimes:

(A) when ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}>1$; this corresponds to the case where T_fr is very close to T_th, and the cooling & freezing time is determined primarily by the freezing onset time;

(B) when ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}<1$; in this case T_fr is not close to T_th, and the cooling & freezing time is determined primarily by the cooling.

These regimes are separated by a point (marked by a circle) where ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}=1={\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right)$. In this point, $T_{\mathit{fr}}={T_0-T}_{*}/\mathit{ln}\left({\displaystyle \frac{N}{\tau }}\right)$, and $T_{\mathit{th}}=T_{\mathit{fr}}{-T}_{*}/{\left[\mathit{ln}\left({\displaystyle \frac{N}{\tau }}\right)\right]}^2$.

With the above values of N=10¹⁵, τ=10^-7 s, and =0.002 s⁻¹, T_* = 400 ℃, one has T_th≈−7.16 ℃ and T_fr(T_th)≈−7.03 ℃ (see green line in Figure 3). With slightly, by 10%, different numerical parameters (=0.0022 s⁻¹, T_* = 360 ℃, while N, τ remain as specified above), we obtain T_th≈−6.45 ℃ and T_fr(T_th)≈−6.34 ℃ (see blue line in Figure 3). These resulting 10% changes in the values of T_th and T_fr have practically no effect—as we shall see below—on the explanation of why onset times and their variabilities are immense at higher freezer temperatures and comparatively small at lower temperatures; however, the freezing onset temperatures T_fr calculated with =0.0022 s⁻¹, T_* = 360 ℃ seem to be closer to the experimentally measured values (see Figure 1 and Footnotes 2 to both Tables). .

Having estimated the dependence of T_fr on T_th, we can find – from Eqs. (3), (11) – the estimate of t_c(T_th), i.e., of the T₁ → T_fr(T_th) cooling time, and the estimate of ${〈t_f〉}_{T_{\mathit{fr}}\left(T_{\mathit{th}}\right)}$, i.e., of the mean time of the onset of the subsequent freezing at the freezer temperature T_th.

The total mean time t_tot = $t_{\mathrm{c}}\left(T_{\mathit{th}}\right)+{〈t_f〉}_{T_{\mathit{fr}}\left(T_{\mathit{th}}\right)}$ of the cooling & freezing process (i.e., the average freezing onset time) is presented in Figure 4.

The value of t_c depends on T_th only slightly, whereas the ${〈t_f〉}_{T_{\mathit{fr}}}\mathrm{dependence\; on\; it}$ is strong. Therefore, a sharp rise in $t_{\mathrm{c}}\left(T_{\mathit{th}}\right)+{〈t_f〉}_{T_{\mathit{fr}}\left(T_{\mathit{th}}\right)}$ at T_th > −7 ℃ depends primarily on a sharp rise of the ${〈t_f〉}_{T_{\mathit{fr}}}$ values (given by Eq. (11)) at higher T_th.

Since the value of T_fr(T_th) − T₀ changes – in our experiments – by a factor of two only (-5 to -10 ℃ at the most), while the T_fr − T_th may vary by orders of magnitude (from ~0.01 ℃ (see lines in Figure 3 at T_th close to -6 ℃) to ~5 ℃ (see lines in Figure 3 at T_th close to -12 ℃)), T_fr − T_th has a much stronger impact, than T₀ − T_fr, on the value of ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}$.

At ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}>1$, which is the case when T_th > −7 ℃, and then T_fr − T_th is very small (see Figure 3), ${\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathrm{fr}}}}\right)$ exceeds ${\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$$}\right.}$ and grows exponentially with increasing T_th and T_fr. This is the case (A), when the expected freezing onset time ${〈t_f〉}_{T_{\mathit{fr}}\left(T_{\mathit{th}}\right)}$ can be exponentially large.

On the contrary, in the case (B), when ${\displaystyle \frac{[{\left(T_0-T_{\mathit{fr}}\right)/T_{*}]}^2}{\left(T_{\mathit{fr}}-T_{\mathit{th}}\right)/T_{*}}}<1$ since T_fr − T_th is not as small (which happens when T_th < −7.2 ℃), ${\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right)<{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$$}\right.}$, and then the expected freezing onset time ${〈t_f〉}_{T_{\mathit{fr}}\left(T_{\mathit{th}}\right)}$ is relatively small – less than the ~500-second delay in cooling the initially hot water.

Now we can estimate the variance of the expected time of the cooling & freezing process. As known, the expected root mean square deviation (RMSD_X) of experimental X values from their mean 〈X〉 is $\sqrt{〈X^2〉-{〈X〉}^2}$, where 〈X²〉 is the mean square of X.

The mean square of the freezing onset time ${\left.t_f\right|}_{T_{\mathit{fr}}}$ is

$${〈t_f^2〉}_{T_{\mathit{fr}}}={\int }_0^{\mathrm{\infty }}{t}^2·P_{t>0}^{+}·\mathit{dt}={\int }_0^{\mathrm{\infty }}{t}^2·d[-\mathit{exp}(-A_{T_{\mathit{fr}}}t\left)\right]={\displaystyle \frac{2}{A_{T_{\mathit{fr}}}^2}}\equiv 2{〈t_f〉}_{T_{\mathit{fr}}}^2,$$

(13)

so that the expected variance of the experimental freezing times is

$${\mathit{RMSD}}_{t_f}=\sqrt{{〈t_f^2〉}_{T_{\mathit{fr}}}-{〈t_f〉}_{T_{\mathit{fr}}}^2}={〈t_f〉}_{T_{\mathit{fr}}}\equiv {\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right),$$

(14)

which is just equal (see Eq. (11)) to the mean freezing onset time ${〈t_f〉}_{T_{\mathit{fr}}}$ – and both drastically grow when T_th (and hence, T_fr) is above −7.2 ℃, see Figure 5 (which, like Figure 4, demonstrates satisfactory agreement between the theoretical curves and the experimental points, especially at =0.0022 s⁻¹, T_* = 360 ℃, see panels 4(b), 5(b)). It should be noted that the fact that the standard deviation of the freezing onset time is equal to the mean freezing onset time is generally typical of the simplest Poisson process (when the expected number of events is zero, just like the number of nucleation events prior to a first-order phase transition).

At the same time, the expected variance of the cooling time t_c(T_fr) is small, despite some (small, because ice-nucleating clusters of waters are large [45,50] around 0 ℃) stochastic deviations of local temperatures from the mean water temperature T_fr [52]. So, the variance of the cooling time t_c can be neglected – unlike the variances in the freezing onset times t_f that can be large and therefore significant.

Thus, the resulting estimate of the time of the whole cooling & freezing process is

$$\left[t_c\left(T_{\mathit{fr}}\right)+{〈t_f〉}_{T_{\mathit{fr}}}\right]{\mathit{RMSD}}_{t_f}\approx \left[{\displaystyle \frac{1}{}}·\mathit{ln}{\displaystyle \frac{T_1-T_{\mathit{th}}}{T_{\mathit{fr}}-T_{\mathit{th}}}}+{\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right)\right]{\displaystyle \frac{\tau }{N}}\mathit{exp}\left({\displaystyle \frac{T_{*}}{T_0-T_{\mathit{fr}}}}\right),$$

(15)

where the exponential terms are dominating when T_th (and hence, T_fr(T_th)) is above −7 ℃.

As a result, one expects (and does observe in experiment, see Figure 1 and Table 1a, 1b) long, highly variable ice-nucleation times when T_th is above −7 ℃, and short, low-variable ice-nucleation times when T_th is below −8 ℃ and especially −10 ℃.

4. Materials and Methods

In each of our many experiments, we used four “initially cool” and four “initially hot” water samples, all placed in the freezer simultaneously. This simultaneous study of water freezing in many vessels allowed us not only to compare freezing events occurring in vessels with different initial temperatures (which has typically been done in other studies of the Mpemba effect), but also to compare processes occurring in vessels with the same initial temperatures (which has not been done in other studies of the Mpemba effect).

The experiments were conducted in the basement of a separate building, where all vibration sources were turned off to prevent their impact on ice freezing. For cooling, a Tuvio FCS15HV1 chest freezer was used; its horizontal lid (see Supplementary Materials 1, Figure S1) minimizes temperature fluctuations at its opening, in contrast to freezers with a vertical door. This freezer allows one to set the temperature from +13 ℃ to −30 ℃ with an accuracy of about 1 ℃. We always used temperatures above −15 ℃ (this was the temperature at the chest freezer bottom). Except for the heat exchanger, the freezer was well insulated. To further stabilize the temperature, about 30 kg of frozen water in big plastic bottles and about 3.5 kg of frozen water in small plastic bottles was placed in the freezer.

To mitigate temperature fluctuations, we constructed an “inner chamber” within the freezer. This chamber is a foam plastic box (see Supplementary Materials 1, Figure S1 and Figure S2) that housed eleven small plastic bottles, each containing 310 mL of a 15% NaCl solution. Before each experiment, these bottles were pre-cooled to about −12 ℃. The gradual thawing of this salt water effectively stabilized the temperature inside the chamber, reducing its fluctuations to less than 0.5 ℃ a day, in contrast to the 2 ℃ fluctuations observed outside the chamber due to the compressor’s cycling.

We initiated each experiment only after the temperature T_th inside the inner chamber, monitored by a thermistor (i.e., a temperature-sensitive resistor connected by a long wire to an external digital thermometer), reached the desired point for this experiment. Within the inner chamber, lightweight wire stand (see Supplementary Materials 1, Figure S1 and Figure S2) was suspended using 18 thin, 100-cm-long wires. Nine thermistors, connected to external electronics, were affixed to the ends of these nine wire pairs. Eight of the thermistors were inserted into eight closed 2-mL test vials through holes in their lids. Each vial was thus suspended by a pair of wires, which minimized vibrations, with a thermistor positioned at the end of the wire pair. Water occupied 80%±1% of the vial volume, and the thermistor measured the temperature in the upper half of the liquid.

The assignment of the same physical vial to the hot/cool groups and to different positions was randomized between runs.

The vials and their contents were weighed with a precision of 10 mg. Four vials containing the initially hot water (+70 ℃) and four containing the initially cool water (+20 ℃) were arranged on the wire stand in an alternating pattern, separated by 4 cm. All vials were placed into the freezer simultaneously. The ninth thermistor was placed in the air between the salt water bottles and the test vials to monitor the ambient temperature.

For temperature measurements, we used nine MF59 100kΩ thermistors, each with a 1-meter-long wire. These thermistors were pre-calibrated using a digital thermometer with a resolution of 0.01 ℃ (LT-300, Termexlab.ru). Eight of these thermistors were used to monitor the water temperature within the test vials, while the ninth measured the air temperature inside the inner chamber containing the small bottles with salt semi-frozen water.

Additionally, we used pre-calibrated digital thermometers (of 0.1 ℃ resolution) with sensors at the ends of the 1-meter-long wires to monitor the temperature inside the inner chamber, the temperature of ambient air, and the temperature in the electrical heater used for preheating the “initially hot” water samples.

These two digital thermometers, DT1 and DT2, were used to monitor the temperature throughout the experiment, especially during the process of water freezing (when the temperature of the water/ice mixture is close to 0 ℃), and after the freezing has been completed (when the temperature of the formed ice approached the temperature inside the inner chamber). These data were used for the final calibration of all thermistors.

We used a custom-built data acquisition module that included a computer running specialized software, which sent 5V pulses sequentially to each of the nine thermistors via the LPT1 parallel port and a custom adapter (see Supplementary Materials 1, Figure S3). Digital voltage readings, representing the temperature, were collected from the Radio Shack 22-168A digital multimeter connected to the RS232 serial port. This multimeter measured the voltage across a single 100 kΩ resistor (shared by all thermistors and connected to ground), the current through which depended on the temperature of the thermistor.

Data regarding voltage and corresponding timestamps were recorded at 43-second intervals (see Supplementary Materials 1, Figure S4); thus, the uncertainty in time measurements was ±22 seconds. The collected voltage data were subsequently converted into temperature data (see Supplementary Materials 1, Figure S5). The maximum expected error in temperature measurements (including the thermistor calibration uncertainty) was ±0.5 ℃. This level of precision is adequate for the objectives of this study, as the temperature increase associated with the onset of freezing exceeds 5 ℃, and the point of our interest was the freezing onset time.

All experiments were conducted using distilled water from the Institute of Protein Research, Russian Academy of Sciences. Before each run, four closed vials were equilibrated at the ambient basement temperature (+15 ÷ +20 ℃), while four others were pre-heated to the temperature of ≈+70 ℃ using small custom-built electric ovens.

The interiors of the test vials and thermistors in the vials were lubricated with “Sylgard 184 silicone elastomer” (without the “Curing Agent”) to adsorb potential foreign ice nucleators from water and to smooth out any micro-roughness on the vial walls.

Upon placing the wire stand with vials into the freezer, a small temporary temperature increase was observed in the upper part of the chamber. However, the air temperature, as monitored by a digital thermometer and the intra-chamber thermistor, returned to the setpoint within 10 minutes—well before freezing onset.

5. Conclusions

We see that the Mpemba effect in water freezing (at least under the given conditions: closed vessels that prevent evaporation, suppress convection, and exclude the entry of alien crystallization centers in water) arises not from any special properties of pure water (as has been hypothesized in many papers on this topic), or of some in-water “motes” [39], but from the general uncertainty of the time of the onset of the first-order phase transition—an uncertainty that becomes especially significant at slightly sub-zero water temperatures.

This fundamental uncertainty (which does not allow one to precisely define the onset time) is well known [52] in the physics of first-order phase transitions, but – oddly enough – it has been never, to our best knowledge, considered in the context of the Mpemba effect in the freezing of liquids, even though it is precisely here that this uncertainty manifests itself perhaps most clearly.

It is this fundamental uncertainty in the freezing onset time (which has been mentioned, but not investigated, in [2,39,53,54]) that produces a kind of randomness in the occurrence of the Mpemba effect, and its poor reproducibility mentioned in many works – which, for many years, “continued to make the Mpemba effect as puzzling as ever” [54] and generated skepticism about its very existence [28].

The large variability in the onset times is due to the specific (Poisson) form of the time-dependent probability of the freezing transitions, $P_{t>0}^{+}~\exp (-A_{T_{\mathit{fr}}}t)$. $P_{t>0}^{+}$ has its maximum at the first moment and then declines. This leads to substantial uncertainty in the onset timing and, at sufficiently high T_th, results in overlap of freezing times for samples with different initial temperatures, thus providing the stochastic Mpemba effect.

The observed stochasticity of the Mpemba effect in freezing of pure water does not exclude other factors [55] that may contribute to the freezing times, e.g., insufficiently clean liquid or vessel walls (we minimized these effects, see Materials and Methods), or ice-nucleating “motes” that can fall in open vials positioned in a freezer. Then, the Mpemba effect, not caused by the liquid itself, may result from some foreign ice nucleators penetrating into the initially hot water [15].

It is important to emphasize that the observed accidental origin of the Mpemba effect in pure water is limited to phenomena such as freezing, which involve first-order phase transitions.

The initiation of freezing is, in essence, a catastrophe [56], or “a single-event phenomenon”: the emergence of a single ice nucleus can trigger the freezing of the entire liquid in the test tube. Therefore, the significant uncertainty in the timing of a stochastic ice nucleation can occasionally demonstrate the Mpemba effect in some experiments. This occurs due to the overlap of broad ranges of ice initiation times for initially hot and initially cool samples. This overlap (see Figure 2c,d)) can mask the delay in cooling of the hotter water sample and leading to the wrong perception that “hot water [deterministically, always] freezes faster than cool water”, or, according to a modern, more general formulation of the Mpemba paradox, that “the state that is initially farthest from its equilibrium state attains the latter at the earliest time” [13].

We do not want to create the impression that our finding of the stochasticity of the Mpemba effect in freezing refutes the very existence of deterministic Mpemba effects altogether. Colloids [8,9,11], as well as simple mechanical systems [40,41], and specifically various oscillators (see Supplementary Materials 2 and [57,58]), exhibit a deterministic Mpemba effect, where “the state that is initially farthest from its equilibrium state (always) attains the latter at the earlier [57] (or at least the same [58]) time”, and where the origin and driving force of this effect (inertia—”kinetic memory”—and a sharper than quadratic increase in potential upon deviation from equilibrium) are clear (as we show in Supplementary Materials 2).

In essence, all discussions regarding the Mpemba effect in various systems boil down to attempts to identify a hidden, stable ordering parameter (or parameters)—e.g., currents—in addition to such familiar quantities as temperature, volume, and pressure—much like what has been done for the Kovacs effect [31,59,60,61]. This effect was observed when the polymer glass—having been brought into a non-equilibrium state, where its macroscopic thermodynamic parameters, nevertheless, coincide with their equilibrium values—continues to evolve and exhibits a temporal “hump” before transitioning to a fully equilibrium state.

The discussions regarding the Mpemba effect now apply particularly to anomalous relaxation behavior of open and isolated nonequilibrium quantum systems [62,63,64,65,66,67,68], whose “resource” is asymmetry - thermal or athermal deviation from symmetry or thermal equilibrium. Most of these studies employ a quantum - theoretic approach based on the resource theory of asymmetry and quantum relative entropy, and bear little relation to first-order thermodynamic phase transitions of the freezing type.

Ironically, the classic Mpemba effect, a phenomenon first found just in the freezing of water sixty years ago, proved to be only stochastic in this system, while other systems?lacking a first-order phase transition (colloids, quantum, etc.)?can exhibit a fully deterministic Mpemba effect.