Molecular Structure, Dynamics, and Vibrational Spectroscopy of the Acetylene:Ammonia (1:1) Plastic Co-Crystal at Titan Conditions

The Saturnian moon Titan has a thick, organic-rich atmosphere, and condensed phases of small organic molecules are anticipated to be stable on its surface. Of particular importance are crystalline phases of organics, known as cryominerals, which can play important roles in surface chemistry and geological processes on Titan. Many of these cryominerals could exhibit rich phase behavior, especially multicomponent cryominerals whose component molecules have multiple solid phases. One such cryomineral is the acetylene:ammonia (1:1) co-crystal, and here we use density functional theory-based ab initio molecular dynamics simulations to quantify its structure and dynamics at Titan conditions. We show that the acetylene:ammonia (1:1) co-crystal is a plastic co-crystal (or rotator phase) at Titan conditions because the ammonia molecules are orientationally disordered. Moreover, the ammonia molecules within this co-crystal rotate on picosecond timescales, and this rotation is accompanied by the breakage and reformation of hydrogen bonds between the ammonia hydrogens and the {\pi}-system of acetylene. The robustness of our predictions is supported by comparing the predictions of two density functional approximations at different levels of theory, as well as through the prediction of infrared and Raman spectra that agree well with experimental measurements. We anticipate that these results will aid in understanding geochemistry on the surface of Titan.


