Article
pubs.acs.org/JPCB

Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

Decoupling of Translational Diﬀusion from the Viscosity ofSupercooled Water: Role of Translational Jump Diﬀusion
Shivam Dueby,† Vikas Dubey,† and Snehasis Daschakraborty*
Department of Chemistry, Indian Institute of Technology Patna, Bihar 801106, India

Downloaded via INDIAN INST OF TECH PATNA on August 14, 2019 at 10:19:19 (UTC).
See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information*

ABSTRACT: Some experiments have witnessed gradual decoupling of viscosity fromthe translational self-diﬀusion of supercooled water with decreasing temperature. This
indicates the breakdown of the Stokes−Einstein equation in supercooled water. Whilesome theoretical and computer simulation studies indicated the jump translation of
the molecules as a probable origin of the above decoupling, direct quantitative
evidence is still lacking. Through a molecular dynamics (MD) simulation study, alongwith careful consideration of translational jump motion, we have found the mostdeﬁnite proof of increasing relevance of translational jump diﬀusion in the above
decoupling phenomena. By separating the jump-only diﬀusion contribution from theoverall diﬀusion of the water, we obtain the residual diﬀusion coeﬃcient, which
remains strongly coupled to the viscosity of the medium at the supercooled regime.
These new ﬁndings can help to elucidate many experimental studies featuringmolecular transport properties, where strong diﬀusion−viscosity decoupling ispresent.

1. INTRODUCTION
There are intriguing properties of supercooled water, including
a strong decoupling between its viscosity and the diﬀusion of
the molecules. Some experimental studies1−3including thatby Dehaoui et al.4have directly or indirectly revealed anincreasing decoupling of viscosity η from the translational
diﬀusion coeﬃcient D of water upon cooling. This indicates a
breakdown of the Stokes−Einstein (SE) relation (D ∝ T/η) as
the temperature decreases below 2.1 Tg, where Tg is the glasstransition temperature. In contrast, the rotational diﬀusion Dr
remains coupled with η for a wider range of temperature,
which implies the validity of the Stokes−Einstein−Debye
relation at a broader range of temperature.4 Similar decoupling
between D and η was also reported earlier in other molecularglass-forming liquid.5−13 The SE equation is obeyed at a
suﬃciently high temperature but severely breaks down below T≈ 1.3Tg. On the contrary, the rotational diﬀusion of themolecular glass-forming liquid and the medium viscosity
remain hydrodynamically coupled with each other even at
the temperature very close to Tg.7,8,14
Deeply supercooled liquids have spatially heterogeneous
dynamics, which have been conﬁrmed by various experiments5,6,15−18 and computer simulation studies.19−23 Somecomputer simulation studies have indicated that the emerging
spatiotemporal heterogeneities in supercooled water and othersupercooled liquids have a strong connection with theincreasing violation of the SE relation with decreasing
temperature.24−31 The violation of the SE relation insupercooled water is also explained by the crossing of the
Widom line, arising from the liquid−liquid critical point(LLCP).32 In this context, it is important to mention the© XXXX American Chemical Society

Goldstein’s hypothesis on the eﬀect of temperature on themechanism of translation of liquid molecules.33 The theory
states that the liquid molecules translate in the rugged freeenergy landscape via hopping from one free energy well toanother via crossing the saddle points at low temperature.
However, the same molecules diﬀuse freely by the Brownian
motion at a higher temperature when the thermal energy
becomes comparable to or higher than the heights of thebarriers of the rugged free-energy landscape. This hypothesis
was supported by many theoretical and computer simulation
studies,34−40 which indicated the presence of translational
jump of molecules and its connection with the breakdown of
the SE relation in the supercooled liquid. Recent works haveshown that the rotation-assisted translational movement ofsolvent water around a solute induces translational jumpdiﬀusion of a tracer from one solvent cage to another insupercooled water.24
Even though the prior studies have implied the pivotal roleof the translational jump of molecules for the breakdown of theSE relation in supercooled water, a more direct quantitative
connection in between the two is still lacking. This work is anattempt of MD simulation to establish a direct connection
between the jump diﬀusion DJump and the violation of the SE
relation. Here, we directly evaluate DJump of the watermolecules and subsequently separate the same from the total
diﬀusion D of the water molecules. The above decomposition
allows us to check the coupling of the viscosity with diﬀerent
Received: February 22, 2019
Revised: July 22, 2019
Published: July 30, 2019
A

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B
components of D and thereby gain further insights into therole of the translational jump diﬀusion of the water molecules
for the increasing breakdown of the SE relation with decreasing
temperature.
The organization for the remainder of this paper is asfollows. The method and model, including the simulation
details, are detailed in Section 2. Section 3 details thesimulation results and discussions. Concluding remarks are
oﬀered in Section 4.

experimental and simulated dynamical quantities. Figure S1(Section S1) of the Supporting Information presents thecomparison between the present simulated densities withearlier experimental47−49 and simulated densities.50,51 The
comparison of the present simulated D and η values with thosereported in an earlier simulation study53 is shown in Figure S2(Section S1) of the Supporting Information. The comparison
between the experimental4,52 and the current simulated
diﬀusion and viscosity is presented in Section 3.

2. METHOD AND MODEL2.1. Simulation Details. We have used the GROMACS
package41 for performing the MD simulations. We haveconsidered 2000 water molecules, which are placed in a cubic
simulation box with periodic boundary conditions employed inall sides. The water molecules are modeled by the TIP4P/2005
force ﬁeld,42 which is one of the most successful water models
in predicting various thermodynamic, structural, and dynamical
properties at a wide range of temperatures including the deeply
supercooled regime.
We have simulated eight diﬀerent temperatures: 210, 220,
230, 240, 250, 260, 280, and 300 K at 1 bar pressure. The ﬁrst
ﬁve temperatures represent the supercooled state of water
because TIP4P/2005 water freezes at 252 K.42 Each of thesesystems is equilibrated for 50 ns time, during which thetemperature is kept constant at the respective desired
temperatures using the Nosé−Hoover thermostat43,44 andthe pressure at 1 bar using the Parrinello−Rahman barostat.45
We have chosen 0.5 ps times for both the pressure andtemperature coupling constant. We extend each of thesimulations for a suﬃciently long time for obtaining the
production trajectory for analyses, which aggregates to a totalof 2.5 μs simulation time. The production trajectory durations
ttraj for the simulations at diﬀerent temperatures are detailed inTable 1.

3. RESULTS AND DISCUSSION
3.1. Diﬀusion−Viscosity Decoupling. We calculate theself-diﬀusion coeﬃcient D of the water molecules using themean-square displacement (MSD) routelim ⟨|Δr(t )|2 = 6Dt ⟩. Here, the MSD ⟨|Δr(t)|2⟩ is calculatedt →∞using the equation ⟨|Δr(t)|2⟩ = ⟨|rk(t) − rk(0)|2⟩, where rk(t)and rk(0) are the positions of the kth water molecule at time t
and time t = 0, respectively. Figure 1 presents MSD against

Figure 1. MSDs of water molecules as a function of time at diﬀerent
temperatures (210−300 K).

time at all eight diﬀerent temperatures. The characteristic
MSD plots at lower temperatures look similar to those forother supercooled liquids. The intermediate plateau region,
occurring between the short-time ballistic part and the longtime diﬀusive part, is subdiﬀusive with fractional timedependence ⟨|Δr(t)|2⟩ ∝ tα;0 < α < 1.54 This subdiﬀusive
part of MSD appears from the rattling of a molecule inside thesolvent cage, while the former is momentarily trapped inside
the cage.
Figure 2a presents the Arrhenius-type plot for simulated Dvalues, calculated from the slope of the MSD plots at a longtime, as a function of temperature. The simulated D values are
in close agreement with the available experimental values.52
We have also compared in Figure S2 (Section S1) of the
Supporting Information the present simulated D values withearlier simulated values.53 Both the simulated and theexperimental values are also listed in Table 2. (The error bar
analysis for the diﬀusion coeﬃcient has been detailed in
Section S2 of the Supporting Information.) We found thatabove T = 230 K temperature, the simulated data ﬁt quite wellwith the empirical Vogel−Fulcher−Tamman (VFT)-type
relationship, D = D0 exp[−B/(T − T0)], where T is thetemperature, T0 is generally close to the glass-transition
temperature, and D0 and B are other ﬁtting constants. We
obtain the following ﬁtting parameters after regressing the VFT
equation onto the simulated D: D0 = 3.69 × 10−4 cm2/s, B =364.05 K, and T0 = 170.9 K. These values are consistent withthe experimental data.52

