Isotopic Effects on Intermolecular and Intramolecular Structure and

Aug 10, 2018 - Rationalization of their physical origins and the obtained physical insights will help in experimentally searching and monitoring not o...
0 downloads 0 Views 1MB Size
Subscriber access provided by University of Sussex Library

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Isotopic Effects on Intermolecular and Intramolecular Structure and Dynamics in Hydrogen, Deuterium and Tritium Liquids: Normal Liquid and Weakly and Strongly Cooled Liquids Kiharu Abe, Shutaro Yamaoka, and Kim Hyeon-Deuk J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b02596 • Publication Date (Web): 10 Aug 2018 Downloaded from http://pubs.acs.org on August 16, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Isotopic Effects on Intermolecular and Intramolecular Structure and Dynamics in Hydrogen, Deuterium and Tritium Liquids: Normal Liquid and Weakly and Strongly Cooled Liquids Kiharu Abe,†,‡ Shutaro Yamaoka,† and Kim Hyeon-Deuk∗,†,¶ Department of Chemistry, Kyoto University, Kyoto, 606-8502, Japan, Present Address:Graduate School of Medical Life Science, Yokohama City University, Kanagawa, 230-0045, Japan, and Japan Science and Technology Agency, PRESTO, 4-1-8 Honcho, Kawaguchi, Saitama, 332-0012, Japan E-mail: [email protected],Phone:+81-75-753-4021



To whom correspondence should be addressed Department of Chemistry, Kyoto University, Kyoto, 606-8502, Japan ‡ Present Address:Graduate School of Medical Life Science, Yokohama City University, Kanagawa, 2300045, Japan ¶ Japan Science and Technology Agency, PRESTO, 4-1-8 Honcho, Kawaguchi, Saitama, 332-0012, Japan †

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32

Abstract

Differences in properties such as phase-transition temperature and transport coefficients among liquids of different isotopic compositions, hydrogen, deuterium and tritium, should originate from their differently pronounced nuclear quantum effects (NQEs) rather than from any subtle difference in the electronic interaction potentials. Accurate and efficient determination of structural and dynamical isotopic effects in the quantum liquids still remains as one of the challenging problems in condensed-phase physics. With a recently developed non-empirical real-time molecular dynamics method which describes non-spherical molecules with the NQEs, we computationally realized and investigated dynamical and quantum isotopic effects of not only traditionally studied isotopes, hydrogen and deuterium, but also a lesser-known radioisotope, tritium, in broad thermodynamic conditions from normal liquid to weakly and strongly cooled liquids which have been hindered by rapid crystallization in spite of numerous experimental attempts at supercooling. Reproducing the previously reported experimental isotope dependence on the bond length and vibrational frequencies of hydrogen, deuterium and tritium liquids, we further demonstrate that distinctive isotope effects appear in their intermolecular and intramolecular structure and dynamics not only at lower temperature but also at higher temperature, which none has so far been able to obtain quantitative results for realistic systems. Rationalization of their physical origins and the obtained physical insights will help future experimental searching and monitoring intermolecular and intramolecular dynamics and structures of these isotopes not only in normal liquid but also in supercooled liquid.

2

ACS Paragon Plus Environment

Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction Understanding microscopic dynamics of condensed-phase systems exhibiting significant quantum effects is still one of the open problems in condensed matter physics. Although hydrogen and its isotopes are the simplest of all molecular species, nuclear quantum effects (NQEs) are expected to be noticeable, making them mysterious and interesting. Typical values for the well known de Boer parameter, which effectively scales the importance of NQEs, are 0.28 for H2 and 0.20 for D2 , which are one order of magnitude larger than the typical values for classical behavior. 1–4 In fact, the NQEs of liquid hydrogen appear as a broader distribution of radial distribution functions (RDFs), large self-diffusion coefficients even at low temperature, and anomalous temperature dependence of thermal conductivity and shear viscosity. 5–11 One should bear in mind that systems obey Boltzmann statistics and therefore assessing extent of the influence of the NQEs is anything but trivial. Differences in thermodynamic properties such as phase-transition temperature and transport coefficients among liquids of the different isotopic compositions should be ascribed to their pronounced NQEs in the order of of H2 , D2 and T2 as well as the mass difference rather than to any subtle difference in the electronic interaction potentials. 12,13 Static and dynamical structure factors as well as coefficients of self-diffusion in liquid deuterium have been measured by neutron diffraction experiments. 3,14–21 Those experimental results confirmed inadequacy of classical description of nuclei and required treating nuclei quantum mechanically taking into account the NQEs. 22 Several quantum mechanical approaches have been employed to elucidate the anomalous tendency observed in the quantum liquids. Semiquantum methods suitable for numerical simulations of p-H2 and o-D2 ensembles have been proposed based on the imaginary-time path-integral theory. 5,8–11,23–29 These provide either static properties or equilibrium time correlation functions(TCFs) accounting for particularly important NQEs such as zero-point energy effect. 5,6,8–11,23,27 Experimental and computational data of hydrogen and deuterium have been compared in 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

many physical quantities and situations, e.g. static and dynamical structure factors, 1,13,16,30–32 self-diffusion coefficients, 33 nuclear delocalization, 2 high-pressure solids, 34 crystallization processes, 2,35 quantum melting in a small cluster, 36 absorption to graphene, 37 observing quantitatively significant differences. Based on the fact that hydrogen and deuterium experience the same electronic potential, the differences observed and found can only be attributed to the different isotopic behaviors of the two molecules. We recently proposed a simulation method of nuclear and electron wave packet molecular dynamics (NEWPMD) based on non-empirical ab initio intramolecular and intermolecular interactions of non-spherical molecules. 38,39 The NWPs were non-perturbatively treated in the NEWPMD which can be distinguished from the previous NWP approaches in which a potential surface was given in advance by a separate modeling and, in many cases, expanded quadratically around the moving NWP centers to perturbatively take into account the NQEs. 23,40–43 While each of nucleus and electron was approximately described by a time-dependent thawed Gaussian wave packet, the NEWPMD gave the correct structures and transport properties of liquid p-H2 under vapor pressure without any empirical parameters. 39,44 It also successfully reproduced the stable vapor-pressure solid of the hexagonalclose-packed (h.c.p.) lattice structure with the reasonable freezing temperature and lattice phonon frequencies. 45 Even the non-equilibrium state, heat conduction under a temperature gradient, was directly achieved and the reasonable thermal conductivity of p-H2 liquid was reproduced. 46 The NEWPMD method is thus distinguished from most of the previous pathintegral approaches for a condensed phase which require an equilibrium state through an equilibrium partition function as well as an interaction potential surface given in advance by a separate modeling with empirical parameters like the Silvera-Goldman potential. 47 In this paper, we report static and dynamical isotope effects, which stem from different masses and different extent of nuclear quantumness, appearing in intermolecular and intramolecular properties of p-H2 ,o-D2 ,and p-T2 liquids. None has so far been able to computationally or experimentally find their quantitative results for realistic systems. We will

4

ACS Paragon Plus Environment

Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

investigate broad thermodynamic situations focusing on not only normal liquid (NL) but also weakly and strongly cooled liquid (WCL and SCL) which have been hindered by rapid crystallization in spite of numerous experimental attempts at supercooling. 2,48–56 For simplicity, p-H2 , o-D2 and p-T2 will be simply written as H2 , D2 and T2 throughout this paper.

