Ice-Amorphization of Supercooled Water ... - ACS Publications

Article http://pubs.acs.org/journal/aesccq

Ice-Amorphization of Supercooled Water Nanodroplets in No Man’s Land Prithwish K. Nandi,* Christian J. Burnham, Zdenek Futera,† and Niall J. English* School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland S Supporting Information *

ABSTRACT: Elucidating freezing mechanisms of liquid water into ice, especially in “No Man’s Land” (150 K< T < 235 K), carries scientific and technological importance. Indeed, superior predictions of upper-troposphere cirrus-cloud formation and surface-bound ice-fog formation constitute powerful motivations in addition to unravelling long-standing puzzles such as persistent liquid fogs well below frost point and understanding interstellar-space water states, together with advancing cryopreservation technology. Unlocking the secrets of water’s anomalous deep-cooling complexities, such as structural ordering and microscopic nucleation mechanisms, are the subject of lively debate. Exploring nucleation mechanism in No Man’s Land (NML) is technically demanding, owing to rapid nucleation rates with, unsurprisingly, very few reported experimental studies. However, amorphization is a key intermediate stage in NML-based nucleation, and it is also not particularly well understood. In this microsecond long molecular dynamics study, we have explored microstructural processes involved in the amorphization of aggressively quenched supercooled water nanodroplets in the gas phase where surface effects are non-negligible. A dynamically arrested state is observed in these droplets that resembles structurally low-density amorphous polymorphs of ice. Importantly, the curved geometry of the nanodroplets themselves is seen to inhibit amorphization relative to bulk systems under identical thermodynamic conditions. KEYWORDS: Amorphization, ice, No Man’s Land, low-density amorphous ice, supercooled water

1. INTRODUCTION Crystallization of water into ice is a long-standing research conundrum. Given water’s fundamental importance for life and nature, extensive efforts have been made to understand its chemistry, physical properties, thermodynamic phases, and the morphology thereof. However, there exists still many open questions. One such less-explored area is water’s “No Man’s Land” (NML) that typically ranges from ∼150 K, below which amorphous water can be detected, to ∼235 K, above which homogeneous nucleation of ice can be seen.1 In No Man’s Land, the nucleation rate of supercooled water to ice is so rapid that it is hardly possible to probe this regime satisfactorily through experiments.1,2 Understanding the underlying freezing process of supercooled water to ice in NML could prove to be useful to elucidate and even predict many puzzling environmental phenomena like cloud formation in the atmosphere,3−5 understanding reasons behind persistent liquid fogs well below the frost point,6 the formation of ice fogs either on earth’s surface7 or in cirrus clouds,8−10 together with states of water in interstellar space.11,12 The nature of microstructural changes and the associated mechanistic processes underlying arrangement/reorientation of water molecules in the supercooled state to form either crystalline or amorphous configurations is widely debated. The description of very fast ice crystallization is based on classical © XXXX American Chemical Society

volume-based homogeneous nucleation rate theory proposed by Volmer,13 where the nucleus was assumed to be formed mainly within the droplet volume. The theory remained almost unchanged until 2002 when Tabazadeh et al.6 showed that such nucleation could also be a surface-based process. The first attempt to simulate NML freezing of water was reported by Moore and Molinero.2 Using a coarse-grained monatomic water (mW) model and a bulk water system, they identified three distinct stages of ice crystallization at ∼180 K. The first stage (Stage I) is characterized by the random formation of ice nuclei throughout the volume; the second stage (Stage II) is dominated by the growth and consolidation of nuclei; in the third stage (Stage III), large-enough ice crystallites start fusing with each other. Although ice nucleation starts with Stage I, the radial distribution function (RDF) at this point does not exhibit a crystal-signature; rather it resembles that of low-density amorphous ice (LDA). In fact, LDA has often been regarded as a precursor for ice nucleation.14,15 Another important study on supercooled-water crystallization was reported by Malkin et al.,16 who showed via X-ray diffractometry that micron-sized Received: February 10, 2017 Revised: April 10, 2017 Accepted: May 8, 2017

A

DOI: 10.1021/acsearthspacechem.7b00011 ACS Earth Space Chem. XXXX, XXX, XXX−XXX

Article

ACS Earth and Space Chemistry

microseconds-long simulations that for such nanodroplets, the molecular systems enter into dynamically arrested states that we confirmed to resemble those of low-density amorphous ice polymorphs. Moreover, we also show that the confined geometry of the droplet acts as an inhibitor for the amorphization process itself vis-à-vis bulk systems under identical conditions.