Table 1. Production Trajectory Duration ttraj for
Simulations at Diﬀerent Temperatures
T (K)production trajectory length ttraj (ns)

210
220
230
240
250
260
280
300

500
500
400
400
300
200
100
100

We have used the Leapfrog−Verlet algorithm for solving theequations of motions every 2 fs time. The production
trajectory is saved at a regular interval of 100 fs for T < 280
K and 20 fs otherwise. The higher time resolution of the
trajectories at a higher temperature is required for correct
estimation of the viscosity coeﬃcient. A cutoﬀ distance aroundhalf of the box length is applied for Lennard-Jones interactions.
Particle mesh Ewald summation technique and LINCS
algorithm46 are used to handle the long-range Coulombinteractions and to constrain the water O−H bonds,
respectively.
The validity of the above simulation protocol is ensured bygood agreement of various simulated quantities (such as
densities, diﬀusion, viscosity, and so forth) with earlier
B

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. Comparison between the simulated and the experimentally measured transport coeﬃcients.4,52 Arrhenius plots of simulated andmeasured translational diﬀusion coeﬃcient D (a) and shear viscosity coeﬃcient η (b). The simulated data are ﬁtted with the VFT equation (bluedashed line) in both (a,b). (c) Temperature-dependent simulated and experimental4 Dη/T values. Both the theoretical and the experimental values
are normalized to 300 K temperature.

Table 2. Comparison between the Simulated and the Experimental Parameters: Diﬀusion Coeﬃcient (D), Viscosity Coeﬃcient(η), and Dη/Ta
T (K)D/10−5 (D∞/10−5) (cm2/s) simulation

210
220
230
240
250
260
280
300

0.0025 (0.0025)
0.0177 (0.0184)
0.0773 (0.0813)
0.1889 (0.2010)
0.3764 (0.4020)
0.6140 (0.6601)
1.2945 (1.4009)
2.2173 (2.4200)

D/10−5 (cm2/s)
experimentη (cp)
simulation

0.204
0.385
0.650
1.440
2.300

1664
155.1
30.56
10.52
5.194
2.990
1.396
0.793

η (cp)
experiment

Dη/T/10−10 (D∞η/T)/10−10(g cm s−2 K−1) (simulation)

Dη/T/10−10 (g cm s−2 K−1)
(experiment)

12.68
5.392
3.059
1.425
0.892

19.81 (19.81)
12.48 (13.00)
10.27 (10.77)
8.280 (8.81)
7.820 (8.35)
7.060 (7.59)
6.450 (6.980)
5.861 (6.390)

10.80
8.300
7.650
7.330
6.840

a
D∞ and D∞η/T values in the parentheseslisted in the second and sixth columnsare the corrected values of D and Dη/T, respectively, forﬁnite size eﬀect.

Now, we calculate the viscosity coeﬃcient η of water at
diﬀerent temperatures using the Green−Kubo relation via
stress tensor correlation function. Figure 2b presents an
Arrhenius plot of η as a function of temperature. Section S2 of
the Supporting Information details the error bar analysis forthe viscosity coeﬃcient. The simulated η values match very
well with both the experimentally measured4 and previously
simulated25,26 η values at all the temperatures. See Figure S2bof the Supporting Information, where the present simulated ηvalues are compared with earlier simulation.25,26 Table 2
summarizes both the simulated and the experimental η values.4
The visible agreement between the presented simulation andexperimental η values validates the simulation force ﬁeld andthe protocols adopted in the present work. Above 230 K
temperature, the simulated data ﬁt quite well with theempirical VFT equation, η = η0 exp[B/(T − T0)]; η0 = 4.21
× 10−2 cp, B = 378.4 K, and T0 = 171.4 K. However, the same
equation does not ﬁt very well with the simulated η values
below 230 K temperature. However, we do not ﬁt the
Arrhenius equation only with the two temperatures: 210 and
220 K. Therefore, we cannot conﬁrm the fragile-to-strong
liquid transition55−59 below the 230 K temperature in thepresent study. It is worthwhile mentioning that in a recent
study by Saito et al.,60 the fragile-to-strong liquid transition has
been observed for TIP4P/2005 water at 190 K temperature.
Therefore, we do not expect the transition to be visible above
210 K temperature, the lowest temperature in our study.
We now check the validity of the SE relation using thepresent D and η values. The simulated and experimental Dη/T
values for the whole temperature range are listed in Table 2.
Figure 2c exhibits the normalized simulated and experimental
Dη/T values as a function of temperature, which should beconstant if the SE relation holds strictly. Note that theexperimental Dη/T values are obtained from the measured D52
and η values.4 Normalization of the experimental and thesimulated Dη/T values are done for the value at T = 300 K. It
should be noted that the SE equation already breaks down at
300 K temperature. However, normalization of the Dη/T
values to 300 K temperature does not result in a serious error
in the calculation because the breakdown of the SE equation is
not much intense at 300 K temperature.2,4 Figure 2c clearlyshows the gradual deviation of the normalized simulated Dη/Tfrom unity as the temperature decreases from the room
temperature. This indicates an increasing violation of the SE
equationwhich reaches ∼80%on cooling down the system
to 210 K. This is consistent with the available experimental
result.4 Another parameter for quantifying the agreement
between the simulated and experimentally measured Dη/T
values is the phenomenological exponent ω, which can becalculated by ﬁtting with the equation, D ∝ η−1+ω.9,61−64 We
obtain the value of ω = 0.12. For hard spheres and othersimilar model liquids, ω = 0.2261 in the SE equation
breakdown regime. For conﬁned water, the value of ω =0.23, which was calculated via the equation D ∝ T−1+ω.2 The
experiment by Dehaoui et al.4 showed that the ﬁtting of thedata using D ∝ (η/T)ξ gives the value of ξ = −1 at hightemperature (T > 350 K) and ξ = −0.8 at low temperature (T< 350 K). Fitting of our simulated data with the above
equation gives ξ = −0.85. Interestingly, the study of Mallamace
et al.65 suggested a putative universal crossover of ξ from −1 to−0.85 upon cooling for nine diﬀerent glass-forming liquids.
Therefore, our simulated SE breakdown parameters are inclose agreement with the experimental ones.
Although 2000 water molecules (suﬃciently large systemsize) are considered in the present simulation, the simulated Dvalues may still suﬀer from ﬁnite size eﬀect because of theC

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. (a) Non-Gaussian parameter α2(t) as a function of time for all the temperatures studied. Times t* of α2(t) maxima are listed in Table 3.
theo(b) The self-part of the van Hove correlation functions Gsimu
S (r,t*) (solid line) and the corresponding ideal Gaussian distribution (GS (r,t*))(dashed line) at time t*, when α2(t) is maximum at T = 210 K. The similar plots for other temperatures are presented in Figure S3 of the
Supporting Information. (c) Distance traversed by one water molecule, λj(t,t*), in the jth trajectory segment. The horizontal line indicates the
cutoﬀ distance r2. The red circles, which are above the cutoﬀ distance r2, represent the jump trajectory segments.