INTRODUCTION
Titan, the largest moon of Saturn, is of great interest to planetary science and prebiotic chemistry communities because of its similarities to Earth 1-5 . Like Earth, Titan has a thick atmosphere, but it is composed mainly of nitrogen and methane, and atmospheric photochemistry creates a rich inventory of organic molecules on Titan [5][6][7][8][9] . Titan's low surface temperature of approximately 94 K -a major difference between Titan and Earthleads to condensation of organic molecules from the atmosphere and onto the surface. The surface of Titan has stable liquids composed mainly of hydrocarbons, which undergo seasonal rainfall cycles analogous to water-based cycles on Earth 10 . In addition to organic liquids, organic solids can form through condensation or surficial processes at Titan surface temperatures. Indeed, data from the Cassini-Huygens mission has provided strong evidence to suggest that Titan's surface is dominated by organic solids [11][12][13][14] . These organic crystals, also referred to as cryominerals, could play a major role in geology, geochemistry, and even prebiotic chemistry on Titan's surface, similar to the importance of Earth's minerals in the terrestrial analogs of these processes [15][16][17] .
Titan's cryominerals are an exciting class of crystalline materials that exist as pure phases or multicomponent solids composed of more than one type of organic molecule -co-crystals 15,16,18 . These molecular co-crystals are typically soft because they are held together mainly by relatively weak intermolecular interactions, such as hydrogen bonding, π-π, or van der Waals interactions. These different types of intermolecular ina) rick.remsing@rutgers.edu teractions act at different energy scales, such that molecular crystals often display rich phase behavior consistent with thermal excitations disrupting various interactions at different temperatures. For example, molecular crystals are translationally and orientationally ordered at low temperatures, resulting in a crystal phase. However, as the temperature is increased, orientational degrees of freedom can become disordered without the solid melting. As a result, the molecular solid is translationally ordered but orientationally disordered, resulting in a 'rotator' or 'plastic crystal' phase 19 .
This complex phase behavior of cryominerals must be quantified to gain an understanding of their role in surface processes on Titan. For example, a plastic crystal typically displays different mechanical, thermodynamic, and chemical properties than those of the perfect crystal. Crystal plasticity usually softens elastic constants and increases mechanical flexibility 20 . Anomalies in the density, thermal expansivity, specific heat, and electrical properties of molecular solids have also been directly connected to the presence of plastic phases 16,21 . The dynamic disorder within a plastic phase can facilitate transport through a crystal via molecular paddlewheel mechanisms 19,[22][23][24][25][26] . Seismic discontinuities across Earth's mantle have also been related to the plastic crystalline phase within deep Earth minerals 27 . Together, these effects may ultimately influence metamorphism, fracturing, crack propagation, and erosion rates on Titan. Connecting these molecular level insights to the large-scale surface features constitutes an essential step in constructing a foundational understanding of the geological evolution and possible prebiotic chemistry of Titan. 16,20 In this work, we investigate a model Titan cryomineral, the acetylene:ammonia (1:1) co-crystal. The acetylene:ammonia (1:1) co-crystal is a prime candidate to display a plastic phase because both components exhibit plastic crystal phases [28][29][30][31][32] . Moreover, the component molecules that make up the acetylene:ammonia (1:1) cocrystal have been confirmed in Titan's atmosphere and on the surface by instruments onboard Cassini [33][34][35][36][37][38] . Acetylene may be the second most abundant photochemical product in Titan's hazy atmosphere 6,39 and has been linked to evaporite deposits 33,40 and equatorial dunes 12,39,41 on Titan's surface. Along with its pure form, acetylene is also known to co-crystallize with many small organic molecules forming a variety of potential Titan cryominerals 15,16,42 . Like acetylene, ammonia has been detected in Titan's atmosphere, although measurement difficulties have sparked a debate on these results 33,38,43 . Titan's subsurface ocean is predicted to be a mixture of water and ammonia, which can erupt from the surface via cryovolcanism, outgassing, and geysering events that move ammonia from Titan's interior to its surface 33,[44][45][46][47] . The co-existence of acetylene with ammonia, either in Titan's atmosphere or via deposition of ammonia-rich slurry on top of solid acetylene rich deposits, can quickly lead to the formation of a cocrystal that is stable under Titan's surface conditions and to methane/ethane fluvial and pluvial events 33 . As a consequence, the acetylene:ammonia (1:1) co-crystal is likely to be found near geochemically and biologically interesting sites such as cryovolcanic cones, ammonia eruption pits, and surface flows, in addition to Titan's stratosphere. 15,33,48,49 Here, we explore the structure and dynamics of the acetylene:ammonia (1:1) co-crystal using ab initio density functional theory-based molecular dynamics simulations at T = 30 K and T = 90 K, where the temperatures correspond to a crystal and plastic crystal, respectively. To identify static signatures of the plastic crystal phase, we quantify the translational and orientational structure of the co-crystal and identify the orientational disorder in the plastic crystal phase. Through an examination of translational and rotational dynamics, we show that the acetylene:ammonia (1:1) co-crystal exhibits rapid dynamic orientational disorder in the plastic phase. We connect the observed orientational disorder to the dynamics of N-H· · · π type hydrogen bonding within the co-crystal. We end with a discussion of the vibrational fingerprints of co-crystal formation while also reflecting on the accuracy of our results through comparison of experimentally-measured and computed vibrational spectra. We conclude by highlighting the commonality of plastic phases in Titan's potential cryominerals and put our results within the broader context of developing a thorough understanding of Titan's minerals.