water droplets crystallize homogeneously from the supercooled state; the adopted structure was composed of randomly stacked layers of hexagonal and cubic sequences, known as stackingdisordered ice I. Such an ice form was reported previously for water in confined geometries, such as pores or capillaries,17,18 and in samples recrystallized from high-pressure ice phases, as opposed to in water droplets at ambient conditions. Indeed, such confined nanosized pores/capillaries were chosen mainly to reduce the probability of fast nucleation. A great deal of experimental efforts in this vein have been reported.19−33 In particular, a recent and interesting study of Manka et al.34 reports nanometer-scale droplets (radii ranging from 3.2 to 5.8 nm) in a supersonic nozzle. This led to a homogeneous ice nucleation rate of the order of 1023 cm−3 s−1 at temperatures in between 202 and 215 K (inferred by way of pressure-trace and small-angle X-ray scattering measurements and Fourier transform infrared spectroscopy). In addition to various (quasi-) crystalline forms (e.g., hexagonal, cubic, stacking-disordered), ice also exists in amorphous forms, primarily high- and low- density amorphous (HDA and LDA) states.1,35−37 Enhancing our understanding of amorphous ice states is crucial, as it often helps to explain the anomalies of supercooled water in terms of a critical point above which the two metastable amorphous phases of ice (i.e., LDA and HDA) become indistinguishable. The phase behavior described by Poole et al.38 has several other implications that include understanding of the development of random tetrahedral network upon supercooling and even explaining properties of water in confined geometries where bulk liquidlike behavior could be expected at comparatively higher temperatures. Owing to its deep implications in understanding supercooling anomalies, there have been many reports in the literature of microscopic studies of amorphous ice too, including many important contributions like the ice-polyamorphism studies at low temperature from constant-pressure simulations,39 reversible pressure-induced LDA−HDA35,36 and HDA−Ih transformations,33−37 and pressure-induced LDA− HDA−VHDA transformations, where VHDA stands for very high-density amorphous ice.40 Thus, exploration of thermodynamic and kinetic behavior of supercooled water in No Man’s Land is important, and most previous computational studies involve bulk systems. In a recent molecular dynamics study by Li et al.,41 the dependence of ice nucleation rate on the size of the nanodroplet was calculated in detail in No Man’s Land and a thermodynamic model was also proposed to describe this size dependence. However, in the present work, using molecular dynamics simulations, we aim to understand underlying ice-amorphization mechanisms (as opposed to outright nucleation per se) in No Man’s Land in a “confined” (spherical) geometry of nanoscale water droplets in vacuo at atomistic level, where surface effects (due to Laplace pressure) are expected to be non-negligible. In part, the experimental relevance of the aerosol-based experiments of Manka et al.34 lies in the efficient controlling of the fast nucleation process generally observed in the NML region (i.e., retardation to the extent of practical experimental observability, in contrast to bulk systems); this acts as further motivation to study nanodroplets. However, in the present work we stress that we have little interest in the results of the work of Manka et al.34 per se; rather, we wish to mimic the aerosol-based premise of the work of Manka et al.,34 given its experimental relevance, and also to compare droplet versus bulk-system amorphization behavior. We show from our

2. METHODS All MD simulations were performed using GROMACS42 under periodic boundary condition with TIP4P/2005 potential model for water.43 The system consists of a spherical water droplet of radius ∼21.5 Å consisting about 1400 water molecules. The droplet was placed at the center of a cubic periodic box of size ∼245 Å so that the interactions between its periodic images are essentially nonexistent, mimicking the vacuum state of ref 34. Each simulation run was carried out at constant temperature and constant volume. A Nosé−Hoover44 (NH) thermostat was applied to impose the desired system temperature with coupling constant of 1 ps. Sudden cooling at the beginning of the simulation exhibits typical oscillatory behavior of system temperature during the first 7−8 ps of the simulation (cf. Figure S1, Supporting Information). After this period, the system temperature oscillated around its prescribed value. We are aware that the initial oscillations are largely affected by the NH thermostat, however, the NH algorithm was chosen because it can reproduce canonical ensemble and it is time reversible (in contrast to velocity scaling, weak coupling or stochastic approaches). The system was propagated in time using the leapfrog MD integrator with time step of 3 fs. The relatively long time step was justified by constraining all bond lengths and angles in the system by a LINCS algorithm of fourth order.45 As there is minimal movement of molecules in the simulations, we updated the Verlet neighbor list every 10 steps, which we found sufficient by checking NVE dynamics, where no jumps in potential/total energies were observed, coupled with good energy conservation. Although our system is aperiodic, we decided to perform the simulations in a large cubic periodic box with a side length of 245 Å and evaluated the electrostatic interaction by smooth particle-mesh Ewald (SPME).46 This decision was made on the basis of our benchmark runs when we observed better computational efficiency than for direct real-space sum of Coulombic contributions. The real-space cutoff in PME was set to 8 Å and 8 × 8 × 8 mesh was used to evaluate the long-range contribution in reciprocal space. A relatively small number of mesh points were justified by a large simulation cell and reasonably good conservation of total energy during NVE benchmarks. A systematic study was done for four different temperatures: 150, 175, 200, and 225 K. The initial configurations were taken from a single MD trajectrory of water at 300 K, and the configurations were quenched instantaneously to the target temperature. Each simulation was run for a duration of at least 1 μs, where trajectories were sampled every 10 ps for analysis. Oxygen−oxygen radial distribution function for these trajectoV ries was computed as g (r ) = 4πr 2N2 ⟨∑i ∑j ≠ i δ(r − rij)⟩, where δ is the Dirac’s delta function, and rij is the separation vector between the oxygen atoms of the ith and jth water molecules in the system. The underlying ring structures were analyzed following the methodology outlined in refs 47−51. Hydrogen bonding connectivity between water molecules was determined B


Article


Figure 1. Snapshot of a thin slice of width of 4 Å at the middle of the water nanodroplet at (a) the beginning of the simulation and (b) after 1 μs at the end of the simulation. The final configuration in (b) shows clearly the ordered ringlike structures. A few of the ringlike structures are shown in (b) schematically by connecting the oxygen sites (red balls) with black dotted lines.

following the Luzar−Chandler criteria:52 (i) interoxygen separation is less than 3.5 Å and (ii) the angle between the O−O axis and the one of the O−H bond is less than 30°. The n-membered ring was then identified by finding n distinct water molecules that are connected sequentially by means of hydrogen bonds. If two rings are “short-circuited” by a common side, we counted both rings, and the large encompassing ring, which one could imagine by disregarding the common side, was not taken into count. The mean square displacement of the center of mass of the water molecules was calculated using the standard multiple-time origin method.53 The average tetrahedral-order parameter (q-parameter) associated with the water molecules of the droplet under freezing condition was calculated using the following relation:54 q=1−

3 8

3

4