occurrence. The most straightforward method of identifying ajump occurrence is via calculating the displacement of amolecule from its position at the beginning of the trajectory t =0.24,71−75 Despite the simplicity of the method, this approach
fails in the quantitative analysis of jump diﬀusion coeﬃcientbecause of the following reason. This method correctly
identiﬁes only those jumps, where the initial and the ﬁnal
positionsthe initial position is the position of the molecule
just before the jump occurrence and the ﬁnal position refers to
the new position of the molecule after the jump occurrence
of the jumping molecule are collinear with the position at t = 0.
The inaccuracies of identiﬁcation of the jump occurrence andestimation of the jump length are very high when the above
three coordinates deviate from linearity.24 The methodcompletely loses its hold when the three coordinates form aright angle triangle. In that case, the jump displacement shows
either no peak at all or a very small peak with an intensity too
low to detect among the thermal noise. This method,
therefore, underestimates the contribution of jump diﬀusion
to the overall diﬀusion.24 A more quantitative method is,therefore, indispensable.
A quantitative method38 was used earlier to estimate thejump diﬀusion of a model glass-forming liquid. The algorithm
is based on the understanding that the ﬂuctuation of a caged
particle’s position from the average position on a given timescale is of the order of the Debye−Waller (DW) factor.
Therefore, this method considers a particle as a jumping one,which has the ﬂuctuation higher than the DW factor. The
caged particles are those which have lower ﬂuctuations thanthe DW factor. This method successfully estimated a single
particle’s jump diﬀusion coeﬃcient, which was eventually
utilized to calculate the diﬀusion constant of the glass-forming
liquid. However, it is not clear whether the jump diﬀusion
coeﬃcientobtained from the above methodis only due tohydrodynamic interactions between periodic images of thesimulation boxes. Therefore, we correct the D values fromsystem size eﬀects. It has been found in theoretical studies thatfor an inﬁnite liquid, the self-diﬀusion coeﬃcient D∞ can beobtained using the equation54,66,67 D∞ = D + 2.837kBT/6πηL.
Here, D is the self-diﬀusion coeﬃcient without the correction(as reported in Table 2), kB is the Boltzmann constant, T is thetemperature, η is the viscosity of the medium, and L is thesimulation box length. The calculated D∞ values are listed inTable 2, where the D∞ values are only slightly diﬀerent from D.
This indicates that 2000 water molecules constitute a
suﬃciently large system size. We have also calculated thesystem size corrected D∞η/T values for diﬀerent temperatures.
The D∞η/T values are also listed in Table 2. As expected,
D∞η/T and Dη/T values are very close to each other.
Last, Figure 2 and Table 2 clearly show that neither theexperiment nor the present simulation shows abrupt changes in
diﬀusion coeﬃcient D, viscosity coeﬃcient η, and Dη/T as thetemperature crosses the simulated Widom line (obtained for
TIP4P/2005 water)68 at T ≈ 230 K for 1 bar pressure. This isconsistent with an earlier simulation study for the E3B3 watermodel.69,70 However, simulation studies at higher pressuressuggested that the diﬀusion coeﬃcient shows an abrupt change
near the Widom line only if the pressure is closer to thetheoretical LLCP.69,70 Therefore, the key reason for notpredicting the rapid changes in diﬀusion coeﬃcient D, viscosity
coeﬃcient η, and Dη/T by the present simulation is mainly theconsideration of 1 bar pressure, which is much less than thesimulated critical pressure for TIP4P/2005 water (T = 193 K,
p = 1350 bar).68
3.2. Analysis of Jump Diﬀusion. Now, we turn our focus
on the quantitative jump analysis for estimating the jumpdiﬀusion coeﬃcient DJump of the water molecules at all thesimulated temperatures. One of the most crucial steps of theanalysis is the correct identiﬁcation of a translational jump
D

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 3. Temperature-Dependent Time Scale t*at Which
the Non-Gaussian Parameter α2(t) Is the Maximumandthe Crossing Distances between Simulated and Theoretical
Self-Part of the van Hove Correlation Functions

large amplitude displacements, which are responsible for the
breakdown of the SE equation in the supercooled liquid.
An alternative method,13 based on oscillator model theory(ﬁrst applied in understanding particle diﬀusion through acylindrical nanopore),76 was used to calculate the jumpdiﬀusion coeﬃcient (neck diﬀusion) of a small solute particle
in a viscous liquid. This method successfully captured theimportance of jump diﬀusion of particles in a glass-forming
liquid. However, it is not clear whether these methods areapplicable for estimating the jump diﬀusion coeﬃcient of abigger solute or solvent in a more complex environment.
We have used here a diﬀerent method, which is based on thetechnique, previously developed by Raptis et al.77,78 and later
modiﬁed by Araque et al.79 This method basically examined
the radii of gyration Rg(t,Δt)of diﬀerent segments of themolecular trajectories in a three-dimensional positional
coordinate space. Rg(t,Δt) for the particular trajectory segmentof length Δt (or n number of time steps) is calculated using thefollowing equation.

T (K)t* (ps)r1 (nm)r2 (nm)

210
220
230
240
250
260
280
300

560
88.6
22.0
9.00
4.30
3.00
1.60
1.10

0.086
0.092
0.098
0.104
0.105
0.111
0.115
0.120

0.215
0.230
0.225
0.240
0.245
0.260
0.270
0.300

water molecule in the jth trajectory segment of t* length isλj(t,t*), which can be calculated using the following
formula77−79
λj(t , t *) = 2R g(t , t *)(4)n

R g(t , Δt ) =1∑ [ri(t ; Δt ) − rCM(t ; Δt )]2
n i=1

We note that the translational jump occurrences are notubiquitous in all these trajectory segments. At this point, we
employ an eﬃcient method to identify the translational jump
segments correctly. The method79 uses the self-part of the van
Hove correlation function Gsimu
S (r,Δt), which is calculatedusing the following equation80,81(1)

In eq 1, ri(t;Δt) and rCM(t;Δt) are, respectively, the positionof one molecule at the ith time frame and the center of mass ofthe trajectory segment of length Δt. rCM(t;Δt) is calculated by
the following equation.

GSsimu(r , Δt ) =nrCM(t , Δt ) =1∑ ri(t ; Δt )
n i=1

(2)

k=1

t0

In eq 5, rk(t0) is the position of the kth molecule at t0 time,rk(t0 + Δt) is the position of the same particle at t0 + Δt time,and Nw is the total number of water molecules. Gsimu
S (r,Δt) isaveraged over both Nw and all possible initial t0 times. The
theodeviation of Gsimu
S (r,Δt) from the Gaussianity (GS (r,Δt) =
80,81
2
3/2
2
2
[3/2π⟨r (Δt)⟩] exp[−3r /2⟨r (Δt)⟩])is maximum attheotime Δt = t*. Both the Gsimu
S (r,t*) and GS (r,t*) are plotted in
Figure 3b for T = 210 K and in Figure S3 of the Supporting
Information for all the temperatures. GSsimu(r,t*) crosses
Gtheo
S (r,t*) at two characteristic r values, r1 and r2, by which
we deﬁne the cage and the jump trajectories, respectively. At
the smaller r limit (r < r1), the actual displacements of thewater molecules are less than the expected value obtained from
Gtheo
S (r,t*). Therefore, the jth trajectory segment of length
λj(t,t*) < r1 is a cage trajectory. On the other hand, at a larger r
limit (r > r2), the actual displacement of the water molecule islarger than the displacement predicted by Gtheo
S (r,t*). Table 3
presents the numerical values of r1 and r2 for the diﬀerent
temperatures. Therefore, the trajectory segments, where thedistance traversed by the molecule λj(t,t*) > r2, is categorized
as the jump trajectory.79
Before moving forward, we make an important point. For
characterizing a jump trajectory segment, we have adopted theabove method, which is also detailed in ref 79. The key
diﬀerence between the present method and the method by
Raptis et al.77,78 is as follows. In both the methods, λj(t,Δt) is
calculated by eq 4. However, there is one issue regarding thecorrect identiﬁcation of the peaks (out of all the peaks of
λj(t,Δt) plotted as a function of t), which are responsible forthe actual jump. In order to resolve the issue, the method by
Raptis et al. used a weighted jump length. Therefore, slightly
big jumps contribute less to the total diﬀusion. This method,(3)

where ⟨r (t)⟩ is the MSD of the kth molecule ⟨r (t)⟩ = ⟨|rk(t)− rk(0)|2⟩ and ⟨r4(t)⟩ is written as ⟨r4(t)⟩ = ⟨|rk(t) − rk(0)|4⟩.
However, we have also checked the dependence of Δt value on
the result and found that a change of Δt value by a maximum
30% does not appreciably change the results. The above
analysis has been described in Section S3 of the Supporting
Information. Also, the peak of α2(t) corresponds to the highest
heterogeneous dynamics of the molecules at time t*.80,81
Figure 3a presents α2(t) against time in the current range oftemperatures (T = 210−300 K). The peak intensity of α2(t)increases with decreasing temperature. This indicates anincrease of dynamical heterogeneity because of the decreasein temperature. Also, the t* value, listed in Table 3, increases
from ∼1 to ∼560 ps on decreasing the temperature from T =300 to 210 K. Therefore, the required length Δt of the
trajectory segmentthe input of eqs 1 and 2increases withdecreasing temperature.
We divide the trajectories of each water molecule into
multiple segments each of length t*. We calculate Rg(t,t*) for
all the trajectory segments. The distance traversed by one2Nw∑ δ(r − |rk(t0) − rk(t0 + Δt )|)(5)

As the diﬀusion of the water molecules increases with thetemperature, consideration of the same Δt for all thetemperatures can lead to unreasonable results. For example,
Δt = 10 ps trajectory segment at 210 K temperature spreads inmuch smaller space than that at 300 K. In order to avoid thisapparent inconsistency, we choose diﬀerent Δt for diﬀerent
temperatures. It is observed that consideration of Δt as thecharacteristic time t*when the non-Gaussian parameterα2(t) is maximumworks perfectly well under any situation.58
We calculate α2(t) using the following equationα2(t ) = 3⟨r 4(t )⟩/5⟨r 2(t )2 ⟩ − 1