SIMULATION DETAILS
We performed Born-Oppenheimer molecular dynamics (BOMD) simulations using density functional theory (DFT) within the QUICKSTEP 50 electronic structure module of the CP2K code [51][52][53] . QUICKSTEP employs a dual atom-centered Gaussian and plane wave (GPW) 50,54 basis approach for representing wavefunctions and electron density, leading to an efficient and accurate implementation of DFT. We used the molecularly optimized (MOLOPT) Goedecker-Teter-Hutter (GTH) triple-ζ single polarization (TZVP-MOLOPT-GTH) Gaussian basis set 55 for expanding orbital functions, along with a plane wave basis set with a cutoff of 500 Ry for representing the electron density. The core electrons were represented using GTH pseudopotentials [56][57][58] . Exchange-correlation (XC) interactions were approximated using the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) to the exchange-correlation functional 59 , as implemented in CP2K. To account for long-ranged dispersion interactions, we used Grimme's D3 van der Waals correction (PBE+D3) 60,61 . To assess the accuracy of the density functional approximation, we also performed simulations using the r 2 SCAN meta-GGA 62 combined with the appropriate non-local rVV10 correction for the longrange dispersion effects (r 2 SCAN+rVV10) 63,64 . We determined a planewave cutoff of 700 Ry to be sufficient for describing the system with r 2 SCAN+rVV10.
We modeled a 2 × 2 × 2 supercell of the acetylene:ammonia (1:1) co-crystal using the structure reported by Boese et al. (Cambridge Structural Database ID: FOZHOS) 65 . We equilibrated the system for at least 10 ps in the canonical (NVT) ensemble at constant temperatures of 30 K and 90 K using the canonical velocity rescaling (CSVR) thermostat 66 . For the calculation of dynamic properties, equations of motion were propagated in the microcanonical (NVE) ensemble for approximately 50 ps using a velocity Verlet integrator with a timestep of 1 fs. Maximally localized Wannier functions (MLWFs) were obtained on-the-fly, minimizing the spreads of the MLWFs according to the general and efficient formulation implemented in CP2K [67][68][69] . The centers of the computed MLWFs were used to define molecular dipoles, and polarizability tensors which were then used to extract IR and Raman signatures from the generated BOMD trajectories using the TRAVIS analyzer [70][71][72] . The IR and the Raman spectra were computed from a 20 ps long trajectory with MLWFs obtained at every step. The matrix elements of the full Raman polarizability tensor were extracted from the displacements of the MLWF centers recomputed by applying an electric field of strength 0.0005 a.u. along the crystal x, y, and z directions respectively. Each simulation with an electric field also resulted in production runs of 20 ps, such that the Raman spectrum is estimated from a total of 80 ps of simulation (four 20 ps trajectories).
The IR spectra were also computed for pure solid acetylene and solid ammonia. Solid acetylene undergoes a phase transition into an orthorhombic phase at 133 K. However, to the best of our knowledge, the crystal structure of orthorhombic solid acetylene has not been determined. Instead, we used the crystal structure for the orthorhombic phase of dideuteroacetylene (C 2 D 2 ) (Cambridge Structural Database ID: ACETYL07) as a FIG. 1. Snapshot illustrating the zig-zag geometry of the acetylene:ammonia (1:1) co-crystal with C-H· · · N and N-H· · · π type hydrogen bonding interactions shown as dotted lines. The hydrogen bonds are colored according to the donor atom. Views from both (a) crystal z and (b) crystal y axes are shown.
close approximation, and we replaced the deuterium atoms with hydrogens in our simulations 73 . For solid ammonia, we similarly used the crystal structure for cubic trideuteroammonia (ND 3 ) (Cambridge Structural Database ID: 34244) 74 .

