Communication: Existence and control of liquid polymorphism in methanol under shear

It is well known that a mixture of two compounds can separate into two liquid phases. However, such a phenomenon had long been thought to be impossible for a pure substance. About 20 years ago,¹ it was suggested that the existence of two liquid forms of water could account for the anomalous properties of supercooled water.^2–4 The idea was that, at low temperature, below the line of homogeneous ice nucleation, water could exhibit two liquid polymorphs, of different densities and structures, with a critical point terminating the coexistence line between the two liquids. This hypothesis provided an explanation for unexpected experimental results for the isothermal compressibility,² heat capacity,³ and thermal expansion coefficient⁴ of supercooled water. Since then, the concept of liquid polymorphism has become increasingly popular^5–23 and liquid polymorphs have been reported for a wide range of systems, including elemental liquids, like phosphorus¹⁹ and silicon,^20–22 as well as molecular liquids, such as water^8–10,15,18 and alcohols.²³ However, since the liquid-liquid coexistence line is located deep in the supercooled region, the formation of the low temperature liquid polymorph is pre-empted by rapid crystallization.⁶ Therefore, studying and understanding the impact of liquid polymorphism has remained extremely challenging.

Since thermodynamics favors crystallization, it is necessary to study the system under nonequilibrium conditions to be able to observe the formation of the low temperature liquid polymorph. A couple of strategies have been designed to achieve this in experiments. Choosing a set of conditions such that the system is controlled by kinetics rather than by thermodynamics provides a first route to explore liquid polymorphism. For instance, by confining water to droplets with diameters of 1-10 mm under high pressures (0.1–1 GPa), Mishima and Stanley⁷ succeeded in kinetically suppressing crystal nucleation and provided an estimate for the liquid-liquid critical point. A second possible route consists in applying a short perturbation that pushes the system out of equilibrium and in analyzing the relaxation pathway. This method led to the transient formation of the low temperature liquid polymorph of silicon, which was achieved by subjecting a silicon sample to laser pulses of femtosecond duration.²¹ While these studies have greatly advanced our understanding of this phenomenon, they also impose a limit on the set of conditions and time during which liquid polymorphs can be examined. In this work, we propose to generate nonequilibrium steady states of the liquid polymorphs of methanol under shear. This approach has several advantages, as it can be applied to arbitrarily large systems (e.g., here on the bulk) and the liquid polymorphs can be stabilized and studied for arbitrarily long times (as long as shear is applied to the system), and also provides insight into how liquid polymorphism can be triggered and controlled by tuning the value of the applied shear.

We use nonequilibrium molecular dynamics (NEMD) simulations to study the response of methanol when subjected to a constant shear stress.²⁴ The molecular model we use for methanol^25,26 is a coarse-grained model based on a continuous soft-core attractive potential.^27,28 It has been shown to exhibit two metastable liquid polymorphs, of the same free energy, below the melting point. While one of the liquid forms (form L₁) is stable at temperatures above melting, the second liquid polymorph (form L₂) never becomes stable and is only observed below melting and for a very short time as it undergoes rapid crystallization.²⁵ We carry out simulations on systems of 4096 methanol molecules for densities ranging from 0.7 g/cm³ to 0.89 g/cm³ at a temperature of 168 K (close to the liquid-liquid critical point of the model). Using the system of reduced units proposed of Huš and Urbic,²⁵ these conditions correspond to reduced densities ρ^* between 0.22 and 0.28 and to a reduced temperature of 0.52 in the system of reduced units. We use this system of reduced units throughout the rest of this work. Simulations are carried out using an in-house code. We model methanol as a linear rigid body, integrating with a five-value Gear predictor-corrector algorithm the equations of motion for the center-of-mass and for the rotation of the molecule.^29,30 From a practical standpoint, we use a time step of 5 × 10⁻⁴ in reduced units. The NEMD method we use is based on the SLLOD equations of motion, developed by Evans and Morriss,²⁴ and a feedback mechanism is employed to maintain the shear stress constant as discussed in prior work.^31,32 Heat dissipation is accounted for through the use of a Gaussian thermostat that maintains the streaming kinetic energy of the centers-of-mass constant, which has been shown to provide results in very good agreement for the range of shear rates studied here (i.e., below 1 in reduced units).²⁹ For each NEMD simulation, we start by carrying out a first run of 1 × 10⁶ time steps to allow the system to reach a steady state. We check that the properties (elements of the pressure tensor, interaction energy, and viscosity) have converged and that we obtain the linear velocity profile expected for the liquid undergoing a shear flow. We then perform a production run of 2 × 10⁶ time steps over which average properties are calculated.