Computational Details The NEWPMD method describes nuclei by floating and thawed Gaussian WPs via the time-dependent Hartree approach, and electron wave packets (EWPs) by the perfect-pairing valence bond theory that appropriately treats the Pauli exclusion energy and intermolecular dispersion energy. 39,44 The width of each EWP linearly changes depending on the bond length of each molecule, which was deduced from the calculation of a single molecule in the same frame work of the NEWPMD approach. 38 We adopt these thawed EWPs to calculate all of the intramolecular potentials while the EWPs for the intermolecular potentials are fixed with the most stable width. The electronic part of the NEWPMD program is exactly the same regardless of the three isotopes while the atomic masses of H2 , D2 and T2 are set as 1.00794, 2.01410 and 3.01605, respectively. All the H2 , D2 , and T2 systems are composed of 640 molecules in a three-dimensional(3D) cubic simulation box with a periodic boundary condition. We examine three distinctive thermodynamic conditions: NL at relatively high temperature (25 K), WCL at slightly lower temperature than the freezing point of 18.7 K (18 K), and SCL at much lower temperature (5 K). Density of the NL of H2 , D2 , and T2 at 25 K was set as the vapor-pressure density of D2 liquid at 25 K: the molar volume is 25.17 × 10−6 m3 /mol. Densities of the WCL and the SCL were set by extrapolating the saturated vapor-pressure line of D2 down to 18 K and 5 K, respectively: 10,12 the molar volume is 23.20 × 10−6 m3 /mol at 18 K and 22.50 × 10−6 m3 /mol at 5 K which are both larger than the deuterium solid molar volume. 12 We can of course set different densities depending on the isotopes. However, the current purpose is extracting

5

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

and elucidating the pure isotope effects which appear even in the same thermodynamic condition and which are free from the density difference. The different density directly leads to different structures, and thus naturally causes different dynamical properties like different diffusion coefficients, making it difficult to distinguish difference stemmed from the pure isotope effects from difference stemmed from simple structural difference. We removed such trivial difference by matching the thermodynamic conditions of H2 and T2 liquids with the conditions of D2 liquid. We started cooling and equilibration runs from a h.c.p. crystal structure of the above molar volume where molecules have random orientation. We made atomic center momentum degrees of freedom influenced by a heat bath set by the velocity scaling thermostat and then by the Berendsen thermostat for hundreds of picoseconds at a given temperature. All integrations of the equations of motion were performed by the velocity-verlet method with the time step 0.1 fs for the cooling and equilibration runs. After the careful equilibration runs, we carried out the NVE (microcanonical) simulations for hundreds of picoseconds with the time step of 0.5 fs. We confirmed that the temperature is stably fluctuating around the aimed temperature in the all NVE simulations (See Fig. S1). The corresponding views of the 3D cubic simulation boxes containing 640 molecules are shown in Fig. S2.

Results and Discussions Isotopic effects on liquid structures Figures 1(a)-(c) characterize liquid structures of the three isotopes through RDFs of the NL, WCL and SCL,

g(r) =

hn(r)i , 4πr2 drn0

6

ACS Paragon Plus Environment

(1)

Page 7 of 32

25 K

5K

18 K 2

RDF

(a) H2 D2 T2

(b) H2 D2 T2

1

(c) H2 D2 T2

2

1 1

0 0

2

4

6

8

10

12

14

0 0

2

4

6

8

10

12

14

0 0

2

4

6

8

10

12

14

o

r (A) (e)

(d) Distribution

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(f)

0.03

0.03

0.03

0.02

0.02

0.02

-1

0

1

-1

0 cos θ

1

-1

0

Figure 1: RDFs of (a) the NL, (b) the WCL and (c) the SCL of H2 , D2 and T2 . Angular distributions of (d) the NL, (e) the WCL and (f) the SCL of H2 , D2 and T2 . As graphically defined, the intermolecular angle θ is a relative angle formed by a pair of molar vectors of two molecules. Only the lower temperature induces the isotopic difference. where n(r) represents the number of NWPs between shells of radii r and r + dr, and n0 is the number density of the whole system. h· · ·i means an ensemble average over different initial times and NWPs. No long-range periodic structure was obtained in the current RDFs unlike the RDFs of the H2 solid. 45 The whole RDF shapes are rather close to the RDFs of the normal H2 liquid, 44 supporting the present liquidity. Achieving the WCL and SCL should be largely owed to the careful equilibration process in the current simulation: setting unsuitable density for crystallization, rapidly cooling the whole system and using a noncommensurate simulation cell. These protocols made it possible to frustrate formation of a crystal even at 5 K. 57 The RDFs of the NL and WCL exhibit little difference regardless of H2 , D2 or T2 . As far as the thermodynamic condition is exactly same and as far as the electronic potentials created by the EWPs are identical, the RDFs of the higher-temperature NL and WCL, where

7

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

difference of dynamics depending on the isotopes does not significantly affect the structure, should be similar regardless of the isotopes. This is a simple equilibrium probability distribution over the free energy landscape achieved by the relatively narrow NWP width which will be shown in Fig.7. In fact, taking into account the probability of the hydrogen nuclear delocalization, which is broadest among the current isotopes, in calculating the RDF does not significantly affect the result. 39 The small difference of 0.01 ˚ A among the H-H, D-D and T-T bond length also contributed to the current similar RDFs regardless of the isotopes as will be shown in Figs.6(a)-(c). The RDFs at the lower temperature have higher peaks and lower valleys than the RDFs at the higher temperature, indicating more structural ordering at the lower temperature. The shoulder appearing in the left part of the first RDF peak becomes clearer at 18 K than at 25 K, and the first RDF peaks are almost double-split at 5 K. The shoulder or split reflects a diatomic structure of H2 , D2 and T2 molecules in the first solvation shell. Because the most stable solvation structure of a pair of hydrogen molecules is a T-shape configuration, its nearest atoms appear as the shoulder or split at the lower temperature reflecting lower thermal modulations of molecules in the first solvation shell. 44 (See Fig. S3) The isotope effect can be clearly seen in the double-split RDFs of the SCL at the lower temperature. While the width of the first RDF peak is similar regardless of H2 , D2 , or T2 , the peak position was shifted toward the shorter distance in the H2 SCL than in the D2 and T2 SCL. The closer first peak of the H2 SCL should be related to the faster freezing of H2 molecules than the D2 freezing observed by Raman spectroscopy. 2 As will be shown in Fig.2, lighter H2 molecules have larger mobility and can fit into the structure of the first-solvation shell more easily and smoothly. On the contrary, the first peak height of the D2 and T2 RDFs is higher than the height of H2 , indicating that the H2 SCL is more liquidized compared to the D2 and T2 SCL in spite of their same electronic potentials and the same thermodynamic condition. The softer structure of the H2 SCL is also supported by the structural order parameter defined as the area of the left sub-peak in the split first

8

ACS Paragon Plus Environment

Page 9 of 32