(

∑ j = 1 ∑k = j + 1 cos ψjk + N

tion of q is calculated as

1 3

which shows that single- and double- precision GROMACS versions produce almost identical results. In order to understand the effect of curvature of water molecules in a volume occupied by spherical nanodroplets, we also carried out the simulations for bulk water at four different temperatures used here with 1 μs long simulation times. For the bulk simulations, we used the NPT ensemble (with a target pressure of 1 bar and featuring the same Nosé−Hoover thermostat setup) while all other parameters remained identical.

3. RESULTS AND DISCUSSION The subtle differences between liquid water and its various crystalline/amorphous forms are mainly manifested by changes in its local orientational order. Although water molecules in a hexagonal ice (Ih) crystal are oriented in perfect tetrahedral order, in other forms (e.g., amorphous ice) such perfect tetrahedrality can never be realized. The initial perception of change of orientational order can often be well recognized by visualization of structure. In Figure 1, we have shown the snapshots of the cross sections of the nanodroplet configurations at the very outset (cf. Figure 1a) and end of the 1 μs simulation (cf. Figure 1b) at 150 K. The formation of various ring structures is clearly visible when compared to the initial liquid nanodroplet of water (see below). This observation indicates clearly that at 150 K water molecules undergo substantial structural ordering vis-à-vis its liquid counterpart at 300 K. A quantitative measure of spatial ordering is given by the RDF. In Figure 2a, RDF plots are shown for the core region (r < 15 Å) of the final nanodroplet configuration obtained after 1 μs run at the four different temperatures studied along with the initial liquid-like droplet at 300 K. In Figure 2b,c, the temporal evolution of the RDF plots are shown for the core (r < 15 Å) and surface regions (15 Å < r < 22 Å) of the nanodroplet at 150 K. As the temperature is lowered from 300 K (i.e., the liquidlike case), the intensity of the first, second, and third RDF peaks increases considerably, signifying more local structural ordering than the liquid (cf. Figure 2a). At 150 K, the intensity of the first minimum touches zero, resembling crystal-like behavior. A slight shift of the second peak occurs as the temperature is decreased from 300 K toward 150 K. Noticing that the maximum change occurred at 150 K relative to liquid water (at 300 K), we analyzed the droplet at 150 K further by dividing it into two zones (see above). Although the position of

2

) . The standard deviaN

N ∑i = 1 qi2 − (∑i = 1 qi)2 N (N − 1)

. The q parameter,

although useful for preliminary assessment of enhanced tetrahedrality, is rather insensitive to distinguish among various forms of crystalline, amorphous ice, and liquid water. Numerous approaches exist, for example, Steinhardt order parameters,55 CHILL algorithm,56 the Báez−Clancy approach57 or the method outlined by Reinhardt et al.58 In this work, we have used the CHILL algorithm that is proven to be successful to distinguish between cubic ice, hexagonal ice, LDA, and liquid water. In this method, very similar to that of Reinhardt et al.,58 the local order around each water molecule i is described by a local orientational bond order parameter vector q̇l(i) having 2l + 1 components defined as (eq 4 of ref 56): 1 4 qlm(i) = 4 ∑ j = 1 Ylm( rij)̂ . The distribution of the alignment of local orientational structure, as defined by eq 5 of ref 56, a ( i , j) =

ql(̂ i) . ql(̂ j) |ql(̂ i)| . |ql(̂ j)|

, for l = 3 and l = 4 can distinguish between

various polymorphs of ice. It is important to mention here that by default GROMACS is compiled in single precision and uses default update of neighbor lists and approximations to speed up electrostatics. A previous study by Wong-Ekkabut et al.59 showed that this might potentially lead to less physically reasonable results. In order to check the accuracy of our simulations, we have also repeated them using a double-precision GROMACS version with a shorter 2 fs time step. A comparison of results is provided in the Supporting Information (cf. Figure S2a−c), C


Article


except for distinctive ice Ih third and fourth peaks at ∼5.2764 and ∼6.27 Å, respectively. Earlier MD studies,2 using the mW water model65 of Moore and Molinero, showed that in NML ice crystallization the third peak does not evolve until the process enters Stage II, as described in ref 2 (and also mentioned in the Introduction), wherein ice crystallites start thickening. In the structural characterization of ice/water, the significance of the position and intensity of the second peak in the O−O RDF at ∼4.5 Å is well-known, as this corresponds to the second nearest-neighbor position in the tetrahedral coordination, 8/3 a where a ∼ 2.8 Å.63,66 Therefore, here the increased relative intensity of the second RDF peak and the shift toward the value of 4.5 Å with lowering temperature essentially points to the increment of the degree of the tetrahedrality supported by directional hydrogen bonds adopted by water molecules. The comparison of the O−O pair distribution function between the nanodroplet and bulk water is presented in Figure S3 (Supporting Information), which exhibits no significant difference between the peak positions. In MD analysis, tetrahedrality is often measured by the qparameter (cf. Methods), the average value of which ranges from 0 (ideal gas) to 1 (perfect tetrahedral solid). The evolution of q with time for four different temperatures is shown in Figure 3a. As temperature is lowered, the value of q approaches toward 1, and fluctuations in the underlying value also diminishes considerably. As the temperature gets lower, q also requires a longer time to converge, unsurprisingly. The converged q value for the final configuration is shown in Figure 3b at the four different temperatures studied here. The corresponding values of q parameter for the bulk case at the end of 1 μs is also shown in Figure 3b. For the nanodroplet, q at 225 K is ∼0.76, while for bulk water it is ∼0.79, agreeing well with the value (∼0.8) reported by Wong et al.35 at ∼225 K. As temperature is lowered from 225 to 175 K, the value of q increases, indicating enhanced tetrahedrality. Interestingly, for all these three temperatures (175, 200, and 225 K), bulk water shows more tetrahedral order than for the water molecules in the nanodroplet. This is also indicative of inhibition of orientational ordering in the nanodroplet as compared to the bulk case. This inhibition is probably a consequence of surfaceinduced Laplace pressure onto the core of the droplet (will be discussed in the following section in more detail). As temperature is lowered from 175 K toward 150 K, q saturates to ∼0.9 for the nanodroplet and for the bulk the q-value decreases slightly. It is worth mentioning that the saturated values for both the nanodroplet as well as for the bulk at 150 K agree well with the reported value of q for LDA ice (∼0.89) at 80 K with the same water model.35 LDA ice with other water models also leads to identical/comparable q values (0.9 with ST2 at 80 K and 0.82−0.84 with SPC/E at 77 K and −500 < P < 400 MPa).35 The decrease in q’s rate of change from 175 to 150 K is consistent with the observation of Wong et al.35 who showed that q barely changes upon cooling from 200 K to temperatures less than 140 K. The convergence of q at ∼0.9 points toward arrested reorientational dynamics of water molecules that prevent the system reaching a perfect tetrahedral order corresponding to ice Ih (which is also the stable form of ice for TIP4P/200567) within microsecond time scales. The signature of a system with arrested glassy dynamics is primarily manifested in its mean square displacement (MSD), a plot of which is shown in Figure 4 for six different temperatures ranging from 100 to 225 K. The plateau regions in the MSD for