1
Nw

2

E

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 4. Four representative jump and cage trajectories of length t* for one water molecule at two representative temperatures: 210 K (top panels)
and 250 K (bottom panels). The color-mapping denotes the direction of the time evolution of a trajectory.

therefore, considers all the peaks of the weighted λj(t,Δt) as“jumps”. The deﬁnition of jump translation is therefore slightly
confusing. On the other hand, the present method79 uses eq 4
only (without using the weighted parameter) to calculate
λj(t,Δt) and distinguish a trajectory segment as the “jump”
only when it meets the criteria, λj(t,Δt) > r2, where r2 is thetheosecond crossing point between Gsimu
S (r,Δt) and GS (r,Δt).
simu
theo
The fact that GS (r,Δt) deviates GS (r,Δt) at the most for
Δt = t* (t* is the time when α2(t) is maximum) makes t* theabsolute choice for Δt. Clearly, the present method is inaccordance with the classical picture of the jump translation,
which is believed to have emerged from the spatio-temporal
heterogeneity of the medium, the time scale of which can bequantiﬁed by t* and the deviation of Gsimu
S (r,Δt) from
80,81
24
Gtheo
(r,Δt).
Ina
recentarticle,twooftheauthors haveS
shown that the jump diﬀusion of a solute in the supercooled
water emanates from dynamic heterogeneity of the medium via
correlated translational/rotational solvent dynamics. Interestingly, the approximate jump analysis in ref 24, via visual
identiﬁcation of large sudden jumps, provided similar results tothe present work. This validates the method of identiﬁcation ofa jump occurrence in the present work.
Figure 3c exhibits an example of λj(t,t*) as a function of time
for one water molecule at T = 210 K temperature as arepresentative temperature. A peak, whose intensity crosses the
distance r2, represents a translational jump trajectory.
Conversely, we characterize a trajectory segment as the cage
trajectory, which has λj(t,t*) value less than the cutoﬀ r1. Here,the cage trajectory refers to the rattling motion of a molecule
inside the solvent cage plus the overall translation of the wholesolvent cage with the tagged molecule in the center. Figures 4
and S4 of the Supporting Information present several examplesof the cage and jump trajectory segments. It is evident from therepresentative examples that the above protocol for identifying
the jump and cage trajectories is working as expected because
the identiﬁed jump trajectory segments consist of suddenchange(s) of the position of the molecule in between initial
cage to ﬁnal cage and identiﬁed cage trajectory segments are
actually rattling trajectories.
As we have correctly identiﬁed the jump trajectory segments,
we now calculate the jump diﬀusion coeﬃcient DJump using thefollowing equation77−79

DJump =

1
νJump λ Jump 26(6)

νJump is the frequency of the translational jump occurrence of
one water molecule. νJump is calculated using the following
equationnJump
νJump =
n w t traj
(7)
λJump2 is the average square jump length, which is obtained byaveraging λj2(t,t*) over the number of jump occurrences nJump
in the entire system using the following equation2

λ Jump =lim

1nJump →∞ nJump

nJump∑ λj 2(t , t *)|for λ (t ,t*) > r
j

j=1

2

(8)

Note that eq 6 can be obtained from the random walk
problem. Because a jumping water molecule can travel a
distance λj in any random directions, and these jumps are
statistically independent of each other, the current problem
eﬀectively turns into the random walk problem. The derivationof eq 6 has been presented in Section S4 of the Supporting
Information.
We have listed, in Table 4, the numerical values of nJump,
νJump, λJump2, and the percentage of jump trajectory segments,
nJump
PJump = n t / t * × 100 (the values of ttraj are presented in
w trajTable 1) for diﬀerent temperatures. Using the above numericalvalues in eq 6, we calculate DJump of the water molecules at all
the temperatures. The values are listed in Table 5.
Figure 5a exhibits the percentage contribution of DJump to
the overall diﬀusion D of the water molecules 100(DJump/D)%
as a function of temperature. The contribution increases withdecreasing temperature and reaches approximately 50% at T =210 K. On the other hand, the percentage of jump trajectory
segments PJump increases from ∼0.12% of the total trajectory
segments at T = 300 K to slightly more than 5% at T = 210 K.
Therefore, approximately 1 jump in 20 trajectory segments
contributes 50% of the overall diﬀusion of the water molecules
at 210 K temperature. In other words, a small fraction of thejump trajectory contributes to a large fraction of the overall
diﬀusion at T = 210 K.
F

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B
Table 4. Average Number of Jump (nJump), Jump Frequency(νJump), Average Jump Displacement (λJump2), and thePercentage of the Jump Trajectory Segment (PJump) at
Diﬀerent Temperatures

T (K)

nJump

jump
frequencyνJump (ns−1)

210
220
230
240
250
260
280
300

93 000487 894
1 418 488
1 843 861
1 631 136
1 040 193
529 883
213 386

0.093
0.488
1.773
2.305
2.719
2.600
2.649
1.067

λJump2(Å2)

percentage of jump
trajectory segment
nJump
PJump = n t / t * × 100

8.12
8.01
7.13
7.62
7.67
8.47
8.94
10.6

5.21
4.32
3.90
2.07
1.17
0.78
0.42
0.12

w trajFigure 5. (a) Percentage contribution of the jump diﬀusion
coeﬃcient DJump to the total diﬀusion coeﬃcient D of the watermolecules as a function of temperature. Temperature-dependent
coupling of viscosity with (b) translational jump diﬀusion coeﬃcient
DJump, (c) overall translational diﬀusion coeﬃcient D and the residual
diﬀusion coeﬃcient DRes (D − DJump), and (d) system size eﬀect
corrected diﬀusion value D∞ and DRes∞ (D∞ − DJump∞). The error barsof Dη/T in (c,d) are the standard errors. (See Section S2 of the
Supporting Information for details of error analysis.)

We now check the validity of the SE equation using the
calculated DJump and the simulated η of the medium. Table 5
lists DJumpη/T at diﬀerent temperatures. Figure 5b displays thenormalized DJumpη/T values as a function of time where we seethat DJumpη/T increases rapidly with decreasing temperature.
The value of the latter increases by ∼200 times as thetemperature is decreased from T = 300 to 210 K. Therefore, as
expected, DJump is completely decoupled from the viscosity
almost at all the temperatures. This is expected because the SE
equation was not developed for large displacement but only forthe small step diﬀusive process. It is evident from Figure 5 thatsimilar to D and η, abrupt change of DJump is also not observed
as the temperature crosses ∼230 K temperature (the positionof the Widom line for TIP4P/2005 water at 1 bar). This isagain due to the lower simulated pressure (1 bar) condition in
the present case than the simulated critical pressure for
TIP4P/2005 water (T = 193 K, p = 1350 bar).68
Note that DJump has been calculated by considering only thelarge amplitude displacements (λj(t,t*) > r2) of the watermolecules. In other words, DJump is that component of the total
diﬀusion D, which considers only some long step displacements. Therefore, we can write the total diﬀusion coeﬃcient D
as a sum of two components. The ﬁrst component DJump comesfrom the jump translation exclusively. On the other hand, thesmaller step displacements collectively can give rise to anotherdiﬀusion coeﬃcient. We call it residual diﬀusion coeﬃcientDRes, which can be calculated from the following equation

DRes = D − DJump(9)

Note that similar decompositions of the total diﬀusion
coeﬃcient are ubiquitous in the literature, although in diﬀerent
aspects. The total diﬀusion coeﬃcient of a small solute particle
in the molecular glass-forming liquid was written as a sum ofneck diﬀusion (analogous to DJump) and the diﬀusion
coeﬃcient emanating from the SE equation.13 Therefore, thephysical meaning of DRes relies on the correct identiﬁcation of
the jump trajectory segments because the rest of the trajectory
segments contribute to DRes. Therefore, DRes is that componentof D, which completely excludes molecular jump translation
and which emanates from the small step diﬀusion process as instandard Einstein’s Brownian motion. This immediately
suggests that DRes should be coupled to the viscosity of themedium and thereby follows the SE equation.
We have calculated DRes for all the temperatures and listed
their numerical values in Table 5. We now check the validity of
the SE equation by checking the eﬀect of the temperature on
the numerical value of DResη/T, which are presented in Table

