Acetonitrile cluster solvation in a cryogenic ethane-methane-propane liquid: Implications for Titan lake chemistry

Discussions of extraterrestrial life often focus on whether liquid water is available because life on Earth is intimately intertwined with the chemical properties of water. However, proponents of the notion that water may not be unique in this regard suggest that Saturn’s moon Titan is a prime candidate for testing this hypothesis.¹ Titan’s surface temperature of between 90 and 95 K (varying with latitude and season) is too cold for liquid water, but data from the Cassini orbiter and Huygens surface probe indicate that alkane liquids on the surface form lakes and seas.² Fluvial features have also been observed, suggesting the existence of a hydrological system based on methane,³ the atmosphere’s second largest component after molecular nitrogen.⁴ Chemistry in Titan’s atmosphere driven by solar-ultraviolet radiation and energetic particles in Saturn’s magnetosphere leads to the formation of a variety of organics, including nitrogen-bearing molecules.⁵ Their solvation, possibly as aggregates, in Titan’s seas could be the first step towards the development of an alternative form of life in which the conventional roles of hydrophobic and hydrophilic molecules are reversed.⁶ A case of particular interest is acetonitrile (CH₃CN, abbr. ACN), which has been detected in gaseous form in the atmosphere and possibly as ice on the surface.⁷ Gaseous ACN molecules may aggregate into ordered clusters via dipole interactions, precipitate, and drop to the surface or, alternatively, may aggregate during fluvial transport of condensed ACN deposits on the surface. Our intent is to explore in detail the solvation characteristics of a species known to be produced in Titan’s atmosphere in the gas phase and inferred to be part of the aerosol that rains down onto the surface and into the lakes. The question of the extent to which differing vapor pressures might lead to aerosols of relatively pure photochemical products remains an open one in the Titan and Cassini literature. For example, Anderson et al.⁸ review the arguments for and against the presence of an essentially separate aerosol layer of cyanoacetylene ice in Titan’s atmosphere. Very recently Ennis et al.,⁹ made laboratory IR spectra of pure ACN ice for comparison with Cassini data. While the match was poor, the authors note that one cannot distinguish between phases comprised of a mix of ACN with other species versus tiny grains of segregated pure ice. In the latter case, upon sedimentation into the lakes, these agglomerations might well separate into the individual pure ice through chemical and physical (e.g., turbulent stirring) effects. Such aggregates, while of pure or of mixed constituents, might provide unique microscopic polar environments capable of facilitating chemistry normally impossible in the otherwise nonpolar alkane lakes.

Under Earth ambient conditions, it is well known that ACN is highly soluble in water, dissolving into single molecules,¹⁰ at least at dilute concentrations. In contrast, at room temperature it is insoluble in hexane and similar organic liquids.¹¹ However, the solvation behavior of ACN in (liquid) solutions of small alkanes, specifically at low temperatures at near Earth ambient pressures, has not previously been studied. Here, we explore the solvation behavior of ACN in small alkanes at cryogenic conditions using classical molecular dynamics (MD) simulations. Such simulations are well suited for this purpose because they allow the exploration of solvation chemistry at Titan’s cryogenic temperatures (∼90–95 K) and near its ambient pressure (∼146.7 kPa). Specifically, we used simulations to compute the solvation characteristics of ACN clusters in pure methane and ethane solutions as well as in an ethane-methane-propane (EMP) mixture that is a possible composition of some of the smaller lakes of Titan.¹² 

We used a potential of mean force (PMF) to quantify the energetics of solute molecules along the pathway from vapor to liquid. This was computed by integrating the mean force of the solvent acting on a solute particle along a reaction coordinate between the two end states. The integration was accomplished using the Umbrella Integration method,¹³ a derivative of the “weighted histogram analysis method” (WHAM).¹⁴ This method allowed us to investigate energetics through the interface, rather than only the solvation difference between the bulk vapor and bulk liquid states.¹⁵ The results suggest that small alkane molecules in liquids at cryogenic temperatures surround small ACN clusters such that ACN aggregates do not separate into individual molecules in this environment.

We used the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) code for the classical MD simulations¹⁶ and integrated the equations of motion with the velocity Verlet algorithm. Long-range Coulombic interactions were computed using the smooth particle mesh (PPPM) method.¹⁷ Intermolecular and intramolecular interactions were modeled using optimized potentials for liquid simulation—all atom (OPLS-AA) force fields¹⁸ with a combined Coulomb and Lennard–Jones interaction that used a cutoff of 1.5 nm for the nonbonding energy between two atoms. The corresponding force equations and parameters for intermolecular OPLS-AA charges and Leonard-Jones interactions, as well as those for the intramolecular bonds, angles, and dihedrals, were taken from the literature.^18–20 