Figure 2. O−O radial distribution function (a) at four different temperatures, (b) for the core region of the nanodroplet (r < 15 Å) at 150 K, and (c) for the surface layer of the droplet at 150 K. As temperature is lowered, the second peak shifts toward a value of 4.5 Å, showing an increment in tetrahedral order. However, the RDF agrees quite well with experimental RDF data for LDA ice.

the first peak remains unaltered, progressing radially outward, the shift of the second-peak position is quite pronounced. Indeed, there is strong similarity in RDF plots, especially at 150 K, with those for LDA ice: X-ray crystallographic data of the O−O RDF exhibits peaks at 2.79, 4.56, and 8.60 Å,60−62 comparable with respective positions here (i.e., 2.73, 4.45, and 6.68 Å in the core and 2.73, 4.39, and 6.67 Å in the surface region). Moreover, there is rather striking similarity of these RDF plots with TIP4P/2005 results in ref 63 (featuring nonequilibrium MD mimicking of evaporative cooling to support experimental study). Certainly, the first two peak positions of the O−O RDF plot for ice Ih are also comparable, D


Article


One of the main anomalies exhibited by water upon cooling toward its freezing point is the decrease of its density. The various polymorphs of ice show different characteristic equilibrium density at ambient condition. In Figure 5a, we

Figure 3. Orientational order parameter q (a) as a function of time at four different temperatures and (b) as a function of temperature for the final nanodroplet configuration obtained at the end of 1 μs simulation run. Clearly, as temperature is lowered, q increases, indicating enhanced tetrahedral order. However, the value saturates at ∼0.9 for the lowest studied temperature (i.e., 150 K). This value of q agrees quite well with other reported LDA values (see main text). The q-values for the bulk water is also plotted in (b), which shows enhanced tetrahedrality as compared to nanodroplets under identical thermodynamic conditions.

Figure 5. Radial density-profile plots (a) for four different temperatures at the end of 1 μs simulationand (b) at 150 K at different stages of the simulation. The density profile shows a very distinct density layering, noticeably at 150 K, and for 175 and 200 K, albeit less prominent. The core density at 150 K approaches a value of ∼0.94 that agrees quite well with experimental LDA density. The density for bulk water at identical thermodynamic conditions is also shown by dotted horizontal lines in (a). In (b), we can observe core densities to be higher than density at ambient conditions (1 g·cm−3) that is associated with the Laplace pressure induced by the surface onto the core of the nanodroplet.

have shown the nanodroplet’s radial density variation outward for four different temperatures for the final configuration after 1 μs. The corresponding densities for bulk water at these temperatures are also shown in this plot by the horizontal dotted line. Although the bulk densities (∼0.94 g·cm−3) for 150, 175, and 200 K are almost identical, resembling the density of LDA ice,35,68 for 225 K the density is slightly higher (∼0.96 g·cm−3). Interestingly, for the nanodroplet, this density profile shows very distinct “layering”, noticeably at 150 K, and for 175 and 200 K, although less prominent (cf. Figure 5a). At 225 K, the system can be regarded as liquid-like as the density matches with the equilibrium liquid density at ambient condition (1.0 g·cm−3) almost up to a distance rc = 0.93Rinitial, where Rinitial is the radius of the nanodroplet at the start of the simulation. This specific density layering is very much consistent with the density layering observed for the case of crystallization in nanoscale Si droplets.69 On the basis of this

Figure 4. Mean square displacement plot for the oxygen atoms at six different temperatures. Although the plot shows a glassy nature for 100, 125, and 150 K, for 175, 200, and 225 K, its shows a subdiffusive behavior. There is an onset of subdiffusive behavior observed at 175 K.

100, 125, and 150 K indicate “trapped” dynamics redolent of glasses, whereas at 175 K the system exhibits the onset of subdiffusive dynamics, which serves to shorten the characteristic cage-disintegration time of diffusion in any liquid. E


Article

ACS Earth and Space Chemistry density layering, the droplet can be broadly divided into bulk and surface parts. In Figure 5b, we have shown how this density layering and the subsequent surface region evolves as the simulation proceeds at 150 K. Strikingly, during the initial stages of the simulation (100 ps), in the core region the density starts decreasing from its initial value and reaches a value of ∼0.94 g· cm−3, agreeing well with the experimental LDA-ice density at P = 0.1 MPa and T = 77 K68 and with the TIP4P/2005 density (∼0.94 g·cm−3) as reported by Wong et al.35 at ∼150 K. It has also been shown that this TIP4P/2005-based LDA density is independent of cooling rate, and the form adopted here is very much consistent with LDA forms obtained experimentally via hyperquenching.70 Density has been recently identified as a key controlling parameter that affects the nucleation of tetrahedral liquids. The underlying mechanism is termed as a density-wavemodulated process, wherein the nanodroplet surface region (as described above) provides a confining environment by acting as a “hard wall” that results in the formation of “ripple-like” density standing waves due to the reflection of the expanded volume against this hard wall,69 akin to the formation of a stacking fault. Importantly, here, although we have observed enhancement of tetrahedral order of the water molecules in the droplet upon cooling (cf. Figure 3), no sign of crystalline ice (either hexagonal or cubic) is evident (cf. Figure 2); rather, the structure resembles LDA ice. In ice polymorphs (e.g., Ih and Ic), the second-nearest neighbors are known to occupy eclipsed and staggered conformations.71 It would seem that, in the present study, the density-wave driven fluctuation mechanism69 is not sufficient to “push” the structure into adopting this ordering. Instead, the system attends a dynamically arrested state wherein the reorientational ordering is slow enough to label it a “frustrated” process. As mentioned earlier in the Methods section, although the qparameter can indicate the enhanced tetrahedrality of a tetrahedral liquid, it is unable to distinguish between its crystalline, amorphous and liquid phases. In the present study, we have used the CHILL algorithm56 to identify the exact phase attended by the water molecules in the NML-bound nanodroplet. In Figure 6, we have plotted the distribution of the alignment of the local orientational structure, a, for the nanodroplet (cf. Figure 6a) as well as for the bulk water (cf. Figure 6b). We have also presented this distribution for liquid water at 300 K (cf. Figure 6a) which resembles with the liquidlike distribution as reported by Moore et al.56 For the