Table 5. Diﬀerent Components of Diﬀusion Coeﬃcients: Diﬀusion Coeﬃcient Due to Jump (DJump) and Residual Diﬀusion(DRes) of the Water Molecules, the Percentage Contribution of DJump to Total Diﬀusion Coeﬃcient D, DJumpη/T Values, andthe DResη/T at Eight Diﬀerent Temperaturesa
T (K)DJump/10−5(DJump∞/10−5)(cm2/s)DRes/10−5 (DRes∞/10−5)(cm2/s)100(DJump/D)%

DJumpη/T/10−10 (DJump∞η/T/10−10)(g cm s−2 K−1)

DResη/T/10−10 (DRes∞η/T/10−10)(g cm s−2 K−1)

210
220
230
240
250
260
280
300

0.0013 (0.0013)
0.0065 (0.0068)
0.0211 (0.0222)
0.0293 (0.0312)
0.0348 (0.0372)
0.0368 (0.0396)
0.0394 (0.0426)
0.0189 (0.0206)

0.0012 (0.0012)
0.0112 (0.0116)
0.0562 (0.0591)
0.1596 (0.1698)
0.3416 (0.3648)
0.5772 (0.6205)
1.2551 (1.3583)
2.1984 (2.3994)

52.00
36.72
27.30
15.51
9.245
5.993
3.044
0.852

10.301 (10.301)
4.5825 (4.7635)
2.8035 (2.9485)
1.2843 (1.3666)
0.7230 (0.7722)
0.4232 (0.4550)
0.1964 (0.2126)
0.0500 (0.0545)

9.5086 (9.5086)
7.8960 (8.2085)
7.4673 (7.8538)
6.9958 (7.4439)
7.0971 (7.5798)
6.6378 (7.1362)
6.2576 (6.7719)
5.8111 (6.3423)

The values in parenthesis are the ﬁnite size eﬀect corrected values of the jump diﬀusion coeﬃcients DJump∞ (second column), the residual diﬀusion
coeﬃcients DRes∞ (third column), DJump∞η/T (ﬁfth column), and DRes∞η/T (sixth column).
aG

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 6. Eﬀect of the Jump (DJump), Residual Diﬀusion (DRes), and DResη/T on the 10% Change of r2 Second Crossing Distanceatheo
Between Gsimu
S (r,Δt) and GS (r,Δt) at Eight Diﬀerent Temperatures
T (K)r2 (Å)

νJump (ns−1)

λJump2 (Å2)

DJump/10−5 cm2/s

DRes/10−5 cm2/s

DResη/T/10−10 cm2/s

210

r2 − 10%r2
r2 + 10%r2 − 10%r2
r2 + 10%r2 − 10%r2
r2 + 10%r2 − 10%r2
r2 + 10%r2 − 10%r2
r2 + 10%r2 − 10%r2
r2 + 10%r2 − 10%r2
r2 + 10%r2 − 10%r2
r2 + 10%

0.132
0.093
0.066
0.775
0.488
0.307
3.116
1.773
1.004
4.631
2.305
1.140
6.064
2.719
1.216
6.195
2.600
1.088
6.912
2.649
0.999
3.342
1.067
0.329

6.95
8.12
9.36
6.81
8.01
9.36
6.00
7.13
8.41
6.40
7.62
9.00
6.40
7.67
9.12
7.08
8.47
10.0
7.40
8.94
10.6
8.82
10.6
12.7

0.00150.00130.00100.00880.00650.00480.03120.02110.01400.04930.02930.0171
0.06500.03480.01850.07310.03680.01830.08540.03940.01760.04920.01890.00690.00100.00120.00150.00890.01120.01290.0461
0.0562
0.06330.1396
0.15960.17180.31140.3416
0.35790.54090.57720.59571.20911.25511.27692.1681
2.19842.21047.92389.5086
11.8866.27457.89609.0945
6.12537.46738.4106
6.1191
6.99587.5306
6.4696
7.09717.43576.2204
6.63786.85066.0282
6.25766.36635.73105.8111
5.8428

220

230

240

250

260

280

300

a

We also compared the jump frequency and average square jump length for all eight temperatures.

5. Figure 5c plots the normalized DResη/T as a function oftemperature. Note that the normalization has been done withrespect to the value at T = 300 K temperature. Veryinterestingly, unlike the Dη/T, which continuously increases
with decreasing temperature, the DResη/T value remains almost
constant.1−4 This indicates that once the translational jump
only diﬀusion is separated out from the overall diﬀusion, theresidual diﬀusion perfectly couples with the viscosity and thus
follows the SE equation. Therefore, the jump diﬀusion of themolecules is the fundamental origin for the observed
decoupling of the molecular diﬀusion from the viscosity ofthe medium. This is one of the key results of this work as it
categorically proves the concept that the origin of the wellknown diﬀusion−viscosity decoupling in supercooled water(liquid) is the translational jump diﬀusion of molecules. Twoof us have previously shown the mechanism of thesetranslational jumps in great details and the crucial role of thesynchronization between the translational and the rotationalmotion of the solvent water molecules for inducing these jump
events.24,75
Similar to the correction of the system size eﬀect for D
values, listed in Table 5, we also correct the DJump and DRes
values for inﬁnite liquid. For this, we have scaled the values of
DJump and DRes by the same correction factor D∞/D. The
corrected DJump∞ and DRes∞ are listed in Table 5. Figure 5d
plots the normalized DRes∞η/T as a function of temperature,
where the normalization is done for the D∞η/T value at T =300 K temperature. Similar to the normalized DResη/T value,
the normalized DRes∞η/T remains almost constant except aslight increase at 210 K, which may come from larger error atthat temperature.1−4

Last, before concluding, we discuss a couple of important
points. First, we discuss how the key results of the present
workmainly the DJump valuesare aﬀected by various input
parameters and jump criteria. As we consider a trajectory
segment as a jump one if λj(t,t*) > r2. Therefore, DJump must
depend on the r2 value. However, the r2 value is not an
independent parameter but the second crossing distance
between GSsimu(r,Δt) and GStheo(r,Δt), which are actuallydependent on the Δt value. Therefore, the only input
parameter in the jump analysis is the numerical value of Δt.
Once we choose the Δt value, the criterion for the jump
motion is already set. In the present study, we have consideredΔt value as the time (t*) at which the α2(t) is maximum. The
above consideration of t* as Δt is as per the previous work by
Araque et al (ref 79 of the MS). This can also be justiﬁed fromthe fact that t* is the time when the displacement distributionof molecules becomes substantially non-Gaussian in nature.
Therefore, at time t*, the dynamics of the molecules areheterogeneous at the most. However, one can check the eﬀect
of Δt value on the results, particularly DJump. We have varied
the Δt value ranging from 0.7t* (30% less than t*) to 1.3t*(30% more than t*) in a regular interval of 0.1t*(10% of t*)
and performed the complete jump analysis for each of them.
The results are presented in Figure S7 of the Supporting
Information. We see that the change is nominal if the Δt value
is increased/decreased by 30% from the t* value. Therefore,
the present method is not very much dependent on the choiceof the Δt value until it is not changed by a very large extent.
Now, we discuss the eﬀect of the r2 value on the key resultsof the present work. Evidently, DJump must depend on the r2
value. Note that the r2 value is not an independent parameterbut actually determined by the second crossing distancetheo
between Gsimu
S (r,Δt) and GS (r,Δt), which ﬁnally depend on
H

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B
the Δt value. However, one may ask the eﬀect of the results on
the subtle deviation of r2 value, which may arise due to small
error of the Gsimu
S (r,Δt). This is particularly important for thehigher temperatures, when the second crossing point between
theo
Gsimu
S (r,Δt) and GS (r,Δt) is not sharp. Therefore, we checkthe dependence of DResη/T on a slight variation of the r2 value.
We calculated the DResη/T with the r2 values changed 10%
from the original data, reported in Table 3. We havesummarized the results in Table 6. The results suggest thatthe key result DResη/T is not changed appreciably, particularly
at the higher temperature, even though various quantities(such as νJump, λJump2, DJump, and so forth) are slightly aﬀected.
This is because the contribution of DJump to the total diﬀusion
is much less signiﬁcant at high temperature.
Finally, we discuss a very crucial point. Even though eq 6 canbe derived from the random walk problemdescribed in
Section S4 of the Supporting Informationone can stillquestion the validity of the equation in the present system
where the solvent is spatiotemporally heterogeneous at thesupercooled region. Also, eq 9 has been written on thejustiﬁcation that while the DJump emanates from the large
amplitude jump motion of the water molecules, the smaller
step displacements collectively can give rise to anotherdiﬀusion coeﬃcient, DRes. Therefore, the total diﬀusion
coeﬃcient is decomposed into two components: DJump andDRes. We validate the above point by evaluating DJump and DRes
independently using the MSD route instead of using eqs 6−9.
For this, we separate the full trajectory into two parts: thetrajectory portion containing only the jump events and theleftover portion of the trajectory containing only smaller step
displacements. The conditions, which set the criteria for thejump and the smaller step trajectory segments, are λj(t,t*) > r2
and λj(t,t*) < r2, respectively. Now, because the number ofjump trajectory segments is much less than the residual
trajectory segments, visible from Table 4, presumably most of
the times two successive jump trajectory segments are
separated by one or more residual trajectory segment(s).
This introduces discontinuity of the jump trajectory, which
gives erroneous DJump value. To solve this issue, we haveconnected the last point of one jump trajectory with the ﬁrst
point of the next jump trajectory segment such that they areseparated by time. Therefore, during the extraction of the jump
trajectory from the original trajectory, the residual trajectory
segments are replaced by the same phase space variables of thelast point of the previous jump trajectory segment. Thus,
during the residual trajectory segments, the particle is at rest.
By this way, we take the jump frequency into account. We
extract the pure residual trajectory in the same way as above,
but this time during the jump trajectory segments the particle
is at rest. Now, we calculate the MSDs from the above two
trajectories (jump and the residual) separately using thefollowing equations