peak divided by the area of of the right sub-peak. Such order parameter can directly reflect stiffness of a T-shape coordination formed by nearest-surrounding molecules. The calculated structural order parameters are smaller in the order of H2 , D2 and T2 , i.e., 0.259 for H2 , 0.363 for D2 and 0.364 for T2 . The NEWPMD method can describe molecular orientational structures as well as the radial structures. Figures 1(d) and 1(e) show that the angular distributions of the NL and the WCL are almost uniform regardless of the isotope species. The flat distributions mean that molecules are free rotors in the whole system and guarantee the liquidity of the NL and WCL. 34 While no difference among the isotopes is observed in the NL and WCL, the SCL again exhibits the clear isotope effect. As Fig.1(f) shows, molecular angles become more ordered at the lower temperature in the order of T2 , D2 and H2 , which is in harmony with the larger structural order parameters in the order of T2 , D2 and H2 . The SCL has metastable and jammed structures, and the less translational dynamics due to the heavier mass leads to the less uniform angular distributions. The more liquidized structure of lighter molecules indicates that the quantum melting occurs also in the current SCL not only in a small cluster or in a high-pressure solid. 34,36

Isotopic effects on diffusive dynamics

20 o2

MSD (A )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) 25K H2 D2 T2

10

(b) 18K H2 D2 T2

0.5

(c) 5K H2 D2 T2

10

0

0

1

2

3

4

0 5 0

2

4 6 Time (ps)

8

10

0.0 0

50

100

150

200

250

Figure 2: MSDs of (a) the NL, (b) the WCL and (c) the SCL of H2 , D2 and T2 . The MSDs exhibit the isotope difference even at the higher temperature. The mean square displacement (MSD), D(t) = h|R(t) − R(0)|2 i, effectively character9

ACS Paragon Plus Environment

300

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

izes time-dependent diffusive dynamics. 58–60 Here, R(t) denotes the position vector of a NWP center at time t. 42,43 The MSDs of the NL and WCL shown in Figs.2(a) and 2(b), respectively, almost linearly depend on a time. The nonlinear MSDs appearing within the initial 1 ps in the WCL case correspond to initial dynamics confined in a stiff first-solvation cage. The estimated self-diffusion coefficients, obtained by taking limt→∞ D(t)/6t (i.e. Einstein’s relation), are listed in Table 1. The self-diffusion coefficient of D2 at 25 K, 0.505 ˚ A2 /ps, is reasonable taking into account the fact that the current diffusion coefficient would increase by approximately 20% in a infinite system as was seen in the normal H2 liquid cases. 12,18,18,23,39,44,61 Although the RDFs and angular distributions did not exhibit any isotope difference at 25 K and 18 K, the MSDs and self-diffusion coefficients become larger in the order of H2 , D2 and T2 . The highest mobility of H2 essentially stems from its lightest mass. There are two counter properties to decelerate H2 molecules compared to D2 and T2 molecules: (1) H-H bond is longer than D-D and T-T bonds as will be explained in Figs.6(a)-(c); a self-diffusion coefficient becomes greatly smaller with larger bond length of molecules in general. 39,44 (2) Hydrogen wave packets (HWPs) are broader than deuterium wave packets (DWPs) and tritium wave packets (TWPs) as will be explained in Fig.7; larger nuclear volume suppresses diffusion. The lightest mass of an H2 molecule overcame these disadvantages and contributed to the largest self-diffusion coefficients of H2 ; the effects of the longer H-H bond and the broader HWP are negligible for the self-diffusion coefficients. We note that, if the masses of the isotope species were the same, T2 would have highest mobility, which is in harmony with the previous report that H2 and D2 molecules with classical nuclei have larger self-diffusion coefficients than H2 and D2 molecules with quantum nuclei. 62 Simply multiplying the H2 diffusion coefficient by the square root of the inverse ratio of the masses predicts diffusion coefficients close to the directly calculated values listed in Table 1; e.g. the self-diffusion coefficients of D2 at 25 K and at 18 K were obtained by the mass scaling as 0.478 ˚ A2 /ps and 0.172 ˚ A2 /ps, respectively. However, the temperature dependence of the self-diffusion coefficients is isotope-dependent now; as Table 1 shows, the ratio of the

10

ACS Paragon Plus Environment

Page 11 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

self-diffusion coefficient at 25 K to the self-diffusion coefficient at 18 K is larger in the order of T2 , D2 and H2 . The isotope-dependent ratios at the two temperatures were caused by the WCL. The heavier molecule, T2 , is more influenced in the WCL compared to the NL case; the T2 mobility is more suppressed in the WCL compared to the mobility of lighter H2 and D2 molecules, leading to the slightly isotope-dependent ratios. The MSD of the SCL shown in Fig.2(c) has a few characteristic stages and was not on a simple line; the diffusion dynamics is dominated by mesoscale cluster deformations in the ˚2 corresponds to SCL. 60 The initial instant increase of the MSD up to approximately 0.2 A phonon motions inside a mesoscale lattice cluster which persistently remains in the SCL and should be called local-β relaxation. 60 The following plateau stage continues for approximately 100 ps, indicating metastable phonon modes formed in a mesoscale lattice cluster. Mesoscale clusters slowly deform and collapse as appearing in the slow increase of the MSD after the long-time plateau stage, which corresponds to α-relaxation. 60 These transitions between the plateau and relaxation stages become more significant in the order of H2 , D2 and T2 . Mesoscale clusters composed of lighter H2 molecules can deform more significantly because of the higher mobility. Table 1: Self-Diffusion Coefficients Estimated by Linear Fitting of the MSD with the Unit 2 ˚ A /ps. Isotope H2 D2 T2 a

25 K 0.680(1.10 × 10−4 )a 0.505(2.58 × 10−5 ) 0.395(1.21 × 10−4 )

18 K 0.243(4.21 × 10−5 ) 0.179(1.88 × 10−5 ) 0.131(2.40 × 10−5 )

Ratio 2.80 2.82 3.02

The numbers written in the parenthesis are the standard error of the mean.

Isotopic effects on microscopic angular dynamics We further discuss how an isotopic effect appears in microscopic angular dynamics of nonspherical molecules. Figures 3(a)-(c) show power spectra obtained by Fourier transformation ∫

of a time-dependent intermolecular angle θ(t), ∼ h| dtθ(t) exp (−iωt) |2 i. All the power spec11

ACS Paragon Plus Environment

The Journal of Physical Chemistry

10 (b) 18K H2 D2 T2

5

(c) 5K H2 D2 T2

2

240

5 1

200 180 160 140 120 100

0

200

400

600

0

200

600 0

400

200

-1

Frequency (cm )

400

(d) H 2 D2 T2

220 -1

(a) 25K H2 D2 T2

Frequency (cm )

Intensity (arb. unit)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

600

80

5

θ 10

15 20 Temperature (K)

Figure 3: Power spectra of the intermolecular angle θ in (a) the NL, (b) the WCL and (c) the SCL of H2 , D2 and T2 . (d) Librational frequencies are blueshifted in the order of H2 , D2 and T2 regardless of the temperature. tra have a peak near 200 cm−1 which can be attributed to librational dynamics of molecules rather than to rotational motions. 44,60 Temperature dependence of the peak frequency and intensity demonstrates the current librational picture; a rigid T-shape solvation structure at the lower temperature provide a sharper librational potential leading to the higher peak frequency of librations regardless of the isotopes.(See also Fig.3(d)) At the higher temperature, the T-shape configuration is softened, resulting in the redshifted librational motions. In addition, the power spectral intensity for the higher temperature is larger than the intensity at the lower temperature due to the stronger thermal fluctuations of librations. Note that the rotational region of Raman spectra is almost independent of the temperature, 55 and, further, the current temperature is too low to cause larger-amplitude and unidirectional rotational motions. The librational power spectra are clearly isotope-dependent at all the temperature in spite of the same structures at 25 K and 18 K shown in Figs.1(a) and 1(b). The librational peaks as well as the high-frequency tails become totally blueshifted and broader in the order of H2 , D2 and T2 . Smaller moments of inertia of the lighter molecules leads to the higher-frequency librational peaks. The broader HWP compared to the DWP and TWP also contributes to the blue shift since the broader NWPs can hinder librational motions due to the additional repulsive force between the broader WPs, which can be called a quantum isotopic effect. Actually, as Fig.3(d) shows, the librational frequency becomes more blueshifted with decreasing the temperature from 25 K to 18 K in the order of H2 , D2 and T2 . At 5 12