Figure 6. Plots showing the distribution of the alignment of local orientational structure, a, as defined by eq 5 of ref 56, a(i , j) =

ql(̂ i) . ql(̂ j) | ql(̂ i)| . | ql(̂ j)|

, for l = 3. In (a), we show the distribution for

the water molecules in the nanodroplet, and in (b) the same is shown for a system of bulk water. In (a), we have also shown the distribution for liquid water at T = 300 K. The distribution at a temperature lying in No Man’s Land resembles the distribution for LDA as given in ref 56. However, there is an important difference between the bulk system and the nanodroplet. One can observe that for nanodroplet systems, as T is lowered the systems gradually makes a transition from liquid to LDA. However, at 150 K the LDA-like distribution is quite prominent; for higher temperatures, the distribution is intermediate between liquid-like (300 K) and LDA-like (150 K). On the contrary, for bulk systems the distributions at 150, 175, and 200 K overlap with each other. From this plot, it is quite evident that surface effects in the nanodroplet play a crucial role in comparative inhibition in the formation of LDA ice.

nanodroplet (cf. Figure 6a), as the temperature is lowered from 300 K toward 150 K, the liquid-like distribution corresponding to 300 K gradually transforms to a LDA-like distribution, very much identical to the distribution for LDA ice as reported in ref 56. For the bulk water also, one can notice that the distribution corresponds to the LDA state. The only difference is that, while for the nanodroplet the change from liquid-like to LDA-like distribution with temperature is gradual, in case of bulk, except for 225 K, the distributions at 150, 175, and 200 K overlap with each other. It is quite convincing to notice that in nanodroplets even LDA formation is inhibited and these can be correlated to the surface effects. In any event, these plots confirm that, in the present study, water molecules in a nanodroplet in the No Man’s Land attain a configuration corresponding to low-density amorphous ice. F


Article


However, cases at higher temperature, where water molecules exhibit subdiffusive mobility, also exhibit strong frustration against crystallization due the presence of comparable fraction of five-membered rings, as compared to the corresponding sixmembered fraction (cf. Figure 7).71−75 It is important to mention here that in our studies we did not observe any transition to any of the crystalline polymorphs of ice (Ic or Ih), unlike the recent works as reported by Moore and Molinero,2,15 Russo and Tanaka,71 Sellberg et al.,63 and HajiAkbari et al.76 Although Moore and Molinero2,15 revealed the underlying microscopic mechanism of water crystallization mediated by an intermediate LDA state using a coarse-grained mW potential model of water, Haji-Akbari and Debenedetti76 proposed a direct approach for the rate calculation of homogeneous nucleation of ice using a TIP4P/Ice molecular model. A more recent work by Russo and Tanaka71 used TIP4P/2005 model and proposed a two-state model to describe ice crystallization, where the system of molecules initially attends a high degree of translational order in the second shell in addition to forming five-membered rings of molecules attached together by means of hydrogen bonds. They especially emphasized the “mixed-character” of ice crystallization, where the enhanced translational order favors crystallization but the five-membered rings inhibit them. In the present study, instead of observing crystallization (although not a firm goal per se), we find a transition to an LDA ice in NMLbound nanodroplets acting as a strong inhibition to crystallization. Two limiting factors have played a role here: the size of nanodroplet and the simulation time. As mentioned earlier in the present study and by previous work by Li et al.,41 in nanodroplets surface plays an important role by exerting a positive Laplace pressure onto the core of the droplet and inhibits crystallization; therefore, we find it worthwhile to explore size effects on nucleation, as well as using longer simulation times. Therefore, we ran a simulation for a water nanodroplet with almost 2500 water molecules (in place of 1400 water molecules used for other systems in this paper) at 190 K for 10 μs. The O−O pair distribution function, as well as the distribution of the alignment of the local orientational structure, a, for the nanodroplet, as calculated using the CHILL algorithm, is presented in Figure S4a,b in the Supporting Information. Notably, this also indicates LDA ice. Although LDA has often been associated with an intermediate state for transition to crystalline ice by others15,63,76 in bulk systems mainly, from the present simulations we can conclude that for nanodroplets where surface effects are non-negligible, crystallization is actually inhibited where water molecules remain dynamically arrested in a metastable state like LDA. Such dynamical arrest is enhanced by surface-induced Laplace pressure, and the abundance of five-membered ring structures that require high entropic cost to break “locally favored structures”, as explained by Russo et al.,71 to achieve crystallinity.

The evolution of the nanodroplet orientational ordering (cf. Figure 1), can be described quantitatively by “ring-analysis” (see Methods). Shown in Figure 7 are normalized and averaged

Figure 7. Distribution of various ring configurations at four different temperatures and that for the liquid configuration at 300 K. As temperature is lowered, there is a noticeable change in the number for the pentagonal, hexagonal, and heptagonal configurations. The pentagon distribution saturates at 175 K, the heptagonal one does so at 200 K. Although there is a steady increase of hexagonal rings as temperature is lowered, the system shows frustration for crystallization due to the presence of the pentagonal rings, which require a high entropic cost for ring rupture and reorientation of their neighbors to adopt crystallization-inducing hexagonal rings.