in between two jump trajectory segments, the displacement
does not add up. The same is true for calculating the MSDRes
using eq 11.
The jump, the residual, and the total MSD plots arepresented in Figure S8 for three representative temperatures.
Clearly, both the jump and the residual MSD components arelinearly proportional to the time. To get the DJump and DRes
values, we calculate the slopes of the jump and the residual
MSDs at a long time and compare with the values obtained via
eqs 6−9 for all the temperatures in Table S2 of the Supporting
Information. The values obtained from the MSD route are invery close agreement with the values calculated using eqs 6−9.
With the above DJump and DRes values, we have plotted in
Figure S9 of the Supporting Information the percentage
contribution of the DJump to the overall diﬀusion D of the watermolecules 100(DJump/D)%, DJumpη/T, and DResη/T at diﬀerent
temperatures and compared them with the data (obtained
from eqs 6−9) already presented in Figure 5a−c. Clearly, theresults are in close agreement with each other. This further
emboldens the validity of eqs 6−9 and the physical signiﬁcanceof the DRes, which emanates from the collective small step
displacements of the water molecules. An even closer
agreement of the results at lower temperatures may be foundby analyzing a longer production trajectory length.
For calculating DRes, instead of removing the jump trajectory
segments from the full trajectory, we can also replace thesejump trajectory segments by the small step diﬀusion channel(residual diﬀusion) during the time taken by the jumps. Thiscan give a stronger proof that DJump is responsible for theviolation of the SE relation. If the jumps last a fraction f of thetotal time interval ttraj, the residual trajectories last a total time
(1 − f)ttraj, contributing to the MSD for an amount MSDRes =6DRest. If the jumps are replaced by residual trajectory
segments, the new MSD would be MSDRes(Modiﬁed) =6DRes(eff) t = 6DRes/(1 − f)t, which is equivalent to an eﬀective
residual diﬀusion coeﬃcient DRes(eff) = DRes/(1 − f). f can becalculated from PJump by the following equation f = PJump/100.
PJump values for diﬀerent temperatures are listed in Table 4. We
have plotted, in Figure S10 of the Supporting Information,
DRes(eff)η/T as a function of time. Clearly, the diﬀerence
between DRes(eff)η/T and DResη/T is very small. This is because
f is very small at all temperatures, and the eﬀect is maximum(∼5%) at 210 K.

4. CONCLUSIONSIn conclusion, we have presented an MD simulation analysis
for quantitative estimation of the role of translational jumpdiﬀusion on the increasing decoupling of viscosity from thetranslational diﬀusion with decreasing temperature. By carefulconsideration of the translational jump trajectories of the watermolecules, we have calculated the jump-only diﬀusion
coeﬃcient DJump. As the temperature is decreased, thecontribution of DJump to the overall diﬀusion increases. By
subtracting DJump from D, we obtain the residual diﬀusion
coeﬃcient DRes. While DJump intensely decouples from theviscosity η, DRes stays coupled fairly strongly with η at all thetemperatures studied. This is a clear proof for the growing
importance of translational jump diﬀusion of the molecules forthe observed breakdown of the SE equation with decreasing
temperature.
These new ﬁndings can help in elucidating many
experimental studies featuring molecular transport properties
in a more complex chemical and biological environment, where

MSDJump = ⟨|Δr(t )|⟩2 = ⟨|rk(t ) − rk(0)|2 ⟩|for λj(t , t *) > r2
(10)
2

2

MSDRes = ⟨|Δr(t )|⟩ = ⟨|rk(t ) − rk(0)| ⟩|for λ j(t , t *) < r2(11)Here, rk(t) and rk(0) are the positions of the kth watermolecule at time t and time t = 0, respectively. Therefore, theMSDJump is calculated only for the jump trajectory segments,
set by the criteria λj(t,t*) > r2 (for the jth trajectory segment of
the kth particle). If a residual jump trajectory segment comesI

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

supercooled water. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 12974−
12978.
(3) Kumar, P. Breakdown of the Stokes-Einstein relation insupercooled water. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 12955−
12956.
(4) Dehaoui, A.; Issenmann, B.; Caupin, F. Viscosity of deeply
supercooled water and its coupling to molecular diffusion. Proc. Natl.
Acad. Sci. U.S.A. 2015, 112, 12020−12025.
(5) Cicerone, M. T.; Ediger, M. D. Enhanced translation of probe
molecules in supercooled o-terphenyl: Signature of spatiallyheterogeneous dynamics? J. Chem. Phys. 1996, 104, 7210−7218.
(6) Ediger, M. D. Spatiallyheterogeneousdynamics Insupercooledliquids. Annu. Rev. Phys. Chem. 2000, 51, 99−128.
(7) Fujara, F.; Geil, B.; Sillescu, H.; Fleischer, G. Translational androtational diffusion in supercooled orthoterphenyl close to the glasstransition. Z. Phys. B: Condens. Matter 1992, 88, 195−204.
(8) Chang, I.; Fujara, F.; Geil, B.; Heuberger, G.; Mangel, T.;
Sillescu, H. Translational and rotational molecular motion insupercooled liquids studied by NMR and forced Rayleigh scattering.
J. Non-Cryst. Solids 1994, 172−174, 248−255.
(9) Debenedetti, P. G.; Stillinger, F. H. Supercooled liquids and theglass transition. Nature 2001, 410, 259.
(10) Banerjee, P.; Yashonath, S.; Bagchi, B. Rotation driven
translational diffusion of polyatomic ions in water: A novel
mechanism for breakdown of Stokes-Einstein relation. J. Chem.
Phys. 2017, 146, 164502.
(11) Guchhait, B.; Daschakraborty, S.; Biswas, R. Medium
decoupling of dynamics at temperatures ∼100 K above glasstransition temperature: A case study with (acetamide + lithium
bromide/nitrate) melts. J. Chem. Phys. 2012, 136, 174503.
(12) Park, Y.; Kim, J.; Sung, B. J. Translation-rotation decoupling oftracers of locally favorable structures in glass-forming liquids. J. Chem.
Phys. 2017, 147, 124503.
(13) Acharya, S.; Nandi, M. K.; Mandal, A.; Sarkar, S.;
Bhattacharyya, S. M. Diffusion of Small Solute Particles in Viscous
Liquids: Cage Diffusion, a Result of Decoupling of Solute-Solvent
Dynamics, Leads to Amplification of Solute Diffusion. J. Phys. Chem. B
2015, 119, 11169−11175.
(14) Cicerone, M. T.; Ediger, M. D. Enhanced translation of probe
molecules in supercooled o-terphenyl: Signature of spatiallyheterogeneous dynamics? J. Chem. Phys. 1996, 104, 7210.
(15) Richert, R. Heterogeneous dynamics in liquids: fluctuations inspace and time. J. Phys.: Condens. Matter 2002, 14, R703−R738.
(16) Cicerone, M. T.; Ediger, M. D. Relaxation of spatiallyheterogeneous dynamic domains in supercooled ortho-terphenyl. J.
Chem. Phys. 1995, 103, 5684.
(17) Kaufman, L. J. Heterogeneity in Single-Molecule Observablesin the Study of Supercooled Liquids. Annu. Rev. Phys. Chem. 2013, 64,
177−200.
(18) Cicerone, M. T.; Zhong, Q.; Tyagi, M. Picosecond Dynamic
Heterogeneity, Hopping, and Johari-Goldstein Relaxation in GlassForming Liquids. Phys. Rev. Lett. 2014, 113, 117801.
(19) Giovambattista, N.; Mazza, M. G.; Buldyrev, S. V.; Starr, F. W.;
Stanley, H. E. Dynamic Heterogeneities in Supercooled Water†. J.
Phys. Chem. B 2004, 108, 6655−6662.
(20) Kob, W.; Donati, C.; Plimpton, S. J.; Poole, P. H.; Glotzer, S. C.
Dynamical Heterogeneities in a Supercooled Lennard-Jones Liquid.
Phys. Rev. Lett. 1997, 79, 2827−2830.
(21) Karmakar, S.; Dasgupta, C.; Sastry, S. Growing length and timescales in glass-forming liquids. Proc. Natl. Acad. Sci. U.S.A. 2009, 106,3675−3679.
(22) Yamamoto, R.; Onuki, A. Dynamics of highly supercooled
liquids: Heterogeneity, rheology, and diffusion. Phys. Rev. A: At., Mol.,
Opt. Phys. 1998, 58, 3515−3529.
(23) Andersen, H. C. Molecular dynamics studies of heterogeneous
dynamics and dynamic crossover in supercooled atomic liquids. Proc.
Natl. Acad. Sci. U.S.A. 2005, 102, 6686−6691.
(24) Dubey, V.; Kumar, N.; Daschakraborty, S. Importance ofSolvents’ Translational-Rotational Coupling for Translational Jump ofstrong diﬀusion−viscosity decoupling is prevalent. Also, a
modiﬁed version of the above methodology can be used for
calculating the rotational jump-only diﬀusion coeﬃcient for anonassociated liquid. This would generalize the existingrotational jump model for liquid water.82,83
Last, we have observed continuous change of various
simulated quantities, total diﬀusion coeﬃcient D, viscosity
coeﬃcient η, Dη/T, jump diﬀusion coeﬃcient D Jump ,percentage contribution of the jump diﬀusion to the overall
diﬀusion 100(DJump/D)%, etc., on the decreasing temperatureof the medium without an abrupt change as the simulated
Widom line is crossed. This is because the pressure ismaintained at 1 bar in the present simulation, which is muchless than the simulated critical pressure for TIP4P/2005 water(T = 193 K, p = 1350 bar).70 The eﬀect of the Widom line
position is generally observed once the simulated condition
reaches close to the LLCP. Further studies (e.g., pressure
dependence on jump diﬀusion coeﬃcient) are required in thisfront.