ACS Paragon Plus Environment

25

Page 13 of 32

K, the H2 liquid structure, which is more softened compared to the liquid structures of D2 and T2 as shown in Figs.1(c) and 1(f), compensates the current quantum isotopic effect, resulting in the similar temperature dependence regardless of the isotopes with decreasing the temperature from 18 K to 5 K. 1.0

1.0

1.0

(a) 25K H2 D2 T2

TCF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(b) 18K H2 D2 T2

(c) 5K H2 D2 T2 0.5

0.5

0

1

2

3

4

5

0.5

0

1

2 3 Time (ps)

4

5

0.0

0

1

2

3

4

Figure 4: TCFs of the intermolecular angle θ in (a) the NL, (b) the WCL and (c) the SCL of H2 , D2 and T2 . The TCFs decay faster in the order of H2 , D2 and T2 regardless of the qualitative change of the decay profiles. Scaled TCFs of θ(t), hθ(t)θ(0)i scaled by the initial value at t = 0, also exhibit isotopic difference.(Fig.4) The TCFs of the NL and of the WCL have the distinguishable shorttime and long-time decays. Their decay timescales obtained by fitting the TCFs with an exponential function are listed in Table 2. Note that the present biexponential decay of the librational dynamics is too complicated to describe by the simple Debye model developed for a rotational motion which exhibits a single exponential correlation decay. A ratio of the shorter timescales, e.g. 208 fs over 152 fs, is close to a ratio of the temperature, 25 K over 18 K, suggesting that the short timescales are related to kinetic dynamics of a molecule inside a first solvation cage. The other ratio of the longer timescales, e.g. 1250 fs over 400 fs, is more sensitive to the temperature and rather close to a ratio of the self-diffusion coefficients 2 2 for 25 K over 18 K, 0.680 ˚ A /ps over 0.243 ˚ A /ps. 39,44 This coincidence indicates that the

longer timescales originated from the diffusive dynamics of a molecule out of a solvation shell. The isotopic effect appears as the faster decay of the TCFs for the lighter molecules reflecting the higher mobility. The TCFs of the SCL have the similar isotopic effect while the decay is much faster and the decay profiles are not simply exponential but oscillating, which 13

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry

again indicates the different mechanism for diffusion appearing in the SCL. Relatively stable angular dynamics inside a mesoscale cluster formed at 5 K resulted in the faster short-time decay in spite of the smaller thermal modulation. The above insights are further supported by the power spectra and TCFs of time-dependent self-orientation of each molecule, φ(t), which behaves quite similarly to the intermolecular angle θ(t) . (See Figs.S4 and S5) Table 2: Decay Timescales for the Librational Dynamics of the Intermolecular Angle θ Estimated from Fig.4. Temperature 25 K Isotope H2 D2 Short Range (fs) 152 200 Long Range (fs) 400 435

T2 H2 250 208 625 1250

18 K D2 278 1429

T2 370 2000

Isotopic effects on low-energy collective dynamics: Boson peak

20 (a) 25K H2 D2 T2

600 Intensity (arb. unit)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32

(c) 5K H2 D2 T2

(b) 18K H2 D2 T2

200

400 10 100 200

0

0

20

40

60

80

100 120 140

0

0

20

40

60

80

100 120 140

0 0

20

40

60

80 100 120 140

-1

Frequency (cm )

Figure 5: Power spectra of the low-energy collective dynamics in (a) the NL, (b) the WCL and (c) the SCL of H2 , D2 and T2 . The power spectra become broader to higher frequency and the peaks appearing at the lower temperature are blueshifted in the order of H2 , D2 and T2 . Figure 5 displays power spectra obtained by Fourier-transforming a real-time coordinate displacement of each atom in the space-fixed frame. Such low-frequency power spectra correspond to low-energy collective dynamics involved in the isotope liquids. 60 The power 14

ACS Paragon Plus Environment

Page 15 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

spectra of the NL have a simple long tail while the power spectra of the WCL and SCL at the lower temperature exhibit a broad low-frequency peak which can be assigned as a boson peak analogous to the peaks observed in various glasses and supercooled liquids. 58,59,63 Since the boson peak is closely related to metastability of a mesoscale lattice cluster, 60 the appearance of the boson peaks in the WCL demonstrates that the WCL is not a normal liquid either like the SCL. In fact, the frequency of the boson peak redshifts with increasing the temperature as shown in Fig.S6; the metastable lattice clusters become more fragile at the higher temperature. There is another peak around 5 cm−1 which can be attributed to translational dynamics of molecules; the peak becomes enhanced in the NL and the WCL with increasing the temperature while it disappears in the SCL at the lower temperature where significant translational motions are almost frozen. The isotopic difference can be seen throughout all the temperature ranges. The lowfrequency power spectra become broader toward higher frequency and the peak positions are blueshifted in the order of H2 , D2 and T2 .(See also Fig.S6) The higher-frequency peaks were induced by smaller mass of lighter molecules just like the librational peaks shown in Fig.3(d). However, the quantum isotopic effect caused by the broader HWP compared to the DWP and TWP does not play an important role in the low-energy collective dynamics; as shown in Fig.S6, the temperature dependence of the peak frequencies does not change depending on the isotopes, which is different from the librational power spectra shown in Fig.3(d). The microscopic WP width of each nucleus does not affect the current collective dynamics involving many molecules. In addition, while the librational peak intensity does not change by the isotopes, the current intensity of the low-energy collective dynamics increases in the order of H2 , D2 and T2 especially in the NL and WCL at the higher temperature. The NL and WCL involve significant translational dynamics and their intensity of the low-energy collective dynamics is influenced by the molecular mobility depending on the isotopes.

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry

0.007

0.07

0.01

(a) H2

25 K

(b) H2

18 K

Distribution

D2 T2

0.000 0.72

0.74

0.76

(c) H2

5K

D2 T2

D2 T2

0.00 0.72

0.78

0.74

0.76

0.00 0.72

0.78

0.74

0.76

0.78

o