distributions of various n-sided rings (n = 3−10) for various temperatures, compared to the liquid configuration at T = 300 K. The evolution of five-, six-, and seven-membered rings are noticeable. Although the number of six-membered rings increases consistently upon cooling, the number becomes saturated for seven- and five- membered rings at 200 and 175 K, respectively. The ring statistics had been previously been used to characterize various amorphous phases of ice by Martonak et al.39,40 For LDA, they had shown that the network is mainly dominated by six-membered rings while a substantial number of five- and seven-membered rings are also present. Similar to this observation, as mentioned above, we also notice that the rings formed by the water molecules in the nanodroplet consists mainly of the six-membered rings with a substantial presence of five- and seven-membered rings. The presence of five- and six-membered rings are crucial to understanding nucleation from water to ice.71 The existence of locally favored structures71 has also been discussed here in understanding the anomalies observed in water. Six-membered rings can be associated with the “staggered” configurations of water molecules that act as “precursor” for water crystallization. On the other hand, five-membered rings are associated with molecular configurations that are centered around the “eclipsed” configurations and are characterized by fluctuations distinct from the ones generally observed in ice Ih crystals.71 The presence of pentagonal rings inhibits crystallization, owing to the high entropic cost of breaking a hydrogen bond and reorienting the neighboring molecules to obtain crystalline orientation. Therefore, water featuring abundant pentagonal rings can easily adopt a supercooled state without crystallizing per se. Here, at 150 K water molecules exhibit dynamically arrested “glassy” behavior (cf. Figure 4) and therefore lack the mobility necessary for reorientation to achieve crystalline order.

4. CONCLUSIONS In summary, we have explored microscopic views of NML water freezing via classical MD in a nanosized water droplet. Indeed, bearing in mind such a droplet’s fast nucleation, as generally observed in this region of water’s phase diagram, MD is an essential tool. In addition, as remarked previously, the experimental relevance of the aerosol-based approach to NML crystallization, such as outlined by Manka et al.,34 acts as a further motivation for the present study (with the Laplace G



■

pressure retarding nucleation and amorphization vis-à-vis bulksystem counterparts). Understanding water’s wider anomalies near its freezing point is still a matter of debate; many plausible theories/hypotheses exist. Here, we have covered almost the entire NML gamut, spanning from ∼235 to 150 K by sampling at four different temperatures ranging from 150 to 225 K with 25 K intervals. At 150 K, we observe adoption of an arrested glassy state resembling LDA ice owing to its trapped dynamics redolent of glass-like materials. At higher temperatures (200 and 225 K), the system also shows LDA-like O−O RDFs, albeit with tetrahedral-order parameter values approaching liquid-like ones, and with MSDs demonstrating subdiffusive behavior rather than glass-like signatures. Therefore, this could well be regarded as deeply supercooled water that did not show any crystallization during the 1 μs long simulations used here. We associate this frustration of crystallization to surface-induced Laplace pressure onto the core of the droplet, as well as to its five-membered ring structures, which require large entropic costs to be overcome to rupture hydrogen bonds and reorient its neighbors to form six-membered rings to make a transition to a crystalline state. In fact, using a 10 μs long simulation, we confirmed that in such nanodroplets “trapping” to the metastable LDA state is favored over any crystalline polymorph. We also notice a considerable inhibition even for amorphization per se in nanodroplets, when compared to the bulk systems subjected under identical thermodynamic conditions.

■

REFERENCES

(1) Mishima, O.; Stanley, H. E. The relationship between liquid, supercooled and glassy water. Nature 1998, 396, 329−335. (2) Moore, E. B.; Molinero, V. Ice-crystallisation in water’s no-man’s land. J. Chem. Phys. 2010, 132, 244504. (3) Smith, R. S.; Kay, B. D. The existence of supercooled liquid water at 150 K. Nature 1999, 398, 788−791. (4) McClure, S. M.; Barlow, E. T.; Akin, M. A.; Safarik, D. J.; Truskett, T. M.; Mullins, C. B. Transport in amorphous solid water films: Implications for self-diffusivity. J. Phys. Chem. B 2006, 110, 17987−17997. (5) Whalley, E. Scheiner’s Halo: Evidence for Ice Ic in the Atmosphere. Science 1981, 211, 389−390. (6) Tabazadeh, A.; Djikaev, Y. S.; Reiss, H. Surface crystallisation of supercooled water in clouds. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 15873−15878. (7) Thuman, W. C.; Robinson, E. Studies of Alaskan ice-fog particles. J. Meteorol. 1954, 11, 151−156. (8) Jensen, E. J.; Toon, O. B.; Tabazadeh, A.; Sachse, G. W.; Anderson, B. E.; Chan, K. R.; Twohy, C. W.; Gandrud, B.; Aulenbach, S. M.; Hemysfield, A.; et al. Ice nucleation processes in upper tropospheric wave-clouds observed during SUCCESS. Geophys. Res. Lett. 1998, 25, 1363−1366. (9) Rogers, D. C.; DeMott, P. J.; Kreidenweis, S.; Chen, Y. Measurements of ice nucleating aerosols during SUCCESS. Geophys. Res. Lett. 1998, 25, 1383−1386. (10) Heymsfield, A. J.; Miloshevich, L. M. Homogeneous Ice Nucleation and Supercooled Liquid Water in Orographic Wave Clouds. J. Atmos. Sci. 1993, 50, 2335−2353. (11) Finney, J. L.; Hallbrucker, A.; Kohl, I.; Soper, A. K.; Bowron, D. T. Structures of high and low density amorphous ice by neutron diffraction. Phys. Rev. Lett. 2002, 88, 225503. (12) Johari, G. P.; Fleissner, G.; Hallbrucker, A.; Mayer, E. Thermodynamic Continuity between Glassy and Normal Water. J. Phys. Chem. 1994, 98, 4719−4725. (13) Volmer, M. Kinetik der Phasenbildung; Steinkopff, Leipzig, 1939. (14) Matsumoto, M.; Saito, S.; Ohmine, I. Molecular dynamics simulation of the ice nucleation and growth process leading to water freezing. Nature 2002, 416, 409−413. (15) Moore, E. B.; Molinero, V. Structural transformation in supercooled water controls the crystallization rate of ice. Nature 2011, 479, 506−509. (16) Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G. Structure of ice crystallized from supercooled water. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 1041−1045. (17) Moore, E. B.; de la Llave, E.; Welke, K.; Scherlis, D. A.; Molinero, V. Freezing, melting and structure of ice in a hydrophilic nanopore. Phys. Chem. Chem. Phys. 2010, 12, 4124−4134. (18) Moore, E. B.; Molinero, V. Is it cubic? Ice crystallisation from deeply supercooled water. Phys. Chem. Chem. Phys. 2011, 13, 20008− 20016. (19) Hare, D. E.; Sorensen, C. M. Raman spectroscopic study of bulk water supercooled to − 33 °C. J. Chem. Phys. 1990, 93, 25−33. (20) Hare, D. E.; Sorensen, C. M. Raman spectroscopic study of dilute HOD in liquid H2O in the temperature range − 31.5 to 160 °C. J. Chem. Phys. 1990, 93, 6954−6961. (21) Mallamace, F.; Broccio, M.; Corsaro, C.; Faraone, A.; Majolino, D.; Venuti, V.; Liu, L.; Mou, C.-Y.; Chen, S.-H. Evidence of the existence of the low-density liquid phase in supercooled, confined water. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 424−428. (22) Mallamace, F.; Branca, C.; Broccio, M.; Corsaro, C.; GonzalezSegredo, N.; Spooren, J.; Stanley, H. E.; Chen, S.-H. Transport properties of supercooled confined water. Eur. Phys. J.: Spec. Top. 2008, 161, 19−33. (23) Gonzalez-Sanchez, F.; Juranyi, F.; Gimmi, T.; Van Loon, L.; Seydel, T.; Unruh, T. Dynamics of supercooled water in highly compacted clays studied by neutron scattering. J. Phys.: Condens. Matter 2008, 20, 415102.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsearthspacechem.7b00011. Plots of temperature equilibration with time, comparison of results between single and double precisions, comparison of oxygen−oxygen radial distribution function (RDF) between the nanodroplet and the bulk water, and the plots for RDF and the distribution of the alignment of local orientational structure for a water droplet at 190 K and from 10 μs long molecular dynamics simulation run (PDF)