RESULTS AND DISCUSSIONS Translational Order and Orientational Disorder at Titan Conditions
The acetylene:ammonia (1:1) co-crystal has a layered structure consisting of antiparallel planes of zigzag chains, as shown in Figure 1 15,33,65,75 . The layered arrangement within the co-crystal is primarily held together via a network of hydrogen bonds and weak van der Waals interactions. We find two types of hydrogen bonds in the acetylene:ammonia (1:1) co-crystal: C-H· · · N and N-H· · · π hydrogen bonds 33,65,75 , shown as blue and red dashed cylinders in Fig. 1, respectively. The C-H· · · N hydrogen bonds are formed by electrostatic attractions between the lone pair electrons of the nitrogen atom and the H atoms of acetylene molecules. Similarly, electrostatic attractions between the H atoms of ammonia and the electron-rich π-system of a neighboring acetylene molecule give rise to N-H· · · π hydrogen bonds. These C-H· · · N and N-H· · · π interactions play a crucial role in determining the structure and dynamics of the acetylene:ammonia (1:1) co-crystal.
We quantify the structure of the co-crystal by computing site-site radial distribution functions (RDFs), g XY (r) 76,77 , where X and Y refer to atomic sites or the midpoint of the C-C bond when representing the acetylene π electrons. The periodic arrangement of carbon and nitrogen atoms within the co-crystal is reflected in the periodic nature of C-C and N-N RDFs (Figure 2). Integration of g CC (r) and g N N (r) over the first peak yields coordination numbers of approximately 4 and 6, respectively, consistent the arrangement of atoms within the co-crystal geometry. Increasing temperature from 30 K a b to 90 K results in broader, lower intensity peaks in the RDFs, consistent with increased motion in the crystal due to thermal excitation. These results are largely independent of density functional approximation, although the r 2 -SCAN+rVV10 functional gives slightly larger and narrower peaks in the C-C and N-N RDFs than those predicted by PBE+D3.
The hydrogen bonds between acetylene and ammonia molecules are reflected in the various RDFs quantifying correlations between hydrogen bonding groups, Fig. 3. The C-H· · · N hydrogen bonds consist of acetylene donating a hydrogen bond to ammonia. Therefore, these hydrogen bonds can be quantified through C-N and (C)H-N RDFs, where (C)H is used to indicate that the hydrogen is bonded to a carbon. The C-N RDF, g CN (r), exhibits a large first peak near 3.5Å (Figure 3a), the integral of which gives a coordination number of two. This peak corresponds to carbon atoms of two different acetylene molecules that are H-bonding with ammonia, highlighted by the snapshots in Figure 1. There is a second and third peak soon after (before 5Å), both of which correspond to carbons that are not involved in H-bonds with the central ammonia molecule. This structure is further emphasized by the (C)H-N RDFs, which exhibit a peak at approximately 2.25Å, suggesting that acetylene hydrogens point directly at the lone pair of the nitrogen to form C-H· · · N hydrogen bonds (Fig. 3c). Integration of the first peak in this RDF yields a coordination number of two, further suggesting that the nitrogen atom is, on average, engaged in hydrogen bonding with two acetylene molecules.
To characterize N-H· · · π hydrogen bonds, we need a site to represent the π system of each acetylene molecule. Here, we represent the location of the π system as a single site located at the midpoint of the C-C triple bond, but note that similar results are found if the centers of maximally localized Wannier functions are used to represent the location of π electrons. The π-N RDF displays a single, well-defined peak just before 4Å, suggesting a single H-bond donated from an ammonia molecule to an acetylene molecule, Fig. 3b, as shown in the snapshots in Fig. 1. The π-H(N) RDF displays three peaks below 4.5Å, corresponding to the three hydrogens of the ammonia molecule, Fig. 3d. The first peak located at distances smaller than 3Å results from H atoms involved in a hydrogen bond with acetylene.
Increasing temperature from 30 K to 90 K weakens hydrogen bonds, as evidenced by the broader and less intense peaks in the RDFs. Importantly, heating the acetylene:ammonia (1:1) co-crystal to 90 K significantly increases the disorder observed in the π-H(N) RDF, Fig. 3d. The deep minima between peaks, especially that between the first and second peak, increase to be comparable in height to some of the peaks. This indicates that ammonia H atoms have a significantly higher probability of occupying these interstitial regions, suggesting that the ammonia molecules are orientationally disordered at 90 K.
We examine the dependence of the H-bond structure on density functional approximation at 90 K. We find that the first peak in the (C)H-N RDF is significantly larger in the r 2 -SCAN+rVV10 system than in the PBE+D3 one, such that the peak height is larger than that at 30 K (with PBE+D3). This suggests that C-H· · · N hydrogen bonds are stronger in the r 2 -SCAN+rVV10 system. In contrast, the π-H(N) RDF is more disordered in the r 2 -SCAN+rVV10 system, suggesting that N-H· · · π hydrogen bonds are weaker and the ammonia molecules are more disordered than in the PBE+D3 system. The weaker N-H· · · π hydrogen bonds produced by the r 2 -SCAN+rVV10 functional may be due to increased charge transfer along the C-H· · · N interaction suggested by the stronger C-H· · · N hydrogen bonds.
The RDFs discussed above are averaged over orientational degrees of freedom and, therefore, cannot directly quantify orientational order in the acetylene:ammonia (1:1) co-crystal. To quantify the orientational structure of ammonia molecules within the co-crystal, we compute the probability distribution of ammonia H atom positions projected on the xz plane of the crystal 78 . At 30 K, the ammonia molecules are orientationally ordered, as reflected by the distinct, disconnected regions of high probability, Fig. 4a. The six distinct peaks in the probability distribution are consistent with the two types of configurations of ammonia molecules shown in Fig. 1a, which are orientationally ordered but point their dipole moments (and hydrogens) in opposite directions.
In contrast, the ammonia molecules show evidence of orientational disorder at 90 K; the probability distribution is continuous along an ellipse in the xz plane, Fig. 4b. The appearance of non-zero probability between the peaks in the distribution at 90 K indicates that ammonia molecules can adopt all orientations in the xz plane, although the six peaks are still the most populated orientations. Similar results are obtained at 90 K using the r 2 -SCAN+rVV10 functional, although the distribution is broader because the ammonia molecules are more disordered, Fig. 4c H· · · π hydrogen bonds in this system. Therefore, in the temperature range on Titan's surface, our results suggest that the acetylene:ammonia (1:1) co-crystal is in an orientationally disordered phase -the acetylene:ammonia (1:1) plastic co-crystal.