Periodic boundary conditions were used. Supercell boundaries were placed at −1.5 nm and 1.5 nm for both the x- and y-directions. The depth in the z-direction depended on the molar volume of the specific liquid but was set to be about four times the cell width to ensure that a bulk state was achieved in runs with liquid/vapor surfaces. To include surfaces, the supercells were extended in the z-direction to create an effective vapor phase of at least 2.0 nm between each of the two (periodic) surfaces and the liquid phase. Thus, the center-of-mass of the alkane liquid was at the origin (0,0,0) and bounded by regions of vapor, such that the periodic boundary in the z-direction only involved the vapor phase.

The most ethane rich of Titan’s liquid bodies is Ontario Lacus, with about 40% ethane; the larger liquid bodies are much more methane-rich.¹² However, smaller lakes may evaporate seasonally, or on longer time scales, and as they do they will become ethane- and propane-rich by vapor pressure distillation. We are interested in the solvation of ACN clusters in such liquids because it is likely to be much higher than in the methane-rich lakes and seas characterized by Cassini. Also, it is in the drying and wetting cycling of such lakes, by analogy with prebiotic systems on Earth, where the most interesting chemistry may occur. We therefore chose an EMP liquid composition corresponding to mole percentages of 82% ethane, 10% methane, and 8% propane. Although nitrogen will be important in more methane-rich systems, it will be negligible for these small lakes having a solubility of <5%. For this nitrogen solubility and the size of the solvation shell, there would only be on average one N₂ molecule present half the time. Cassini radar data show that the more methane-rich compositions of the larger seas may reach N₂ molecular abundances of 10%-15%, but this remains a small enough fraction that, on average, only one N₂ molecule would be present in the solvation shell and the basic interactions would not be affected. Overall, the liquid system contained 1200 alkane molecules which comprised of 984 ethane molecules, 120 methane molecules, and 96 propane molecules. The pure methane and ethane systems each contained 1200 molecules. Initial starting structures of the liquid alkane systems were built by packing each of the molecule types into a fixed box size using Packmol.^21,22 The alkane liquid-vapor systems were brought to equilibrium at 94 K using the NVT (canonical) ensemble with a Langevin thermostat during a 2 ns simulation with a 1 fs time step. A 1 ns NVE (microcanonical) ensemble simulation with a 1 fs time step and a subsequent 2 ns simulation with a 2 fs time step were used to verify that equilibrium had been reached.

Density profiles of the total liquid and of each individual molecular species were determined for the liquid-vapor system by binning the center of mass positions of the molecules in the z-direction (Figure 1). Snapshots of this distribution were collected every 25 steps during a 2 ns NVE production run with a 2 fs time step.

Capillary waves between liquid–vapor interfaces are known to dominate the interfacial thickness, w, determined in MD simulations by fitting interfacial results to a theoretically derived fundamental equation²³ that uses the Gibbs dividing surface, z_G,

The Gibbs dividing surface, z_G, is the z-position at half the total density along the density profile as it drops off from the bulk to the vapor phase. Equation (1) is fit to the density profile along z to find the interfacial thickness w.

Surface tension, γ, was determined as the difference between the pressure tensor components perpendicular and parallel to the interface,^24,25 with a dependence on the length of the simulation box in the direction normal to the interface, z_L,

Umbrella Integration uses restraint dynamics along a trajectory divided into windows. A bias force is applied to the molecule or cluster with the bias potential $Ub=12K(ξ−ξeq)2$⁠, where K is the spring constant and $ξeq$ is the bias’s equilibrium position. The probability distribution, P(ξ), of finding the particle at a position ξ within a specified window leads to the Helmholtz energy (A) of the biased system

The Helmholtz energy without the bias potential requires removing the bias (or restraint force)

where F, the integration constant, is removed by taking the partial derivative of Equation (4) relative to $ξ$⁠,

An expansion of the biased probability distribution about an equilibrium leads to a leading term in the form of a Gaussian distribution

Thus, the average and variance of ξ are obtained from the window of the biased simulation, which is centered at $ξeq$ such that

The effective mean force is then determined by umbrella sampling techniques that use independent bins to combine the simulation windows.¹³ Integration of Equation (7) along the trajectory gives the Helmholtz energy as a function of ξ.