■

Article

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. Telephone: +353 1 716 1646. *E-mail: [email protected]. Telephone: +353 89 950 4210 ORCID

Prithwish K. Nandi: 0000-0003-3458-8853 Christian J. Burnham: 0000-0001-5574-4339 Niall J. English: 0000-0002-8460-3540 Present Address †

(Z.F.) Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS All authors thank Science Foundation Ireland for funding under Grant SFI 15/ERC-I3142. We also thank the anonymous referees for their insightful suggestions to amend the manuscript. H


Article

ACS Earth and Space Chemistry (24) Stan, C. A.; Schneider, G. F.; Shevkoplyas, S. S.; Hashimoto, M.; Ibanescu, M.; Wiley, B. J.; Whitesides, G. M. A microfluidic apparatus for the study of ice nucleation in supercooled water drops. Lab Chip 2009, 9, 2293−2305. (25) Taborek, P. Nucleation in emulsified supercooled water. Phys. Rev. B: Condens. Matter Mater. Phys. 1985, 32, 5902. (26) Butorin, G. T.; Skripov, V. P. Crystallisation of supercooled water. Sov. Phys. Crystallogr. 1972, 17, 322−326. (27) Huang, J. F.; Bartell, L. S. Kinetics of Homogeneous Nucleation in the Freezing of Large Water Clusters. J. Phys. Chem. 1995, 99, 3924−3931. (28) Bartell, L. S.; Chushak, Y. G. In Water in Confining Geometries; Buch, V, Devlin, J. P., Eds.; Springer-Verlag: Berlin Heidelberg, 2003. (29) Wyslouzil, B. E.; Wilemski, G.; Strey, R.; Seifert, S.; Winans, R. E. Small angle X-ray scattering measurements probe water nanodroplet evolution under highly non-equilibrium conditions. Phys. Chem. Chem. Phys. 2007, 9, 5353−5358. (30) DeMott, P. J.; Rogers, D. C. Freezing Nucleation Rates of Dilute Solution Droplets Measured between − 30° and − 40°C in Laboratory Simulations of Natural Clouds. J. Atmos. Sci. 1990, 47, 1056−1064. (31) Stockel, P.; Weidinger, I. M.; Baumgartel, H.; Leisner, T. Rates of Homogeneous Ice Nucleation in Levitated H2O and D2O Droplets. J. Phys. Chem. A 2005, 109, 2540−2546. (32) Kramer, B.; Hubner, O.; Vortisch, H.; Leisner, T.; Schwell, M.; Ruhl, E.; Baumgartel, H. Homogeneous nucleation rates of supercooled water measured in single levitated microdroplets. J. Chem. Phys. 1999, 111, 6521−6527. (33) Duft, D.; Leisner, T. Laboratory evidence for volume-dominated nucleation of ice in supercooled water microdroplets. Atmos. Chem. Phys. 2004, 4, 1997−2000. (34) Manka, A.; Pathak, H.; Tanimura, S.; Wolk, J.; Strey, R.; Wyslouzil, B. E. Freezing water in no-man’s land. Phys. Chem. Chem. Phys. 2012, 14, 4505−4516. (35) Wong, J.; Jahn, D. A.; Giovambattista, N. Pressure-induced transformations in glassy water: A computer simulation study using the TIP4P/2005 model. J. Chem. Phys. 2015, 143, 074501. (36) Winkel, K.; Bowron, D. T.; Loerting, T.; Mayer, E.; Finney, J. L. Relaxation effects in low density amorphous ice: Two distinct structural states observed by neutron diffraction. J. Chem. Phys. 2009, 130, 204502. (37) English, N. J.; Tse, J. S. Reversible pressure-induced crystalamorphous structural transformation in ice Ih. Chem. Phys. Lett. 2014, 609, 54−58. (38) Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E. Phase behaviour of metastable water. Nature 1992, 360, 324−328. (39) Martonak, R.; Donadio, D.; Parrinello, M. Polyamorphism of Ice at Low Temperatures from constant-pressure simulations. Phys. Rev. Lett. 2004, 92, 225702. (40) Martonak, R.; Donadio, D.; Parrinello, M. Evolution of the structure of amorphous ice: from low-density amorphous through high-density amorphous to very high-density amorphous ice. J. Chem. Phys. 2005, 122, 134501. (41) Li, T.; Donadio, D.; Galli, G. Ice nucleation at the nanoscale probes no man’s land of water. Nat. Commun. 2013, 4, 1887. (42) 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. (43) Abascal, J. L. F.; Vega, C. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (44) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (45) 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. (46) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577−8593.