We start by commenting on the response of the system, as a function of density, when subjected to a constant shear stress. As shown in experiments³³ and prior simulation work,^34–37 the Steinhardt order parameter³⁸ Q₆ provides a measure of the amount of crystalline order within a system and is calculated as

where $r^ij$ is the unit vector joining two neighboring atoms i and j, which are less than a distance of 2.1 from each other (corresponding to the first minimum of the pair distribution function for the liquid), $Y6m(r^ij)$ is the spherical harmonics, and N_b(i) is the number of neighbors for molecule i. Q₆ takes values greater than 0.3 for crystal phases and vanishes in the liquid.³⁴ 

At low densities (ρ^* = 0.23) and at equilibrium (shear rate γ^* = 0), methanol is under the L₁ liquid form, with a value for the order parameter Q₆ around 3 × 10⁻². As shown in Fig. 1, the system flows as soon as a shear stress is applied and a linear flow profile develops across the liquid. We add that the value for Q₆ remains very low over the whole range of shear stress. On the other hand, at high densities (ρ^* = 0.27), methanol is solid at equilibrium. This can be seen on the right panel of Fig. 1 with the much larger values found for Q₆, in excess of 0.3. In this case, applying a shear stress above 0.2 leads to a complete shear-induced melting of the solid, allowing us to obtain a flowing liquid. This is confirmed by the drop in Q₆ value by an order of magnitude. We discuss in the next paragraph the results obtained on sheared liquids over the whole density range, before providing a full picture of the nonequilibrium phase behavior of methanol under shear.

We begin by analyzing the steady-state response of the fluid to the applied shear. Figure 2 shows the results obtained for the steady-state shear viscosity as a function of density for a given value of the shear rate. Shear viscosity values are of the same order (between 0.9 and 1.5) across the density range, showing that we have succeeded in obtaining a liquid phase for all densities between 0.22 and 0.28. The plot shows the overall behavior generally expected for a liquid, i.e., an increase of the shear viscosity with density. However, a closer inspection of the results establishes that the rate at which the shear viscosity increases with density exhibits two regimes, with a crossover for a reduced density of 0.25. The shear viscosity increases linearly for a low-density regime (0.22-0.25), which corresponded to the L₁ liquid form at equilibrium. For the high-density regime (densities greater than 0.25), we still have a linear increase in the shear viscosity, but the slope is greater than in the low-density regime. This means that, for densities above 0.25, which corresponded to the domain where rapid crystallization pre-empted the formation of the L₂ liquid polymorph at equilibrium,²⁵ liquid methanol exhibits a different response to the applied shear. An analysis of the results over the whole range of shear rates further shows that the density at which the crossover between the two regimes takes place remains the same (0.25), regardless of the value of the shear rate. This shows that the existence of the crossover is not shear rate-dependent, but that it is intrinsic to the fluid and stems from the change in density. The results suggest the existence of two liquid forms of methanol under shear and that the crossover observed for the shear viscosity can be interpreted as a signature of the transition from one liquid form to the other. They also establish these two liquid forms can be stabilized by shear, since the results presented here are obtained in nonequilibrium steady states. We finally add that the temperature considered here is very close to the liquid-liquid critical point,²⁵ which is consistent with the absence of a noticeable jump in the viscosity at the liquid-liquid transition.³⁹ 

The next step consists in examining the effect of shear on the structure of liquid methanol. The left panel of Fig. 3 shows the pair distribution function between the O(H) sites of methanol at equilibrium and under increasing shear for the liquid polymorph L₁ for a density of 0.23. For the range of shear rates considered here, applying shear does not lead to dramatic changes in the steady-state structure of the liquid, and the structure of the L₁ liquid polymorph of methanol under shear remains essentially the same as at equilibrium. Figure 3 shows that increasing the applied shear results in a slightly less structured fluid, as evidenced by the changes in the height of the second peak and, to a lesser extent, of the first peak. However, the distances for which this function exhibits extrema are not affected by the applied shear, showing that the generic features of the L₁ polymorph are retained in the nonequilibrium steady states generated in this work.