Translational and Orientational Dynamics
The structural analysis above gives no information about molecular dynamics within the co-crystal, and therefore cannot suggest whether or not the orientational disorder is dynamic. Therefore, we also quantify the translational and orientational dynamics of molecules in the acetylene:ammonia (1:1) co-crystal. We first quantify the translational dynamics of the co-crystal by computing normalized velocity auto-correlation functions, where v i (t) is the velocity of atom i at time t and · · · denotes an ensemble average. The power spectrum, C v (ω), is obtained by Fourier transforming C v (t), and corresponds to the vibrational density of states of the system, functionals used here show similar behavior, with slight differences in peak shapes. The ω → 0 limit of C v (ω) is proportional to the diffusion coefficient 79 . Importantly, all three power spectra vanish in this limit, indicating a lack of translational diffusion within the crystal, as expected for translationally ordered crystals.
We probe the orientational dynamics within the cocrystal by computing the rotational time correlation functions C n (t) = P n (n(t) · n(0)) , where n(t) is an N-H bond vector for ammonia and a C-C bond vector for acetylene molecules at time t and P n (x) is the nth order Legendre polynomial. We focus on C 1 (t) and C 2 (t), which can be determined through spectroscopy 80,81 . Correlations in the orientational structure of acetylene molecules persist throughout the duration of the simulations at both 30 K and 90 K, Fig. 6a,c, suggesting that the acetylene molecules within the cocrystal are in a static crystalline phase and do not un-dergo any significant orientational dynamics. Similar long-time correlations occur for ammonia molecules at 30 K, Fig. 6b,d, indicating an absence of any rotational dynamics within the co-crystal at low temperatures. In contrast, C 1 (t) and C 2 (t) for ammonia molecules decay on picosecond timescales at 90 K. This indicates that ammonia molecules exhibit dynamic orientational disorder. We estimate the rotational correlation times τ 1 and τ 2 as the integral of the corresponding time correlation functions, C 1 (t) and C 2 (t), respectively. In practice, we split the correlation function into two parts at t = τ cut to perform the integration: (i) the short time part from t = 0 to t = τ cut , which we numerically integrate (shown for C 2 (t) in the inset of Fig. 6d), and (ii) the rest of the correlation function, which we fit to an exponential function, A 0 e −τ fit /t +A 1 . Then, the rotational correlation time is given by Using τ cut = 0.5 ps, we obtain rotational correlation times of τ 1 = 1.6 ± 0.2 ps for PBE+D3 and τ 1 = 0.77 ± 0.01 ps for r 2 -SCAN+rVV10, and τ 2 = 1.5 ± 0.1 ps and τ 2 = 0.58 ± 0.01 ps for the PBE+D3 and r 2 -SCAN+rVV10 systems, respectively. This fast decay of ammonia rotational correlations suggests that the acetylene:ammonia (1:1) co-crystal is in a dynamic orientationally disordered plastic phase at Titan surface conditions.
Because the ammonia molecules rotate, but the acetylene do not, ammonia rotations are coupled to transient breakage and reformation of N-H· · · π hydrogen bonds. To quantify the dynamics of N-H· · · π hydrogen bonds, we compute the hydrogen bonding correlation function [82][83][84]  where the indicator function is given by and Θ(x) is the Heaviside function. With appropriate geometric definitions for a hydrogen bond, h(t) = 1 if a hydrogen bond exists between a triplet of sites at time t and h(t) = 0 otherwise. For N-H· · · π hydrogen bonds, R(t) and φ(t) correspond to the N· · · π distance and the angle formed by the N-H and N· · · π vectors, respectively; the π-site is taken to be the midpoint of the C-C bond.
To define a N-H· · · π hydrogen bond, we use a radial distance cutoff of R c = 4Å and an angular cutoff of φ c = 30 • . For comparison, we also compute the hydrogen bond correlation function for C-H· · · N hydrogen bonds, where R(t) corresponds to a C-N distance, φ(t) is the angle made by the C-N and (C)H-N vectors, R c = 4Å, and φ c = 30 • . At 30 K, C HB (t) does not decay for either type of hydrogen bond, Fig. 7a,b. The lack of decay indicates that there is no significant hydrogen bond breakage within the co-crystal at low temperatures. At 90 K, C HB (t) also does not significantly decay for C-H· · · N hydrogen bonds, Fig. 7a, indicating that these interactions remain intact for long periods of time. In contrast, C HB (t) for N-H· · · π hydrogen bonds at 90 K decay rapidly, suggesting that hydrogen bond breakage accompanies the rotational motion of ammonia molecules. We determine the hydrogen bond correlations time, τ HB , from C HB (t) following the same approach used for rotational times. We find correlation times for N-H· · · π hydrogen bond relaxation of τ HB = 0.9 ± 0.06 ps from PBE+D3 and τ HB = 0.35 ± 0.01 ps from r 2 -SCAN+rVV10, slightly faster than the corresponding rotational times. This must be the case if hydrogen bond breakage accompanies rotation. Moreover, the faster hydrogen bond dynamics of r 2 -SCAN+rVV10 are consistent with weaker N-H· · · π hydrogen bonds in this system. These faster hydrogen bond dynamics are also consistent with the faster rotational dynamics produced by r 2 -SCAN+rVV10, further suggesting that hydrogen bond breakage accompanies rotational dynamics.