Trajectories to obtain the PMFs used sampling prisms centered 0.1 nm apart, where the strength of the spring constant defined the width of the prism. Simulations employed an NVT ensemble at 94 K (the average high temperature of Titan²⁶) using a Langevin thermostat. The initial equilibration sampling at a sampling position was performed using 400 000 1 fs time steps. The ACN molecule or cluster was “placed” into each trajectory prism center such that each simulation in each prism position was independent of all others. Each of the independent initial trajectory runs was extended for additional 800 000 2 fs time steps from which the final PMF data were sampled every 50 steps. These extended runs were necessary because the EMP liquid at Titan temperatures exhibits slow relaxation (glassy) dynamics. Convergence was confirmed by extending a few trajectories for another 800 000 2 fs time steps.

To compute entropy and internal energy solvation thermodynamics, additional simulations were conducted at 94, 99, and 104 K for the single ACN molecule and for clusters of 4, 9, 12, and 16 molecules. Assuming $ΔU$ and $ΔS$ is independent of T in this small range of temperature, we plotted $ΔA$ vs. T and $ΔA$ vs. 1/T, linearly fit the results, and calculated the slopes, which correspond to the negative entropy (⁠$−ΔS$⁠) and internal energy (⁠$ΔU$⁠), respectively.

We also determined the PMF for a pair of ACN molecules as they were separated from each other in the bulk liquid at 94 K. Both the interactions with the solvent and the pairwise interaction of the two solute molecules were quantified. These simulation trajectories were carried out as described above.

Figure 1 shows the total and component density profiles computed for the EMP liquid comprised of 82, 10, and 8 mol. % of ethane, methane, and propane, the estimated composition of Titan’s lakes. The overall liquid density is 0.66 g/cm³, where ethane, methane, and propane contribute densities of 0.54, 0.03, and 0.09 g/cm³, respectively (see Table I). The component density profiles reveal that methane is enriched and, correspondingly, propane and ethane are depleted, at the surface. This is explained by the relative surface tensions computed for the three components (Table I). Methane migrates preferentially to the surface because it has the lowest surface tension (Table I), whereas propane, with the highest surface tension, is the most depleted. In addition, surface thicknesses were determined for the simulated EMP liquid and were compared to those of the simulated pure liquids of ethane, methane, and propane. Overall, the EMP liquid has a surface thickness and density similar to that of the ethane liquid (Table I).

The PMFs of clusters comprised of 1, 2, 3, 4, 5, 6, 9, 12, and 16 ACN molecules in vapor and liquid are shown in Figure 2 as a function of the position of cluster center-of-mass. The plots start in the vapor phase at 6.5 nm above the liquid surface, contact the liquid surface at roughly 5 nm, proceed into the bulk, and end at 2.2 nm from the center of the system. The average radii of the clusters vary from 0.26 nm to 0.52 nm for the 4 and 16 sized clusters, respectively. The magnitude of the decrease in Helmholtz free energy increases with cluster size, which indicates that the larger clusters will likely not dissociate to smaller clusters or into single molecular entities. The curves shift rightward with increasing cluster size due to earlier contact between the liquid and the clusters of increasing radius. The clusters maintain a constant Helmholtz free energy once they are completely submerged in or surrounded by the solvent in a solvation cavity.

While the thermodynamic results indicate that the non-polar liquid would engulf the clusters, regardless of whether they remain on the surface or sink due to other forces, the clusters provide a polar environment to which other polar particles could attach themselves and, perhaps, even be encased within the clusters themselves.

While we focused on an ethane-rich EMP mixture, recent comparisons of microwave absorptivity between the high-latitude northern sea “Ligeia Mare” and the southern hemisphere lake “Ontario Lacus” suggest that other lake/sea compositions, including some of primarily methane,^27,28 are also found on Titan. To see if such variations would significantly affect our results, we conducted additional simulations to compute the PMFs of a single ACN molecule in pure methane, pure ethane, and the EMP liquid (Figure 3). We found that the energetic benefits in the dilute limit of solvation are similar for all three liquids, as well as to the results in which the addition of ethane to methane leads to small changes.¹⁵ Note, however, that the saturated mole fraction of ACN is significantly higher in ethane than in methane.²⁹ (The lateral shifts between the curves are due to the different molar volumes of each simulated liquid and are irrelevant in interpreting the thermodynamics.) Specifically, the energetic differences were small compared to the numerical variations of the PMF simulations between runs (±0.50 kcal/mole). We interpret this as evidence that ACN aggregates are generally surrounded by a wide range of ethane-methane mixtures, and it indicates that our results are robust against compositional uncertainties.