We compare in the right panel of Fig. 3 the structural features of liquid methanol under shear for a density of 0.23, below the crossover density, and for a density of 0.27, beyond the crossover density. For both densities, the structural features of the fluid are typical of the liquid state, providing further confirmation that the high density system of methanol has undergone complete shear-induced melting. The structural features observed for the liquids at the two densities differ markedly from each other. The first difference lies in the relative height of the first two peaks. At low density, the pair distribution function exhibits a lower maximum for the first peak than that for the second peak, indicating that there are fewer neighbors in the first coordination shell than in the second coordination shell. However, at high density, we observe the reverse, with a higher maximum for the first peak than for the second peak. The second difference is subtler. Looking more closely at the pair distribution function for a density of 0.27, we notice that this function seems to display an inflection around a distance of about 3. To ascertain this point, we calculate the derivative of the pair distribution function, with respect to the distance, and compare the results obtained in the two regimes (bottom panel on the right of Fig. 3). This confirms that the two functions exhibit a different behavior around a distance of 3, with the derivative displaying a dip for a density of 0.27. These results show that methanol under shear can exist in two liquid forms of different structure and that the crossover in transport properties of Fig. 2 serves as a marker of the nonequilibrium liquid L₁-liquid L₂ transition. This establishes the existence of liquid polymorphism under nonequilibrium conditions, when the system is driven away from equilibrium by an external field. Furthermore, the external field provides a means to stabilize, in the steady state, not only the L₁ liquid polymorph but also the L₂ polymorph, which was not stable at equilibrium.

The formation of the L₂ polymorph of methanol is pre-empted by crystallization at equilibrium. However, subjecting methanol to shear allows us to determine sets of conditions for which steady states of the L₁ and L₂ polymorphs can be obtained. This means that the nonequilibrium phase diagram for methanol under shear differs significantly from the equilibrium phase diagram. We present in Fig. 4 the phase diagram for methanol under shear, as predicted by the simulation results. This diagram, plotted in the pressure-shear stress plane, shows the domains where steady states of each of the 3 phases (liquid polymorph L₁, liquid polymorph L₂, and solid) can be obtained. We have seen in the first part of this work that the response of liquid methanol to shear exhibits a crossover from the low density (L₁ polymorph) regime to the high density (L₂ polymorph) regime. Therefore, the location of the liquid L₁-liquid L₂ line can be readily determined by varying the applied shear and by reporting the simulation data (pressure and shear stress) for the crossover points. For the liquid L₂-solid line, we consider high density systems of methanol and determine the minimum shear stress for which a linear flow profile, and thus a complete shear-induced melting, is observed. The nonequilibrium phase diagram in Fig. 4 summarizes the domains where each of the 3 phases can be stabilized. This phase diagram also establishes how the control of liquid polymorphism can be achieved in methanol under shear. For instance, when the L₁ liquid form is subjected to a given shear stress (of the order of 0.2 or greater), increasing the pressure will lead to the formation of a steady state of the L₂ polymorph. The phase diagram also suggests two possible paths to trigger the L₂ → L₁ transition. This can be done through a decrease in pressure at fixed shear stress or through an increase in the applied shear at fixed pressure. The first path for the L₂ → L₁ transition is easily understood in terms of the decrease in density associated with the decrease in pressure at fixed shear. The second path, which is obtained by increasing the shear at constant pressure, results from the compensation of two effects. Increasing the shear stress leads to an increase in the potential energy and thus to the contribution of the potential energy to the pressure.²⁴ To counteract this effect and to maintain the pressure constant, the liquid undergoes dilation, which results in a decrease in density and hence to the L₂ → L₁ transition. Future work will focus on elucidating the nature (discontinuous vs. continuous) of the liquid-liquid transition under shear.

The results presented in this work point to the existence and stabilization of liquid polymorphs in sheared liquids. Furthermore, our findings, obtained here in the case of a coarse-grained model for methanol, show how nonequilibrium liquid-liquid transitions can be triggered and controlled, as indicated in the nonequilibrium phase diagram of methanol under shear. The results presented here also suggest how future research, based on the nonequilibrium approach outlined here, can provide new ways to achieve a better understanding of the properties of deeply supercooled liquids for a wide range of elemental and molecular liquids, including water.

Partial funding for this research was provided through a NSF CAREER Award No. DMR-1052808.