Theoretical Vibrational Spectroscopy
Infrared (IR) and Raman spectroscopy are complimentary vibrational techniques that are well-suited for studying Titan cryominerals because of the sensitivity of molecular vibrations to their environment 15,33,[85][86][87][88] . To provide molecular-scale interpretation of experimental spectra and to validate the accuracy of our models through comparison with experiments, we computed both the IR and Raman spectra of the acetylene:ammonia (1:1) co-crystal. The IR absorption spectra were computed from the Fourier transform of the dipole derivative autocorrelation function defined according to 70 , where µ(t) refers to the instantaneous dipole moment andμ(t) is its time derivative. The unpolarized Raman where k B is Boltzmann's constant, ω in is the frequency of the incident light, is Plank's constant, and T is the temperature. The isotropic (a p ) and anisotropic polariz-abilities (γ p ) are defined according to and where the C ab (t) indicates the time correlation function of the time derivatives of the polarizability, α ab (t), and a and b indicate x, y, z.
We assess the accuracy of our IR results by comparing them to the experimentally measured band shifts and associated frequencies in the pure and co-crystalline phase 75 . The IR spectra computed for the acetylene:ammonia (1:1) co-crystal and the pure acetylene and ammonia crystals agree well with the corresponding experimental spectra, Fig. 8. In addition to the absolute peak positions, the peak shifts that occur upon co-crystallization are in reasonably good agreement with the experiment, Table I. This good agreement between theory and experiment suggests that the density functional approximation used here can describe the acetylene:ammonia (1:1) co-crystal with reasonable accuracy and yield accurate predictions for quantities like rotational correlation times.
The blue shift of the C 2 H 2 ν 5 asymmetric bend and the red shift of the ν 3 asymmetric C-H stretch upon crystallization have been identified as characteristic signatures of co-crystal formation 75 , and these shifts are captured by our predicted IR spectra. The blue shifts of the cocrystal ν 5 bands are slightly stronger in our simulation as compared to the observed experimental frequencies.
Analogously, the red shift of the C 2 H 2 ν 3 asymmetric stretch is more pronounced in our simulations than in experiments because of the large blue shift of the ν 3 vibration in the pure phase with respect to the experimental ν 3 frequency. Our computed spectrum also captures the relative shifts of the NH 3 bands, however, given their small magnitude evident from the experimental shifts, they are rather subtle to observe. Notably, the blue shift in the NH 3 ν 3 vibrational band is accurately reproduced in our spectrum though the relative shift in frequency is overestimated.
In addition to IR, we also compute the Raman spectrum for the acetylene:ammonia (1:1) co-crystal at 90 K, Fig. 9. We compare the computed Raman signatures with the experimental bands reported by Cable et al. 33 . The peak positions predicted by our simulations are in good agreement with, but slightly higher frequency than those obtained experimentally, Table II. This level of agreement further validates the accuracy of our models and suggests that the range of predictions made regarding rotational timescales is reasonable. FIG. 9. Unpolarized Raman spectra computed for the acetylene:ammonia (1:1) co-crystal at T = 90 K. Vertical lines indicate experimentally measured peak positions 33 .