The results presented above demonstrate that larger ACN clusters in the solution are preferred over smaller clusters in the alkane liquid. An additional study was used to determine to what extent larger clusters are stable in the solution. Since it was impractical to directly determine the free energy change for each of the large number of possible dissociation products, we considered the solvent-ACN and ACN-ACN interactions independently to assess the factors that control the cluster dissociation energies in the EMP solvent. This allows the exploration to what extent a solute will dissociate into its molecular components in a solvent. [In principal, the solubility of a solute is considered to be the ability of a solute to dissociate from its aggregates into its molecular constituents.]

The change in the solvent-ACN free energy was determined by calculating the per-monomer solvent-ACN free energies obtained by dividing the total Helmholtz free energies of solvation (Figure 2) by the number of monomers in the cluster (Figure 4). The Helmholtz free energies were further separated into entropic and internal energy contributions by computing the temperature-dependencies. The Helmholtz energy for the monomer, −2.7 kcal/mol, represents the favorable interaction of a fully solvated ACN, and the per-monomer value decreases in magnitude as the cluster size increases as expected due to the decreasing surface-to-volume ratio and resultant reduction in the average amount of solvent-monomer contact. For example, the per-monomer solvation Helmholtz free energy for the 16-mer was −1.5 kcal/mol, representing a net destabilizing effect of 1.2 kcal/mol per monomer due to reduced solvation interactions in the larger cluster.

This solvation effect is countered by the attraction between ACN monomers. To assess this, we computed the PMF of an ACN molecule as it was separated from its partner in an ACN dimer, either in the solvent (Figure 5; middle curve, solid) or vacuum (Figure 5; lower curve, dotted). The difference between them (Figure 5; top curve, dot-dash) represents the solvent rearrangement energy when the dimer is formed. The solvent rearrangement does not favor the formation of a dimer as it steadily increases in energy. The key point is that the cost in Helmholtz energy needed to rearrange the solvent about the pair of ACN molecules as they are brought together is overcome by their interaction energy giving a net (attractive) free-energy change of −2.8 kcal/mol. Examination of the molecular configurations in the simulations showed that this was primarily due to dipole-dipole interactions: The molecules are in a side-by-side anti-parallel alignment for the PMF local minimum at ∼0.35 nm and are in an end-to-end parallel alignment for the local minimum at ∼0.55 nm. Since an ACN n-mer has at least n − 1 pair-wise ACN-ACN interactions, this free energy of interaction will greatly outweigh the difference in per monomer free energies of solvation. Therefore, we anticipate that the stability in the solution will continue to increase with cluster size, and also that the dipole-dipole interactions are strong enough to facilitate the formation of dimers and larger clusters from dissolved monomers. In addition, we expect including polarization effects will lead to slightly less solvation stability, as this would emphasize the difference between polar and non-polar interactions.¹⁵ However, this would not affect the solvation mechanism.

In conclusion, MD simulations and PMF calculations proved useful for investigating the solvation in an EMP model of Titan’s alkane lakes at Titan’s cryogenic temperatures. These methods may be extended to study other molecules of interest for Titan chemistry. We found that the EMP liquid surrounds acetonitrile clusters without dissociating into individual molecules. Furthermore, the relative free energy contributions of ACN solvation and ACN-ACN interactions suggest that the stability of ACN clusters in the bulk alkane liquid increases with cluster size so that individual ACN monomers will form clusters when dissolved in the liquid. Even though the density of ACN (0.786 g/cm³ at 298 K) is greater than that of the liquids considered here, aggregates smaller than the mass needed to reach terminal velocity in the liquid would likely be suspended. Simulations in pure methane, ethane, and propane solutions were essentially the same as those in the EMP mixture, suggesting that these conclusions are robust and largely applicable to other solvent mixtures. This suggests that ACN-facilitated polar chemistry could occur in a variety of possible Titan lake compositions. We suggest that future studies investigate the possibility and time scale of precipitation of solid ACN crystals formed from large clusters. They could also extend the basic model to include a nitrogen vapor pressure approximating that in Titan’s atmosphere.

The authors would like to acknowledge funding from the John Templeton Foundation and the use of the high performance computing facility at the University of Arizona. L.R.C. would like to acknowledge the hospitality of both the Department of Chemistry and Biochemistry at the University of Arizona and of the Department of Astronomy at Cornell University.