■ASSOCIATED CONTENT

S Supporting Information*

The Supporting Information is available free of charge on the
ACS Publications website at DOI: 10.1021/acs.jpcb.9b01719.

■

Comparison between simulated and experimental
density, comparison between simulated and experimental dynamical properties, dependence of the results on
input parameters, comparison between simulated andtheoretical van Hove correlation functions, someexamples of jump and cage trajectories for all thetemperatures, section for error analysis, derivation of eq
1
6 of the main text (DJump = 6 νJump λ Jump 2 ) from random
walk problem, and direct evaluation of DJump and DRes
values using the MSD route (PDF)

AUTHOR INFORMATION

Corresponding Author*E-mail: snehasis@iitp.ac.in
ORCID

Snehasis Daschakraborty: 0000-0002-2694-7142
Author Contributions†

S.D. and V.D. contributed equally.

Notes

The authors declare no competing ﬁnancial interest.

■

ACKNOWLEDGMENTS
We acknowledge Prof. James T. Hynes, Prof. Jannis Samios,Prof. Claudio E. Margulis, Prof. Ranjit Biswas, and Dr. Sandipa
Indra for fruitful discussion and IIT Patna for the computational facility. S.D. and V.D. acknowledge IIT Patna for theirresearch fellowships. We acknowledge the anonymous
reviewers for their constructive criticism and continuous eﬀort
for enhancing the quality of the work.

■

REFERENCES(1) Qvist, J.; Mattea, C.; Sunde, E. P.; Halle, B. Rotational dynamics
in supercooled water from nuclear spin relaxation and molecularsimulations. J. Chem. Phys. 2012, 136, 204505.
(2) Chen, S.-H.; Mallamace, F.; Mou, C.-Y.; Broccio, M.; Corsaro,C.; Faraone, A.; Liu, L. The violation of the Stokes-Einstein relation in
J

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B
a Small Hydrophobic Solute in Supercooled Water. J. Phys. Chem. B
2018, 122, 7569−7583.
(25) Kawasaki, T.; Kim, K. Identifying time scales for violation/
preservation of Stokes-Einstein relation in supercooled water. Sci. Adv.
2017, 3, No. e1700399.
(26) Guillaud, E.; Merabia, S.; de Ligny, D.; Joly, L. Decoupling ofviscosity and relaxation processes in supercooled water: a molecular
dynamics study with the TIP4P/2005f model. Phys. Chem. Chem.
Phys. 2017, 19, 2124−2130.
(27) Galamba, N. On the hydrogen-bond network and the nonArrhenius transport properties of water. J. Phys.: Condens. Matter
2017, 29, 015101.
(28) Giovambattista, N.; Mazza, M. G.; Buldyrev, S. V.; Starr, F. W.;
Stanley, H. E. Dynamic Heterogeneities in Supercooled Water†. J.
Phys. Chem. B 2004, 108, 6655−6662.
(29) Kumar, P.; Buldyrev, S. V.; Becker, S. R.; Poole, P. H.; Starr, F.
W.; Stanley, H. E. Relation between the Widom line and thebreakdown of the Stokes-Einstein relation in supercooled water. Proc.
Natl. Acad. Sci. U.S.A. 2007, 104, 9575−9579.
(30) Furukawa, A.; Kim, K.; Saito, S.; Tanaka, H. Anisotropic
Cooperative Structural Rearrangements in Sheared Supercooled
Liquids. Phys. Rev. Lett. 2009, 102, 016001.
(31) Parmar, A. D. S.; Sengupta, S.; Satry, S. Length-Scale
Dependence of the Stokes-Einstein and Adam-Gibbs Relations inModel Glass Formers. Phys. Rev. Lett. 2017, 119, 056001.
(32) Kumar, P. Breakdown of the Stokes-Einstein relation insupercooled water. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 12955−
12956.
(33) Goldstein, M. Viscous liquids and the glass transition: apotential energy barrier picture. J. Chem. Phys. 1969, 51, 3728−3739.
(34) Sastry, S.; Debenedetti, P. G.; Stillinger, F. H. Signatures of
distinct dynamical regimes in the energy landscape of a glass-forming
liquid. Nature 1998, 393, 554−557.
(35) Mukherjee, A.; Bhattacharyya, S.; Bagchi, B. Pressure andtemperature dependence of viscosity and diffusion coefficients of a
glassy binary mixture. J. Chem. Phys. 2002, 116, 4577−4586.
(36) Bhattacharyya, S.; Bagchi, B. Anisotropic local stress andparticle hopping in a deeply supercooled liquid. Phys. Rev. Lett. 2002,
89, 025504.
(37) Cicerone, M. T.; Zhong, Q.; Tyagi, M. Picosecond dynamic
heterogeneity, hopping, and Johari-Goldstein relaxation in glassforming liquids. Phys. Rev. Lett. 2014, 113, 117801.
(38) Pastore, R.; Coniglio, A.; Pica, M. From cage-jump motion tomacroscopic diffusion in supercooled liquids. Soft Matter 2014, 10,5724−5728.
(39) Pica Ciamarra, M. P.; Pastore, R.; Coniglio, A. Particle jumps instructural glasses. Soft Matter 2016, 12, 358−366.
(40) Bhattacharyya, S. M.; Bagchi, B.; Wolynes, P. G. Facilitation,complexity growth, mode coupling, and activated dynamics insupercooled liquids. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 16077−
16082.
(41) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark,A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J.
Comput. Chem. 2005, 26, 1701−1718.
(42) Abascal, J. L. F.; Vega, C. A general purpose model for thecondensed phases of water: TIP4P/2005. J. Chem. Phys. 2005, 123,
234505.
(43) Nosé, S. A Unified Formulation of the Constant Temperature
Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519.
(44) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space
Distributions. Phys. Rev. A 1985, 31, 1695−1697.
(45) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single
Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52,
7182−7190.
(46) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M.
LINCS: A Linear Constraint Solver for Molecular Simulations. J.
Comput. Chem. 1997, 18, 1463−1472.
(47) Schufle, J. A.; Venugopalan, M. Specific volume of liquid waterto −40°C. J. Geophys. Res. 1967, 72, 3271−3275.