Bond Length (A) 4606 Frequency (cm-1 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32

3354

(d)

2775

(e)

4604

3352

2773

4602

3350

2771

4600

3348

2769

H2 45985

10

15

20

25

3346

(f)

T2

D2 5

10

15

20

25

2767

5

10

15

20

Temperature (K)

Figure 6: Distributions of H-H, D-D and T-T bond length in (a) the NL, (b) the WCL and (c) the SCL. The length range of each figure is the same. The distribution of T-T bond length exhibits split double peaks. (d)-(f) Temperature dependence of vibrational bond frequencies of H2 , D2 and T2 . The frequency range of each figure is the same. The temperature dependence is more steep in the order of H2 , D2 and T2 .

16

ACS Paragon Plus Environment

25

Page 17 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Isotopic effects on intramolecular structures and dynamics The NEWPMD approach includes intramolecular degrees of freedom of a non-spherical molecule, and can describe detailed microscopic structures and dynamics of each molecule even in a condensed phase. Figures 6(a)-(c) show distributions of H-H, D-D, and T-T bond length calculated using the position vector of each NWP center R(t) in the NL, WCL and SCL. The distribution becomes sharper with decreasing the temperature regardless of the isotope species reflecting the smaller thermal fluctuations at the lower temperature. It is remarkable that the shape of the T-T bond length distribution in the SCL has split double peaks because of the following reasons: (1) The metastable structure of the SCL involving deformation and collapse of mesoscale lattice clusters is quite irregular especially in the T2 case as shown in Figs.1(c) and 1(f). (2) The lower mobility of T2 molecules prevents intercommunication between stiff T-shape solvation in the mesoscale cluster and weak solvation out of the cluster region. (3) A T-T Bond is effectively elongated by attraction from EWPs of adjacent T2 molecules in a T-shape shell as the system becomes more condensed at the lower temperature. (See Fig.S7.) These facts result in the larger difference between long and short T-T bond length, rationalizing the split double peaks. The split double peaks reflect that the liquid has two typical structures, a well-coordinated solvation structure inside a solvation shell and a poorly-coordinated structure out of a solvation shell, as was reported in Ref.44 where the H-H vibrational power spectra of the hydrogen liquid required two functions to be fitted. Because the H-H vibrational frequency highly depends on the surrounding solvation structures, the necessity of the two functions itself indicates the existence of the two typical solvation structures in the hydrogen liquid. We also successfully fitted all the current vibrational spectra of H2 , D2 and T2 only with two functions as shown in Fig. S8. 44 There is another isotopic effect in the bond length distributions: The T-T bond length distributions are broader than the H-H and D-D bond length distributions also in the NL and the WCL. Since the liquid structures in the NL and WCL are almost the same independent 17

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

of the isotopes, the current isotopic difference should be caused by another mechanism: The TWPs narrower than the HWPs and DWPs, i.e. the tritium nuclear volume smaller than the hydrogen and deuterium volume, induce stronger attraction of a T2 molecule from surrounding T2 molecules even in the same liquid structures while T2 molecules become more shrunk out of a solvation shell due to a weaker intramolecular repulsive force between smaller TWPs, rationalizing the broader distributions of T-T bond in spite of the same liquid structures. (See Figs.7(a)-(c).) The present quantum isotope features determine the intramolecular structure even at the higher temperature not only at the lower temperature. As shown in Fig.S7, the obtained average bond length as well as the isotope dependence that the bond length is longer in the order of H2 , D2 and T2 are reasonably compared to the experimental results. Note that the average bond length in a liquid phase is generally longer than bond length of an isolated free molecule; because the most stable solvation structure is the T-shape configuration, the bond of a molecule is effectively elongated by the intermolecular attractive interaction from adjacent molecules as the system becomes more condensed. 44,64 The average bond length of isolated H2 , D2 and T2 molecules, which is purely determined by the intramolecular interaction including the anharmonicity and completely free from the intermolecular interactions, was calculated as 0.76429 ˚ A, 0.75163 ˚ A and 0.74659 ˚ A, respectively: the experimental bond length in a gas phase is 0.7510 ˚ A for H2 , 0.7482 ˚ A for D2 and 0.7469 ˚ A for T2 . 38,64 If there was no isotope effect in the bond length, all the bond length would be the same due to the same electronic potentials. The temperature dependence of the bond length is also isotope-dependent; the bond length changes more rapidly with the temperature in the order of H2 , D2 and T2 because the lighter molecule is more easily elongated with a same change of the electronic force compared to the heavier molecule. This isotopic mass effect overcomes the counter property that the heavier molecule achieves the stronger structurization than the lighter molecule in the SCL. Such trend in the isotope-dependent temperature dependence of the bond length is kept throughout all the temperature range.

18

ACS Paragon Plus Environment

Page 19 of 32

Bond frequencies shown in Fig.6(d)-(f) have the similar isotope dependence to the bond length but exhibit the opposite temperature dependence against the bond length; the bond frequencies become smaller with decreasing the temperature. As the bond is more stretched in the more rigid structure at the lower temperature, the bond becomes weakened and the vibrational frequency is redshifted. The bond frequencies were estimated from the peak frequencies of the power spectra obtained by a Fourier transformation of time-dependent bond length.(See Fig.S8) We note that the current bond frequencies of H2 , D2 , and T2 overestimate the experimental vibrational frequencies by approximately 11 % that would not significantly affect other important properties such as radial distribution functions and diffusion coefficients reported here. 22,24,44

Distribution

0.02

(a) 25K HWP DWP TWP

0.03

(b) 18K HWP DWP TWP

(c) 5K HWP DWP TWP 0.04 0.05

0.02

0.03

0.01 0.02

0.01

0.01 0.00

0.04

0.05

0.06

0.07

0.00 0.04

0.05

0.06

0.07

0.00 0.04

0.05

0.06

0.07

o

NWP Width (A)

0.065765 0.04960 0.054754

o

NWP width (A)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

0.04958 0.065755 0.04956

0.054750

0.04954 0.065745 5

(d) HWP 10

15

20

0.054746 5 25

(e) DWP 10

(f) TWP 15 20 Temperature (K)

25

5

10

15

20

Figure 7: Distributions of the HWP, DWP and TWP width in (a) the NL,(b) the WCL and (c) the SCL. The distributions are shown in the same length range. The TWP exhibits anomalously broadened double peaks even at the high temperature. (d)-(f) Temperature dependence of width of the HWP, the DWP and the TWP. While the WP width is larger in the order of the HWP, DWP and TWP, the TWP exhibits more sensitive and qualitatively different temperature dependence. The error bars estimated by the standard errors of the mean indicate that all the average WP width converged almost within the shown digits. Figures 7(a)-(c) show extent of nuclear delocalization depending on the isotopes as distri19

ACS Paragon Plus Environment

25

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

butions of the NWP width. The distributions become sharper as the temperature decreased from 25 K down to 5 K reflecting the poorer thermal modulations at the lower temperature. The isotopic difference appears in the distribution width which is broader in the order of T2 , D2 and H2 at any temperature. Especially, only the distributions of the TWP width have split double peaks even at the high temperature. The anomalously broadened double peaks can be rationalized by the broad distributions of the T-T bond length discussed in Figs.6(a)-(c); the more elongated the bond length is, the more delocalized the NWPs are, and vice versa. 44,60 This is a simple correlation between the bond length and WP width and does not mean that T2 molecules are more quantum compared to H2 and D2 molecules. In fact, the average WP width is larger in the order of H2 , D2 and T2 . (Figs.7(d)-(f)) The temperature dependence for the HWP and DWP width is similar to the temperature dependence of their average bond length shown in Fig.S7 where the bond length is simply correlated to the NWPs as explained above. However, T2 molecules exhibit more sensitive and qualitatively different temperature dependence; the average TWP width is rather anticorrelated to the T-T bond length. Such anomalous temperature dependence stems from the larger difference between the long and short T-T bond length discussed in Figs.6(a)-(c). T-T bond is relatively more elongated in the more condensed structure at the lower temperature compared to T-T bond in the less condensed structure at the higher temperature, and TWPs at edges of the stretched bond are relatively more influenced in the former by adjacent T2 molecules, that is, the TWPs are relatively localized due to a repulsive force from surrounding TWPs composing the stiffer solvation shell. This localization is more enhanced at the lower temperature, especially in the case of T2 molecules which are more structured in the SCL than H2 and D2 molecules, leading to the decrease of the average TWP width at the lower temperature. The frozen translational dynamics at the lower temperature also contributes to the current decrease by keeping TWPs more influenced by surrounding TWPs especially. In contrast, narrowing of the TWPs, which accompanies shrinkage of T-T bond outside a solvation shell, is determined rather by the intramolecular repulsion force between

20

ACS Paragon Plus Environment

Page 21 of 32

two TWPs of each molecule because each molecule outside a shell is freer and much less affected by surrounding TWPs of other molecules, giving little influence on the anomalous temperature dependence of the TWP width. The current split and broad distributions of the TWP width would not be affected even if we adopted a more complicated form of a NWP because each of the TWP width is rather small and each TWP keeps a Gaussian form very well, validating the present use of the Gaussian-type NWP. H2 D2 T2

4 Quantumness

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

3 2 1 5

10

15 20 Temperature (K)

25

Figure 8: Quantumness of H2 , D2 and T2 which measures anomalous quantum property of NWPs. The quantumness grows more rapidly in the order of H2 , D2 and T2 with decreasing the temperature. Figure 8 introduces a relative measure of quantumness defined by a ratio of peak height of a bond length distribution over peak height of a WP width distribution given in Figs.6(a)-(c) and Figs.7(a)-(c), respectively. 60 At moderate conditions, the NWP width is simply related to the bond length; as the bond length increases, the NWP becomes more delocalized, and vice versa. Evaluating relative deviation from this simple relation is important in order to examine the extent of anomalous nuclear delocalization which is not simply correlated to the bond length. This is achieved by comparing the relative quantity, the quantumness, defined by the ratio of the peak height of the bond distribution over the peak height of the NWP width distribution; the quantumness should be relatively constant as far as the nuclear delocalization is purely linked to the bond fluctuation. In this sense, the present quantumness

21

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

is different from other intuitive quantum measures such as NWP width itself which serves as a direct measure of nuclear delocalization. We note that the current quantumness was defined by dividing the original quantumness by the quantumness at 25 K of each isotope; such division enables us to equally compare the quantumness at once under the same condition which is free from the difference among the average WP width depending on the isotopes. The quantumness does not significantly change in the NL and in the WCL, while it is extremely enhanced in the SCL regardless of the isotopes. The rapid growth of the quantumness at 5 K can be explained by the two contrasting facts: On one hand, the NWP is more delocalized with the more elongated bond in the rigid mesoscale cluster. On the other hand, NWPs of molecules out of a cluster keep shrunk since the cluster formation is hindered at 5 K due to the little translational dynamics discussed in Fig.5(c). The two contrasting NWP behaviors cause the anomalously broad distributions of the NWP width in the SCL, leading to the rapid growth of the quantumness at 5 K. We found that the quantumness grows more rapidly in the order of H2 , D2 and T2 and the strongest quantumness appears in the H2 SCL in spite of its less irregular structure compared to the SCL structures of D2 and T2 . (See Figs.1(c) and 1(f)) The more anomalous HWP behavior not simply correlated to the H-H bond length should originate from its broadest NWP width among the current isotopes, and demonstrates that H2 molecules are more anomalously quantum than D2 and T2 molecules and thus have more possibility to achieve superfluidity. 2

Conclusions In summary, by comparing the real-time dynamics of H2 , D2 and T2 molecules including the intramolecular degrees of freedom in the various liquid phases, we computationally found that the significant isotope effects appear as the following differences in the order of H2 , D2 and T2 : (1) The liquid is less structured only at the lower temperature. (2) The molecular mobility is higher. (3) The librational frequencies are totally blueshifted. (4) The frequencies

22

ACS Paragon Plus Environment

Page 23 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of the low-energy collective dynamics are totally blueshifted. (5) The bond length is more elongated but the distribution of the bond length is sharper. (6) The vibrational bond frequencies and the average bond length are more sensitive to the temperature. (7) The WP width is more delocalized but the distribution of the WP width is sharper; especially the TWP distribution has the split double peaks and exhibits the qualitatively different temperature dependence. (8) The quantumness grows more rapidly with decreasing the temperature. It is remarkable that these isotopic difference among H2 , D2 and T2 molecules appears in spite of their same electronic potentials and the same thermodynamic conditions. All the isotopic difference was reasonably explained by the physical difference among H2 , D2 and T2 such as isotope-dependent solvation structure and dynamics including the mesoscale cluster in addition to the isotope-dependent mass, mobility, and nuclear delocalization. Furthermore, reproducing the experimental isotope dependence on the bond length and vibrational frequencies demonstrated that the NEWPMD method can properly describe the isotopic effects even in the condensed phases. We note that, although there is little doubt that the NEWPMD method gives the sensible values for the observables, the way to analyze the physical quantities using XYZ-coordinates of the NWP centers is not fully quantum mechanical and full quantum analyses in the NEWPMD method are an open question. We extracted and elucidated the pure isotopic effects appearing in the intramolecular and intermolecular dynamics of not only traditionally well studied isotopes, H2 and D2 , but also a lesser-known radioisotope, T2 . The predicted properties and obtained physical insights in this paper will definitely help future experimental searching and monitoring intermolecular and intramolecular dynamics and structures of these isotopes not only in normal liquid but also in supercooled liquid. Our work has initiated a quantitative trace of isotope-dependent molecular dynamics even in a condensed phase, and has wide scope for further extensions in various phases and thermodynamic states.

23

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 32

Supporting Information. Additional structural and dynamical data of hydrogen, deuterium and tritium liquids. This information is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgments KHD thanks the financial supports from JST (PRESTO), Japan Association for Chemical Innovation (JACI), Cybermedia Center of Osaka University, and Grant-in-Aids for Scientific Research from Japan Society for the Promotion of Science (KAKENHI), Grant No. 15K05386.

References (1) Bermejo, F. J.; Kinugawa, K.; Dawidowski, J.; Cabrillo, C.;Fernandez-Perea, R. Beyond Classical Molecular Dynamics: Simulation of Quantum-Dynamics Effects at Finite Temperatures; the Case of Condensed Molecular Hydrogen. Chem. Phys. 2005, 317,198–207. (2) Kuhnel, M.; Fernandez, J. M.; Tramonto, F.; Tejeda, G.; Moreno, E.; Kalinin, A.; Nava, M.; Galli, D. E.; Montero, S.;Grisenti, R. E. Observation of Crystallization Slowdown in Supercooled Parahydrogen and Orthodeuterium Quantum Liquid Mixtures. Phys. Rev. B. 2014, 89,180201. (3) Mompean, F. J.; Garcia-Hernandez, M.; Bermejo, F. J.;Bennington, S. M. SingleMolecule Kinetic Energy of Condensed Normal Deuterium. Phys. Rev. B. 1996, 54,970–977. (4) Walker, J. A.; Mezyk, S. P.; Roduner, E.;Bartels, D. M. Isotope Dependence and Quantum Effects on Atomic Hydrogen Diffusion in Liquid Water. J. Phys. Chem. B. 2016, 120,1771–1779. 24

ACS Paragon Plus Environment

Page 25 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(5) Bermejo, F. J.; Fak, B.; Bennington, S. M.; Kinugawa, K.; Dawidowski, J.; FernandezDiaz, M. T.; Cabrillo, C.;Fernandez-Perea, R. Microscopic Structure Factor of Liquid Parahydrogen: Monte Carlo and Molecular Dynamics Simulations. Phys. Rev. B. 2002, 66,212202. (6) Bermejo, F. J.; Kinugawa, K.; Cabrillo, C.; Bennington, S. M.; Fak, B.; FernandezDiaz, M. T.; Verkerk, P.; Dawidowski, J.;Fernandez-Perea, R. Quantum Effects on Liquid Dynamics as Evidenced by the Presence of Well-Defined Collective Excitations in Liquid para-Hydrogen. Phys. Rev. Lett. 2000, 84,5359–5362. (7) Galliero, G.; Boned, C. Thermal Conductivity of the Lennard-Jones Chain Fluid Model. Phys. Rev. E. 2009, 80,061202. (8) Scharf, D.; Martyna, G. J.;Klein, M. L. Structure and Energetics of Fluid paraHydrogen. Low Temp. Phys. 1993, 19,364–367. (9) Yonetani, Y.; Kinugawa, K. Transport Properties of Liquid para-Hydrogen: The Path Integral Centroid Molecular Dynamics Approach. J. Chem. Phys. 2003, 119,9651– 9660. (10) Yonetani, Y.; Kinugawa, K. para-hydrogen ni CMD yuukou. J. Chem. Phys. 2004, 120,10624–10633. (11) Zoppi, M.; Neumann, M.;Celli, M. Microscopic Structure Factor of Liquid paraHydrogen. Phys. Rev. B. 2002, 65,092204. (12) Souers, P. C. Hydrogen Properties for Fusion Energy. ; University of California Press : Berkeley, USA, 1986. (13) Zoppi, M.; Celli, M.;Soper, A. K. Neutron Diffraction Determination of the Thermodynamic Derivatives of the Microscopic Structure of Liquid Parahydrogen. Phys. Rev. B. 1998, 58,11905–11910. 25

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 32

(14) Bermejo, F. J.; Mompean, F. J.; Garcia-Hernandez, M.; Martinez, J. L.; MartinMarero, D.; Chahid, A.; Senger, G.;Ristig, M. L. Collective Excitations in Liquid Deuterium: Neutron-Scattering and Correlated-Density-Matrix Results. Phys. Rev. B. 1993, 47,15097–15112. (15) Celli, M.; Bafile, U.; Cuello, G. J.; Formisano, F.; Guarini, E.; Magli, R.; Neumann, M.;Zoppi, M. Microscopic Structure Factor of Liquid Hydrogen by Neutron-Diffraction Measurements. Phys. Rev. B. 2005, 71,014205. (16) Celli, M.; Magli, R.; Fischer, H.; Frommhold, L.;Barocchi, F. Quantum Mechanical Effects on the Static Structure Factor of Pairs of Orthodeuterium Molecules. Phys. Rev. Lett. 1998, 81,5828–5831. (17) Mukherjee, M.; Bermejo, F. J.; Fak, B.;Bennington, S. M. Microscopic Dynamics in Liquid Deuterium: A Transition from Collective to Single-Particle Regimes. Europhys. Lett. 1997, 40,153–158. (18) O‘Reilly, D. E.; Peterson, E. M. Self-Diffusion of Liquid Hydrogen and Deuterium. J. Chem. Phys. 1977, 66,934–937. (19) Zoppi, M. Microscopic Structure of Liquid Hydrogen. J. Phys. Cond. Mett. 2003, 15,R1047. (20) Zoppi, M.; Bafile, U.; Magli, R.;Soper, A. K. Neutron-Diffraction Determination of the Microscopic Structure of Liquid Deuterium. Phys. Rev. E. 1993, 48,1000–1007. (21) Zoppi, M.; Soper, A. K.; Magli, R.; Barocchi, F.; Bafile, U.;Ashcroft, N. W. Structure Factor of Compressed Liquid Deuterium Close to the Melting Transition. Phys. Rev. E. 1996, 54,2773–2779. (22) Kosmider, A.; Drexlin, G.; Eichelhardt, F.; Michling, R.; Welte, S.;Wurster, W. In-

26

ACS Paragon Plus Environment

Page 27 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

frared Spectroscopy in Liquid Hydrogen Isotopologues for the ISS of ITER. Fus. Sci. Tech. 2011, 60,956–959. (23) Georgescu, I.; Deckman, J.; Fredrickson, L. J.;Mandelshtam, V. A. Thermal Gaussian Molecular Dynamics for Quantum Dynamics Simulations of Many-Body Systems: Application to Liquid para-Hydrogen. J. Chem. Phys. 2011, 134,174109. (24) Herrero, C. P.; Ramirez, R. Molecular Hydrogen in Graphite: A Path-Integral Simulation. Phys. Rev. B. 2010, 82,174117. (25) Hone, T. D.; Pulsen, J. A.; Rossky, P. J.;Manolopoulos, D. E. Comparison of Approximate Quantum Simulation Methods Applied to Normal Liquid Helium at 4 K. J. Phys. Chem. B. 2008, 112,294–300. (26) Liu, J.; Miller, W. H. Test of the Consistency of Various Linearized Semiclassical Initial Value Time Correlation Functions in Application to Inelastic Neutron Scattering from Liquid para-Hydrogen. J. Chem. Phys. 2008, 128,144511. (27) Miller, T. F.; Manolopoulos, D. E. Quantum Diffusion in Liquid para-Hydrogen from Ring-Polymer Molecular Dynamics. J. Chem. Phys. 2005, 122,184503. (28) Nakayama, A.; Makri, N. ForwardBackward Semiclassical Dynamics for Quantum Fluids using Pair Propagators: Application to Liquid para-Hydrogen. J. Chem. Phys. 2003, 119,8592–8605. (29) Poulsen, J. A.; Nyman, G.;Rossky, P. J. Quantum Diffusion in Liquid Para-hydrogen: An Application of the Feynman-Kleinert Linearized Path Integral Approximation. J. Phys. Chem. B. 2004, 108,19799. (30) Smith, K. K. G.; Poulsen, J. A.; Cunsolo, A.;Rossky, P. J. Refinement of the Experimental Dynamic Structure Factor for Liquid para-Hydrogen and ortho-Deuterium using Semi-Classical Quantum Simulation. J. Chem. Phys. 2014, 140,034501. 27

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 32

(31) Smith, K. K. G.; Poulsen, J. A.; Nyman, G.; Cunsolo, A.;Rossky, P. J. Application of a New Ensemble Conserving Quantum Dynamics Simulation Algorithm to Liquid para-Hydrogen and ortho-Deuterium. J. Chem. Phys. 2015, 142,244113. (32) Zoppi, M.; Bafile, U.; Guarini, E.; Barocchi, F.; Magli, R.;Neumann, M. Microscopic Structure and Intermolecular Potential in Liquid Deuterium. Phys. Rev. Lett. 1995, 75,1779–1782. (33) Colognesi, D.; Bafile, U.; Celli, M.; Neumann, M.;Orecchini, A. Hydrogen SelfDynamics in Liquid H2 -D2 Mixtures Studied Through Inelastic Neutron Scattering. Phys. Rev. E. 2015, 92,012311. (34) Geneste, G.; Torrent, M.; Bottin, F.;Loubeyre, P. Strong Isotope Effect in Phase II of Dense Solid Hydrogen and Deuterium. Phys. Rev. Lett. 2012, 109,155303. (35) Kuhnel, M.; Fernandez, J. M.; Tramonto, F.; Tejeda, G.; Moreno, E.; Kalinin, A.; Nava, M.; Galli, D. E.; Montero, S.;Grisenti, R. E. Mixing Effects in the Crystallization of Supercooled Quantum Binary Liquids. J. Chem. Phys. 2015, 143,064504. (36) Mallory, J. D.; Mandelshtam, V. A. Quantum Melting and Isotope Effects from Diffusion Monte Carlo Studies of p-H2 Clusters. J. Phys. Chem. A. 2017, 121,6341– 6348. (37) Garashchuk, S.; Jakowski, J.; Wang, L.;Sumpter, B. G.

Quantum Trajectory-

Electronic Structure Approach for Exploring Nuclear Effects in the Dynamics of Nanomaterials. J. Comp. Theor. Chem. 2013, 9,5221–5235. (38) Hyeon-Deuk, K. ; Ando, K. Intermolecular Diatomic Energies of a Hydrogen Dimer with non-Born-Oppenheimer Nuclear and Electron Wave Packets. Chem. Phys. Lett. 2012, 532,124–130.

28

ACS Paragon Plus Environment

Page 29 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(39) Hyeon-Deuk, K. ; Ando, K.

Quantum Molecular Dynamics Simulation of Liq-

uid para-Hydrogen by Nuclear and Electron Wave Packet Approach.

J. Chem.

Phys.(Communication). 2014, 140,171101. (40) Burghardt, I.; Giri, K. ; Worth, G. A. Multimode Quantum Dynamics Using Gaussian Wavepackets: The Gaussian-Based Multiconfiguration Time-Dependent Hartree Method Applied to the Absorption Spectrum of Pyrazine J. Chem. Phys. 2008, 129,174104. (41) Heatwole, E.; Prezhdo, O. V. Canonical Averaging in the Second Order Quantized Hamilton Dynamics by Extension of the Coherent State Thermodynamics of the Harmonic Oscillator J. Chem. Phys. 2007, 126,204108. (42) Hyeon-Deuk, K.; Ando, K. Semiquantum Molecular Dynamics Simulation of Liquid Water by Time-Dependent Hartree Approach J. Chem. Phys. 2009, 131,064501. (43) Hyeon-Deuk, K.; Ando, K. Quantum Effects of Hydrogen Atoms on the Dynamical Rearrangement of Hydrogen-Bond Networks in Liquid Water J. Chem. Phys. 2010, 132,164507. (44) Hyeon-Deuk, K. ; Ando, K. Correlations of Intra- and Intermolecular Dynamics and Structure in Liquid para-Hydrogen. Phys. Rev. B. 2014, 90,165132. (45) Hyeon-Deuk, K. ; Ando, K. Dynamical and Structural Analyses of Solid Hydrogen under Vapor Pressure. J. Chem. Phys.(Communication). 2015, 140,171102. (46) Abe, K.; Hyeon-Deuk, K. Dynamical Ordering of Hydrogen Molecules Induced by Heat Flux. J. Phys. Chem. Lett. 2017, 8,3595–3600. (47) Silvera, I. F.; Goldman, V. V. The Isotropic Intermolecular Potential for H2 and D2 in the Solid and Gas Phases. J. Chem. Phys. 1978, 69,4209–4213.

29

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 32

(48) Grebenev, S.; Sartakov, B.; Toennies, J. P.;Vilesov, A. F. Evidence for Superfluidity in Para-Hydrogen Clusters Inside Helium-4 Droplets at 0.15 Kelvin. Science. 2000, 289,1532–1535. (49) Khairallah, S. A.; Sevryuk, M. B.; Ceperley, D. M.;Toennies, J. P. Interplay between Magic Number Stabilities and Superfluidity of Small Parahydrogen Clusters. Phys. Rev. Lett. 2007, 98,183401. (50) K¨ uhnel, M.; Fern´andez, J. M.; Tejeda, G.; Kalinin, A.; Montero, S.;Grisenti, R. E. Time-Resolved Study of Crystallization in Deeply Cooled Liquid Parahydrogen. Phys. Rev. Lett. 2011, 106,245301. (51) Kuyanov-Prozument, K.; Vilesov, A. F. Hydrogen Clusters that Remain Fluid at Low Temperature. Phys. Rev. Lett. 2008, 101,205301. (52) Kwon, Y.; Whaley, K. B. Nanoscale Molecular Superfluidity of Hydrogen. Phys. Rev. Lett. 2002, 89,273401. (53) Li, H.; Le Roy, R. J.; Roy, P.-N.;McKellar, A. R. W. Molecular Superfluid: Nonclassical Rotations in Doped Para-Hydrogen Clusters. Phys. Rev. Lett. 2010, 105,133401. (54) Sindzingre, P.; Ceperley, D. M.;Klein, M. L.

Superfluidity in Clusters of p-H2

Molecules. Phys. Rev. Lett. 1991, 67,1871–1874. (55) Sliter, R.; Vilesov, A. F. Temperature Dependence of the Raman Spectra of Liquid Parahydrogen. J. Chem. Phys. 2009, 131,074502. (56) Zeng, T.; Roy, P.-N. Microscopic Molecular Superfluid Response: Theory and Simulations. Rep. Prog. Phys. 2014, 77,046601. (57) Osychenko, O. N.; Rota, R.;Boronat, J. Superfluidity of Metastable Glassy Bulk para-Hydrogen at Low Temperature. Phys. Rev. B. 2012, 85,224513.

30

ACS Paragon Plus Environment

Page 31 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(58) Helfferich, J.; Ziebert, F.; Frey, S.; Meyer, H.; Farago, J.; Blumen, A.;Baschnagel, J. Continuous-Time Random-Walk Approach to Supercooled Liquids. I. Different Definitions of Particle Jumps and Their Consequences. Phys. Rev. E. 2014, 89,042603. (59) Helfferich, J.; Ziebert, F.; Frey, S.; Meyer, H.; Farago, J.; Blumen, A.;Baschnagel, J. Continuous-Time Random-Walk Approach to Supercooled Liquids. II. Mean-Square Displacements in Polymer Melts. Phys. Rev. E. 2014, 89,042604. (60) Hyeon-Deuk, K. ; Ando, K. Distinct Structural and Dynamical Difference between Supercooled and Normal Liquids of Hydrogen Molecules. Phys. Chem. Chem. Phys. 2016, 18,2314–2318. (61) Hone, T. D.; Voth, G. A. A Centroid Molecular Dynamics Study of Liquid paraHydrogen and ortho-Deuterium. J. Chem. Phys. 2004, 121,6412–6422. (62) Miller, T. F.; Manolopoulos, D. E.; Madden, P. A.; Konieczny, M.;Oberhofer, H. Quantum Diffusion in Liquid para-Hydrogen from Ring-Polymer Molecular Dynamics. J. Chem. Phys. 2005, 122,057101. (63) Yamamuro, O.; Tsukushi, I.; Kanaya, T.;Matsuo, T. Neutron Scattering Study of Boson Peak in Glassy 3-Methylpentane. J. Phys. Chem. Sol. 1999, 60,1537–1539. (64) Loubeyre, P.; Jean-Louis, M.;Silvera, I. F. Density Dependence of the Intramolecular Distance in Solid H2 . A. Spectroscopic Determination. Phys. Rev. B. 1991, 43,10191– 10196.

31

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic.

32

ACS Paragon Plus Environment

Page 32 of 32