CONCLUSIONS
We have used ab initio molecular dynamics simulations to characterize the structure and dynamics of the acetylene:ammonia (1:1) co-crystal. At 90 K, the temperature of Titan's surface, the acetylene:ammonia (1:1) cocrystal is a plastic co-crystal. Ammonia molecules within the co-crystal are orientationally disordered and undergo rotational motion on a picosecond timescale. This rotational correlation time should be accessible in experiments through Raman and NMR spectroscopies 80,81,89 . At 30 K, the acetylene:ammonia (1:1) co-crystal is in a crystalline phase (translationally and orientationally ordered), such that a co-crystal-to-plastic co-crystal phase transition occurs between 30 and 90 K. These conclusions agree well with the early experimental hints at the orientational disorder of ammonia molecules within the co-crystalline phase 65 . The rotational dynamics of the ammonia molecules are dictated by the breakage and reformation of N-H· · · π hydrogen bonds within the cocrystal, which we predict to have a lifetime on the order of a picosecond. The accuracy of our predictions is sup-ported by the good agreement between experimentally measured and predicted IR and Raman spectra for the co-crystal.
We expect the acetylene:ammonia (1:1) co-crystal to be one of many Titan cryominerals that exist in semidisordered phases between crystal and liquid. Many of the multicomponent cryominerals that have been predicted are made of molecular components that exhibit orientationally disordered plastic crystal phases 15,16 , suggesting that these co-crystals might also exhibit these disordered phases. Understanding the existence of these disordered phases at Titan surface conditions is an essential step in understanding the geology of Titan from the molecular scale. For example, the dynamic disorder of plastic crystals can lead to significantly different mechanical properties than those of the ordered crystalline phase. In this context, the disorder of these phases must be accounted for when predicting mechanical properties, especially for large-scale modeling of surface geology.
We note that our ab initio simulations do not account for nuclear quantum effects. The quantum mechanical nature of nuclei, especially light nuclei, becomes increasingly important at low temperatures like those on the surface of Titan. Future work will detail the importance of nuclear quantum effects in the molecular structure and dynamics of Titan's cryominerals.

ACKNOWLEDGMENTS
This work is supported by the National Aeronautics and Space Administration under grant number 80NSSC20K0609, issued through the NASA Exobiology Program. We acknowledge the Office of Advanced Research Computing (OARC) at Rutgers, The State University of New Jersey for providing access to the Amarel cluster and associated research computing resources that have contributed to the results reported here. This work used the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS), formerly Extreme Science and Engineering Discovery Environment (XSEDE) 90 , which is supported by National Science Foundation grant number ACI-1548562. Specifically, this work used Stampede2 and Ranch at the Texas Advanced Computing Center through allocation TG-CHE210081.