(48) Kell, G. S. Density, thermal expansivity, and compressibility ofliquid water from 0.deg. to 150.deg.. Correlations and tables foratmospheric pressure and saturation reviewed and expressed on 1968
temperature scale. J. Chem. Eng. Data 1975, 20, 97−105.
(49) Hare, D. E.; Sorensen, C. M. The density of supercooled water.
II. Bulk samples cooled to the homogeneous nucleation limit. J. Chem.
Phys. 1987, 87, 4840.
(50) González, M. A.; Valeriani, C.; Caupin, F.; Abascal, J. L. F. A
comprehensive scenario of the thermodynamic anomalies of waterusing the TIP4P/2005 model. J. Chem. Phys. 2016, 145, 054505.
(51) Sumi, T.; Sekino, H. Effects of hydrophobic hydration on
polymer chains immersed in supercooled water. RSC Adv. 2013, 3,
12743−12750.
(52) Price, W. S.; Ide, H.; Arata, Y. Self-Diffusion of Supercooled
Water to 238 K Using PGSE NMR Diffusion Measurements. J. Phys.
Chem. A 1999, 103, 448−450.
(53) Montero de Hijes, P.; Sanz, E.; Joly, L.; Valeriani, C.; Caupin,
F. Viscosity and self-diffusion of supercooled and stretched waterfrom molecular dynamics simulations. J. Chem. Phys. 2018, 149,
094503.
(54) Sengupta, S.; Karmakar, S. Distribution of diffusion constants
and Stokes-Einstein violation in supercooled liquids. J. Chem. Phys.
2014, 140, 224505.
(55) Shi, R.; Russo, J.; Tanaka, H. Origin of the emergent fragile-tostrong transition in supercooled water. Proc. Natl. Acad. Sci. U.S.A.
2018, 115, 9444.
(56) Angell, C. A. Water II is a “strong” liquid. J. Phys. Chem. 1993,
97, 6339.
(57) Ito, K.; Moynihan, C. T.; Angell, C. A. Thermodynamic
determination of fragility in liquids and a fragile-to-strong liquidtransition in water. Nature 1999, 398, 492.
(58) Gainaru, C.; Agapov, A. L.; Fuentes-Landete, V.; AmannWinkel, K.; Nelson, H.; Köster, K. W.; Kolesnikov, A. I.; Novikov, V.
N.; Richert, R.; Böhmer, R.; Loerting, T.; Sokolov, A. P. Anomalously
large isotope effect in the glass transition of water. Proc. Natl. Acad.
Sci. U.S.A. 2014, 111, 17402.
(59) Loerting, T.; Fuentes-Landete, V.; Handle, P. H.; Seidl, M.;
Amann-Winkel, K.; Gainaru, C.; Böhmer, R. The glass transition in
high-density amorphous ice. J. Non-Cryst. Solids 2015, 407, 423.
(60) Saito, S.; Bagchi, B.; Ohmine, I. Crucial role of fragmented andisolated defects in persistent relaxation of deeply supercooled water. J.
Chem. Phys. 2018, 149, 124504.
(61) Kumar, S. K.; Szamel, G.; Douglas, J. F. Nature of the
breakdown in the Stokes-Einstein relationship in a hard sphere fluid. J.
Chem. Phys. 2006, 124, 214501.
(62) Eaves, J. D.; Reichman, D. R. Spatial dimension and the
dynamics of supercooled liquids. Proc. Natl. Acad. Sci. U.S.A. 2009,
106, 15171.
(63) Charbonneau, B.; Charbonneau, P.; Jin, Y.; Parisi, G.; Zamponi,
F. Dimensional dependence of the Stokes-Einstein relation and its
violation. J. Chem. Phys. 2013, 139, 164502.
(64) Charbonneau, P.; Jin, Y.; Parisi, G.; Zamponi, F. Hopping and
the Stokes-Einstein relation breakdown in simple glass formers. Proc.
Natl. Acad. Sci. U.S.A. 2014, 111, 15025.
(65) Mallamace, F.; Branca, C.; Corsaro, C.; Leone, N.; Spooren, J.;
Chen, S.-H.; Stanley, H. E. Transport properties of glass-forming
liquids suggest that dynamic crossover temperature is as important as
the glass transition temperature. Proc. Natl. Acad. Sci. U.S.A. 2010,
107, 22457.
(66) Yeh, I.-C.; Hummer, G. System-Size Dependence of Diffusion
Coefficients and Viscosities from Molecular Dynamics Simulations
with Periodic Boundary Conditions. J. Phys. Chem. B 2004, 108,
15873.
(67) Tazi, S.; Boţan, A.; Salanne, M.; Marry, V.; Turq, P.; Rotenberg,
B. Diffusion coefficient and shear viscosity of rigid water models. J.
Phys.: Condens. Matter 2012, 24, 284117.
(68) Abascal, J. L. F.; Vega, C. Widom line and the liquid-liquidcritical point for the TIP4P/2005 water model. J. Chem. Phys. 2010,
133, 234502.
K

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B(69) Hestand, N. J.; Skinner, J. L. Perspective: Crossing the Widom
line in no man’s land: Experiments, simulations, and the location of
the liquid-liquid critical point in supercooled water. J. Chem. Phys.
2018, 149, 140901.
(70) Ni, Y.; Hestand, N. J.; Skinner, J. L. Communication: Diffusion
constant in supercooled water as the Widom line is crossed in no
man’s land. J. Chem. Phys. 2018, 148, 191102.
(71) Bhattacharyya, S.; Mukherjee, A.; Bagchi, B. Correlatedorientational and translational motions in supercooled liquids. J.
Chem. Phys. 2002, 117, 2741.
(72) Takeuchi, H. A jump motion of small molecules in glassy
polymers: A molecular dynamics simulation. J. Chem. Phys. 1990, 93,2062.
(73) Müller-Plathe, F.; Rogers, S. C.; van Gunsteren, W. F.
Computational evidence for anomalous diffusion of small moleculesin amorphous polymers. Chem. Phys. Lett. 1992, 199, 237−243.
(74) Tocci, E.; Hofmann, D.; Paul, D.; Russo, N.; Drioli, E. A
molecular simulation study on gas diffusion in a dense poly(etherether-ketone) membrane. Polymer 2001, 42, 521.
(75) Indra, S.; Daschakraborty, S. Mechanism of translational jumpof a hydrophobic solute in supercooled water: Importance ofpresolvation. Chem. Phys. Lett. 2017, 685, 322−327.
(76) Kumar, A. V.; Bhatia, S. K. Mechanisms Influencing Levitationand the Scaling Laws in Nanopores: Oscillator Model. J. Phys. Chem.
B 2006, 110, 3109−3113.
(77) Raptis, T. E.; Raptis, V. E.; Samios, J. New effective method forquantitative analysis of diffusion jumps, applied in molecular
dynamics simulations of small molecules dispersed in short chain
systems. J. Phys. Chem. B 2007, 111, 13683−13693.
(78) Raptis, T. E.; Raptis, V. E.; Samios, J. Quantitative study ofdiffusion jumps in atomistic simulations of model gas-polymer
systems. Mol. Phys. 2012, 110, 1171−1178.
(79) Araque, J. C.; Yadav, S. K.; Shadeck, M.; Maroncelli, M.;
Margulis, C. J. How Is Diffusion of Neutral and Charged Tracers
Related to the Structure and Dynamics of a Room-Temperature Ionic
Liquid? Large Deviations from Stokes-Einstein Behavior Explained. J.
Phys. Chem. B 2015, 119, 7015−7029.
(80) Shell, M. S.; Debenedetti, P. G.; Stillinger, F. H. Dynamic
heterogeneity and non-Gaussian behaviour in a model supercooled
liquid. J. Phys.: Condens. Matter 2005, 17, S4035−S4046.
(81) Rahman, A. Correlations in the Motion of Atoms in Liquid
Argon. Phys. Rev. 1964, 136, A405−A411.
(82) Laage, D.; Hynes, J. T. A molecular jump mechanism of water
reorientation. Science 2006, 311, 832−835.
(83) Laage, D.; Hynes, J. T. On the molecular mechanism of water
reorientation. J. Phys. Chem. B 2008, 112, 14230−14242.

L

DOI: 10.1021/acs.jpcb.9b01719
J. Phys. Chem. B XXXX, XXX, XXX−XXX