(47) Rahman, A.; Stillinger, F. H. Hydrogen-bond patterns in liquid water. J. Am. Chem. Soc. 1973, 95, 7943−7948. (48) Speedy, R. J. Self-replicating structure in water. J. Phys. Chem. 1984, 88, 3364−3373. (49) Speedy, R. J.; Madura, J. D.; Jorgensen, W. L. Network topology in simulated water. J. Phys. Chem. 1987, 91, 909−913. (50) King, S. V. Ring configurations in a random network model of vitreous silica. Nature 1967, 213, 1112−1113. (51) Pirzadeh, P.; Beaudoin, E. N.; Kusalik, P. G. Structural evolution during water crystallization: Insights from ring analysis. Chem. Phys. Lett. 2011, 517, 117−125. (52) Luzar, A.; Chandler, D. Hydrogen-bond kinetics in liquid water. Nature 1996, 379, 55−57. (53) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. (54) Errington, J. R.; Debenedetti, P. G. Relationship between structural order and the anomalies of liquid water. Nature 2001, 409, 318−321. (55) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Bondorientational orders in liquids and glasses. Phys. Rev. B: Condens. Matter Mater. Phys. 1983, 28, 784. (56) Moore, E. B.; de la Llave, E.; Welke, K.; Scherlis, D. A.; Molinero, V. Freezing, melting and structure of ice in a hydrophilic nanopore. Phys. Chem. Chem. Phys. 2010, 12, 4124−4134. (57) Uttormark, M. J.; Thompson, M. O.; Baez, L. A.; Clancy, P. Solid/liquid cluster recognition in heterogeneous systems. Mol. Simul. 1993, 11, 121. (58) Reinhardt, A.; Doye, J. P. K.; Noya, E. G.; Vega, C. Local order parameters for use in driving homogeneous ice nucleation with allatoms model of water. J. Chem. Phys. 2012, 137, 194504. (59) Wong-Ekkabut, J.; Miettinen, M. S.; Dias, C.; Karttunen, M. Static charges cannot drive a continuous flow of water molecules through a carbon nanotube. Nat. Nanotechnol. 2010, 5, 555−557. (60) Bizid, A.; Oumezzine, A. D. M.; et al. Structure of high-density amorphous water. I. X-ray diffraction study. J. Chem. Phys. 1987, 87, 2225−2230. (61) Bellissent-Funel, M. C.; Bosio, J. T. L. Structure of high-density amorphous water. II. Neutron scattering study. J. Chem. Phys. 1987, 87, 2231−2235. (62) Bowron, D. T.; Finney, J. L.; Hallbrucker, A.; Kohl, I.; Loerting, T.; Soper, M. A. K. The local and intermediate range structures of the five amorphous ices at 80K and ambient pressure: A Faber-Ziman and Bhatia-Thornton analysis. J. Chem. Phys. 2006, 125, 194502. (63) Sellberg, J. A.; et al. Ultrafast X-ray probing of water structure below the homogeneous ice nucleation temperature. Nature 2014, 510, 381−384. (64) Narten, A. H.; Venkatesh, C. G.; Rice, S. A. Diffraction pattern and structure of amorphous solid water at 10 and 77 K. J. Chem. Phys. 1976, 64, 1106−1121. (65) Molinero, V.; Moore, E. B. Water Modeled As an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008. (66) Skinner, L. B.; Huang, C.; Schlesinger, D.; Pettersson, L. G. M.; Nilsson, A.; Benmore, C. J. Benchmark oxygen-oxygen pairdistribution function of ambient water from x-ray diffraction measurements with a wide Q-range. J. Chem. Phys. 2013, 138, 074506. (67) Abascal, J. L. F.; Vega, C. Dipole-Quadrupole Force Ratios Determine the Ability of Potential Models to Describe the Phase Diagram of Water. Phys. Rev. Lett. 2007, 98, 237801. (68) Mishima, O.; Calvert, L. D.; Whalley, E. An apparently firstorder transition between two amorphous phases of ice induced by pressure. Nature 1985, 314, 76. (69) Lu, Y.; Bi, Q.; Yan, X. Density-wave modulated crystallisation in nanoscale silicon films and droplets. J. Chem. Phys. 2016, 144, 234508. (70) Mayer, E. New method for vitrifying water and other liquids by rapid cooling of their aerosols. J. Appl. Phys. 1985, 58, 663. (71) Russo, J.; Tanaka, H. Understanding water’s anomalies with locally favoured structures. Nat. Commun. 2014, 5, 3556. I


Article

ACS Earth and Space Chemistry (72) Tanaka, H. Bond orientational order in liquids: Towards a unified description of water-like anomalies, liquid-liquid transition, glass transition, and crystallisation. Eur. Phys. J. E: Soft Matter Biol. Phys. 2012, 35, 1−84. (73) Tanaka, H. Importance of many-body orientational correlations in the physical description of liquids. Faraday Discuss. 2014, 167, 9− 76. (74) Tanaka, H. A simple physical model of liquid-glass transition: Intrinsic fluctuating interactions and random fields hidden in glassforming liquids. J. Phys.: Condens. Matter 1998, 10, L207−L214. (75) Shintani, H.; Tanaka, H. Frustration on the way to crystallisation in glass. Nat. Phys. 2006, 2, 200−206. (76) Haji-Akbari, A.; Debenedetti, P. G. Direct calculation of ice homogeneous nucleation rate for a molecular model of water. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 10582−10588.

J


Ice-Amorphization of Supercooled Water ... - ACS Publications

Recommend Documents