On the Occurrence of Clathrate Hydrates in Extreme Conditions: Dissociation Pressures and Occupancies at Cryogenic Temperatures with Application to Planetary Systems

We investigate the thermodynamic stability of clathrate hydrates at cryogenic temperatures from the 0 K limit to 200 K in a wide range of pressures, covering the thermodynamic conditions of interstellar space and the surface of the hydrosphere in satellites. Our evaluation of the phase behaviors is performed by setting up quantum partition functions with variable pressures on the basis of a rigorous statistical mechanics theory that requires only the intermolecular interactions as input. Noble gases, hydrocarbons, nitrogen, and oxygen are chosen as the guest species, which are key components of the volatiles in such satellites. We explore the hydrate/water two-phase boundary of those clathrate hydrates in water-rich conditions and the hydrate/guest two-phase boundary in guest-rich conditions, either of which occurs on the surface or subsurface of icy satellites. The obtained phase diagrams indicate that clathrate hydrates can be in equilibrium with either water or the guest species over a wide range far distant from the three-phase coexistence condition and that the stable pressure zone of each clathrate hydrate expands significantly on intense cooling. The implication of our findings for the stable form of water in Titan is that water on the surface exists only as clathrate hydrate with the atmosphere down to a shallow region of the crust, but clathrate hydrate in the remaining part of the crust can coexist with water ice. This is in sharp contrast to the surfaces of Europa and Ganymede, where the thin oxygen air coexists exclusively with pure ice.


Introduction
Clathrate hydrates are nonstoichiometric compounds composed of guest molecules accommodated in cages of the host lattice structure (Davidson 1973;Sloan & Koh 2008). The cages are limited mostly to 12-, 14-, and 16-hedron when accommodating only one sort of guest species. Those are known as CS-I, made of 12-and 14-hedra, and CS-II, made of 12-and 16-hedra. Guest molecules are generally apolar or only weakly polarized, such as hydrocarbons and carbon monoxide. Clathrate hydrates are stable only when encaging a number of guest molecules, which depends on temperature and pressure as well as guest species.
Because clathrate hydrates of simple molecules attract considerable attention owing to their existence on Mars, some satellites of the outer planets, and also comets, the thermodynamic condition for their formation/dissociation has been examined, focusing on the existence of clathrate hydrates within the interiors (Lewis 1971;Lunine & Stevenson 1985;Choukroun et al. 2010;Mousis et al. 2010Mousis et al. , 2014Kite et al. 2017). In particular, they are believed to exist in large amounts in the crust of dwarf planets and the satellites of outer planets such as Pluto and Titan (Tobie et al. 2006;Bouquet et al. 2019;Kamata et al. 2019). The existence of clathrate hydrates on Titan is also argued in conjunction with the observation by the Cassini Radar (Lopes 2019). Methane clathrate hydrate can be associated with the outgassing mechanism of methane into the atmosphere of Titan by cryovolcanism (Davies et al. 2016) and with its insulating effect (Kalousová & Sotin 2020). It has also been pointed out that the dissociation of clathrate hydrate rather than the phase transition of amorphous ice is responsible for the observed outgassing from a comet (Luspay-Kuti 2016). The stable regions of clathrate hydrates in the form of CS-I and CS-II are widely spread in temperature (T)-pressure (p) space, roughly T=0 to 300 K and p=10 −12 to 10 8 Pa. Owing to the important role of clathrate hydrates in planetary science, much effort should be devoted to measuring and/or estimating their stability. The formation and dissociation of clathrate hydrate are expected to be very slow at cryogenic temperatures. This makes it difficult to observe thermodynamic conditions associated with the formation and dissociation processes in the low-temperature regime in laboratory experiments while equilibrium can be attained on a geological timescale. It seems desirable to establish a way to estimate the phase equilibrium of clathrate hydrate to cover the thermodynamic conditions mentioned above on the basis of the statistical mechanics available in the low-temperature regime.
A theoretical method was proposed to calculate the water/ hydrate/guest three-phase equilibrium by van der Waals & Platteeuw (1959, hereafter vdWP). The vdWP theory was formulated for clathrate hydrate in which one of the constituents, the host lattice, is made of a fixed number of water molecules, and guest molecules enter or leave the lattice under the condition in which the chemical potential of the guest is the same as that of the coexisting guest fluid. The other equilibrium between the hydrate and water phases is realized by the equivalence of the chemical potential of water in the hydrate phase to that in the aqueous (ice) one. The three-phase coexistence can be estimated from reliable intermolecular interactions with various revisions, resulting in a better agreement with experiments (Tanaka & Kiyohara 1993a, 1993b. Advantages of the vdWP theory have been detailed elsewhere (van der Waals & Platteeuw 1959;Sloan & Koh 2008;Tanaka & Kiyohara 1993a, 1993bTanaka & Matsumoto 2013).
For a complete description of the phase behaviors of clathrate hydrate, the essential knowledge is the locations of the phase boundaries between water (liquid or crystalline ice) and clathrate hydrate, referred to as water/hydrate, and that between clathrate hydrate and the guest fluid (or solid), referred to as hydrate/guest, in T, p, and composition (y) space. The phase behaviors of clathrate hydrates have been examined in the entire p-T-y space by introducing empirical parameters or developing a theory to calculate the chemical potentials of both guest species and host water in clathrate hydrate (Wierzchowski & Monson 2007;Tanaka et al. 2018aTanaka et al. , 2018b. This enables us to evaluate the boundaries in terms of intermolecular interactions. It has been found that the shape of the two-phase boundaries in the phase diagram exhibits a variety of guest dependence as shown in Figure 1. This figure demonstrates that the phase boundaries of CH 4 on the T-y plane are dissimilar to those of C 2 H 6 despite the fact that both guest species form CS-I hydrate. Specifically, the areas of the stable zone enveloped by the two boundaries are remarkably different from each other, and the content of C 2 H 6 varies significantly from one boundary to the other because of the large difference between the affinity of C 2 H 6 in a small cage and that in a large cage, which is in turn due to the size of the C 2 H 6 molecule.
Measurements of the equilibrium properties become difficult with decreasing temperature. The difficulty is further compounded close to 0 K by the crystallization of guest species, which is instead a fluid state in most laboratory experiments. So far, the revised vdWP theory has been applied to clathrate hydrates in thermodynamic conditions accessible in laboratory experiments. In this study, we reveal that the method is applicable in low temperatures by substituting the quantum partition functions for the classical ones and by considering crystalline solids of guests in equilibrium with clathrate hydrates. With the aid of this implementation, it is possible to examine the phase behaviors of various clathrate hydrates in extreme conditions from only intermolecular interactions. To this end, we first examine the thermodynamic stability of clathrate hydrates at 0 K and then show that the stability zone is smoothly connected to the finite temperature region, avoiding an extrapolation from high temperatures.
We calculate the thermodynamic properties of clathrate hydrates of small hydrocarbons and noble gases. Small hydrocarbons are key substances in exploring the interior activity of icy satellites of outer planets such as the outgassing mechanism of CH 4 on Titan (Tobie et al. 2006). Noble gases, specifically the abundance of individual isotopes, are important species in relation to the various stages of accretion of planetary systems. Those substances form clathrate hydrates in the presence of water in water-rich or guest-rich conditions. For example, the former corresponds to the condition on the subsurface of the internal ocean, and the latter is expected on the surface of Titan. It should be emphasized that either equilibrium involving a clathrate hydrate can occur at thermodynamic conditions even far away from the three-phase coexistence. Such an equilibrium is deduced from the phase diagram in the whole thermodynamic space, as suggested in Figure 1.
Because an equilibrium state is expected to be realized on a geological timescale at low temperatures, it is crucial to assess the equilibrium thermodynamic conditions for the formation and dissociation of clathrate hydrates theoretically. An elaborate calculation of the stability of clathrate hydrates is helpful to argue the role of clathrates in the outgassing and inventory of volatiles. This paper is organized as follows. The computational methods are presented in Section 2. Briefly described in Section 2.2 is the revised vdWP theory for the evaluation of the phase behaviors of clathrate hydrates in the entire space of thermodynamic variables. A method to calculate the thermodynamic properties derived from the free energy of clathrate hydrates is also presented. We show its application to various clathrate hydrates to determine the phase diagrams in the lowtemperature and low-pressure region including the conditions on the surfaces of Titan, Europa, and Ganymede in Section 3. We highlight the role of clathrate hydrate on Titan as an example by noting the phase behaviors under cryogenic conditions and by indicating the importance of considering not only the three-phase coexistence but the two-phase coexistence, which occurs exclusively in real systems. The latter has bearing on our inference on the likely scenario in the crust of Titan involving methane clathrate hydrate. Our study is concluded in Section 4 with a few remarks.

Generation of Crystalline CS-I, CS-II, and Ice Ih Structures
We generate 500 hydrogen-disordered configurations for each of the CS-I, CS-II, and low-pressure ice (ice Ih) structures using the GenIce tool to calculate the free energy and other properties (Matsumoto et al. 2018). The CS-I, CS-II, and ice Ih structures consist of 384, 1088, and 1120 water molecules, respectively. Each generated configuration has no net polarization and satisfies the ice rule (Bernal & Fowler 1933;Figure 1. Phase boundaries of CH 4 and C 2 H 6 hydrates in the temperaturecomposition plane at 10 MPa. The solid and dotted lines indicate the water/ hydrate and the hydrate/guest boundaries, respectively, at temperatures higher than 200 K. The water/hydrate/guest three-phase equilibrium is achieved at the intersections between the solid and dotted lines, marked by filled circles. Note that only the loci of the intersections projected onto the p-T plane, i.e., the temperature dependence of the dissociation pressure, are obtained from the original vdWP theory. Pauling 1935). The number of generated structures is larger than in previous works, and we adjust as explained below the water-guest interaction according to the newly generated structures, though the magnitude of the change is minor. The standard deviation in the chemical potential arising from hydrogen disordering is roughly 0.01 kJ mol −1 .

Intermolecular Interaction Potentials
The clathrate hydrates investigated here are composed of water and apolar guest molecules. The water-water interaction is described by the TIP4P/Ice model, which satisfactorily reproduces the melting temperature of ice Ih (Abascal et al. 2005;Espinosa et al. 2014). A TIP4P/Ice water molecule consists of three charged sites corresponding to an oxygen (a little distant from it) and two hydrogen atoms and one Lennard-Jones (LJ) interaction site located on the oxygen atom. The O-H distance and H-O-H angle are fixed to 0.095 72 nm and 104.52°.
Small hydrocarbons, tetrafluoromethane, oxygen, nitrogen, and noble gases are chosen as the guest species, each of which is approximated by a spherical molecule represented by a single LJ interaction site except for the linear C 2 H 6 model. This approximation is rationalized by the good agreement of the calculated dissociation pressures with experimental ones (Tanaka et al. 2018a(Tanaka et al. , 2018b). An effect of the deviation from the sphere is discussed for C 2 H 6 hydrate. The LJ size parameters are taken from a classical textbook (Hirschfelder et al. 1954) and the Optimized Potentials for Liquid Simulations (OPLS) force field model (Jorgensen et al. 1984). The energy parameter for clathrate hydrates is optimized roughly to reproduce the experimental dissociation pressure around 273 K by multiplying a scaling factor, s, to the original energy parameter, ε 0 . Those parameters are tabulated in Table 1. The Lorentz-Berthelot rules are applied to interactions for any unlike pairs. 3. Theory for the Stability of Clathrate Hydrates in the Low-temperature Regime

Condition for Two-phase Coexistence
We consider two two-phase boundaries, one between water and clathrate hydrate, the water/hydrate boundary, and the other between clathrate hydrate and the guest species, the hydrate/guest boundary. Those boundaries form curved surfaces in temperature, pressure, and composition space. Here, we restrict ourselves to a simple clathrate hydrate containing one guest species in the two-phase coexistence condition having only one component axis (mole fraction of guest species, y), except for the N 2 +CH 4 mixed hydrate expected on the surface of Titan. Thus, we treat mostly equilibria between water and clathrate hydrate and between clathrate hydrate and the crystalline solid or gaseous state of the guest species, which meets the cryogenic condition. Because the guest content in crystalline ice Ih is extremely low, the ice Ih phase can be well approximated by that of pure ice without guest impurity (Petrenko & Whitworth 1999). The water/ hydrate boundary emerges when the chemical potential of water in clathrate hydrate, μ c , equals that of ice Ih, μ i , (μ c = μ i ). Similarly, the chemical potential of the guest species can be set to the corresponding value of each pure guest in the solid or fluid phase denoted by m g 0 . Notably, it is a crystalline solid at 0 K at any finite pressure. The hydrate/guest coexistence imposes the condition m m = g g 0 , where μ g indicates the chemical potential of the guest species in clathrate hydrate.

Statistical Mechanical Theory for a Fixed Number of Guests
Clathrate hydrate is the only stable phase in the region enclosed by the two-phase boundaries in the phase diagram ( Figure 1). The boundaries are calculated from the chemical potentials of water and guest species on the basis of statistical mechanics theory.
We assume that the Helmholtz free energy of a clathrate hydrate is composed of three contributions, the free energy of the empty lattice structure, the free energy of guests in each cage, and the mixing entropy arising from the random distribution of occupied and unoccupied cages. Thus, the Helmholtz free energy of a clathrate hydrate, A c (T, V, N w , N g ), made of N w water and N g guest molecules at a fixed volume V and temperature T is constructed as indicates the Helmholtz free energy of an empty clathrate hydrate, and the sum is taken over all possible cage types with the ratio of the number of the cages of type j to that of water molecules denoted by α j and the corresponding occupancy by x j (0 x j 1). The free energy of the guest in a cage of type j is denoted by f j . The occupancy of each cage type is calculated by introducing an undetermined parameter, μ g , as Once the mole fraction, y, of the guest in clathrate hydrate is given, the parameter, μ g , corresponding to the chemical potential of the guest in clathrate hydrate is obtained by solving

Chemical Potentials of Water in Ice Ih and Clathrate Hydrate
The Helmholtz free energies of ice and empty clathrate hydrate consisting of N w water molecules at temperature T and volume V is composed of three terms, the average interaction energy at its minimum energy structure, , and the residual entropy, S r (N w ), arising from hydrogen disordering. Hence, the free energy is given as in which Pauling's approximate residual entropy (Pauling 1935), is adopted for ice Ih, empty CS-I, and empty CS-II.
For small-amplitude vibrations expected in the low-temperature regime, the vibrational free energy is approximated as a collection of harmonic oscillators for ice Ih. It is calculated from the quantum mechanical partition function as (Pohorille et al. 1987) The 500 hydrogen-disordered ice Ih structures are optimized using the steepest descent algorithm. The potential energy in Equation (4) , is the average over the 500 optimized energies. The Hessian matrix is constructed for each optimized structure. The frequency of intermolecular vibrational motion in Equation (6), ν j , is obtained by diagonalization of the Hessian matrix.
The equilibrium volume á ñ V for a given pressure, p, satisfies the following condition, The Gibbs free energy is obtained using á ñ V as

w
The chemical potential of water in ice Ih is calculated (Tanaka 1998) as The Helmholtz free energy of empty hydrate, is similarly decomposed, where the vibrational part of the free energy of the empty hydrate, , is also replaced by that of the harmonic vibrations as given on the right-hand side of Equation (6). The free energy of clathrate hydrate in Equation (1) is rewritten as as a function of pressure (Yagasaki et al. 2016). The chemical potential of water in clathrate hydrate is expressed as a function of T, p, and y as is the chemical potential of the empty clathrate hydrate.

Free Energy of Cage Occupation
Because the coupling between host and guest motions is weak for small guest molecules such as Ne, Ar, and O 2 , the frequency modulation of the host lattice due to the presence of guests is insignificant (Tanaka & Kiyohara 1993b). Then, it is reasonable to assume that the contributions from the host and guest to the free energy are separable, which is the basic idea in the original vdWP theory to estimate the dissociation pressures (van der Waals & Platteeuw 1959). In the opposite case where the guest size is comparable to that of the cage, such as Xe, C 2 H 6 , and CF 4 , a harmonic approximation for the motions of guest molecules works well, where the frequency modulation due to the presence of guests is not negligible (Tanaka & Kiyohara 1993a). We describe below two methods for calculating the free energy of the guest inside each cage (for a guest of an intermediate size such as Kr, N 2 , and CH 4 , the original vdWP theory is still applicable, but the frequency modulation of the host lattice should be taken into account).
The potential surface of a small guest molecule inside a cage, represented against the radial distance, deviates significantly from parabolic and is characterized by a minimum at an offcenter position in each cage (Yagasaki et al. 2016). In order to take account of such a large anharmonic character, the free energy of cage occupation for a small molecule is calculated from the eigenenergies of the guest in the potential field of the surrounding host water molecules. The energy, E j , is actually evaluated by assuming that the potential is spherically symmetric. Thus, it is orientationally averaged and is a function of only the radial distance from each cage center. The Schrödinger equation with this potential is solved numerically for the radial distance while the orientational part is simply described by the spherical harmonics. The obtained eigenenergy, E j , is used to calculate the free energy of cage occupation, f k , in a cage of type k for a small guest molecule considering the degeneracy, Ω j , as This method is adopted for Ne, Ar, Kr, O 2 , N 2 , CH 4 , and Xe. The potential surface of a large guest is approximately parabolic, and its motion is well described by harmonic oscillation. Therefore, the free energy for a large molecule is calculated as the sum of the minimum interaction energy of the guest at the center of the cage and the difference in the vibrational free energy between the occupied and empty hydrate in the same way as Equation (4). That is, the vibrational free energy of the occupied hydrate, , is calculated from a set of harmonic vibrational frequencies of the guest molecules coupled with the host lattice as For a clathrate hydrate having more than one cage type, we must consider the filling order of the cage type, which is easily anticipated by the magnitude of the interaction energy. That is, the free energy of cage occupation in a type k cage is given as where ¢ N g indicates the number of all filled cages with the guests having lower free energy than f k , and N gk denotes the number of such cages . In this procedure, the free energy change of the host lattice due to the occupation of cages (by the frequency modulation) is incorporated in f k . This manipulation is applied to Xe, C 2 H 6 , and CF 4 .
We apply the two methods to Xe hydrate to make sure that a large guest species can be treated by either method. In addition, it is worthwhile to note that the free energy of cage occupation, f j , turns to a function of pressure in the process of conversion to the isobaric system through a procedure similar to Equation (7).

Chemical Potentials of the Guest in the Solid or Fluid Phase
The guest phase can be any state under the conditions we deal with (T 200 K and 10 −12 p 10 8 Pa), either a crystal or fluid. Each guest is approximated as a spherical particle, which naturally leads to the face-centered cubic (fcc) structure for any guest solid whatever the real crystalline structure in the low-temperature regime is, possibly a minor modification having lower symmetry. For the linear ethane model, the crystalline structure is the lowest temperature phase known as ethane III in which the orientations of individual molecules are also specified (van Nes & Vos 1978). The free energy of each guest crystal is calculated in the same way as ice Ih, eliminating the residual entropy term. The free energy in the fluid phase at a high pressure is calculated based on the corresponding Redlich-Kwong equation of state (Redlich & Kwong 1949). The free energy at a low temperature cannot be properly calculated with this equation of state. We use the Redlich-Kwong equation of state for the stability of CH 4 clathrate hydrate only in the high-temperature region above 200 K concerning the three-phase coexistence. The pressure along the three-phase coexistence curve turns out to be about four orders of magnitude lower than the vapor pressure line for liquid (or solid) CH 4 below 200 K where the vapor pressure is low enough to be treated as an ideal gas (Lide 2008). It is also sufficient to consider the guest-rich phase of either the crystalline or ideal gas state in all other calculations for the three-phase coexistence curve below 200 K.

Chemical Potentials on Phase Boundaries
The chemical potential of water in ice Ih given by Equation (9) is the same as that in clathrate hydrate given by Equation (11) when ice Ih coexists with clathrate hydrate. The coexistence conditions for the two phases are determined by The second equality is derived from Equation (2). In a practical application of Equation (15), a hybrid form of the two different expressions may be used when x i →1 and x k →0(i≠k). The boundary for hydrate/guest is obtained by setting where the left-and right-hand sides are the chemical potential of the guest in the clathrate hydrate and in the pure guest phase, respectively.

Summary of the Numerical Calculations
A program map (flowchart) of the actual calculations is given in Figure 2. We encode all of the programs in the figure.
The computational requirement is solely the size of the memory to perform the diagonalization of the mass-weighted Hessian.

Stability at 0 K
The structure of the clathrate hydrate we handle is restricted to either CS-I or CS-II. Thus, the number of cage types in each structure is two, with the larger cage (14-hedron for CS-I and 16-hedron for CS-II) suffixed by l and the small cage (12hedron) by s. The magnitude relationship between f l and f s relies on the guest species. We may expect f l <f s for large guest species such as C 2 H 6 and CF 4 , but the relation depends on the volume of clathrate hydrate for small guests such as CH 4 and Ar.
Here, we consider an extreme condition, T→0 K. The magnitude of the free energy of cage occupation relative to the chemical potential of the guest species places a limit on the composition value depending on the chemical potential of the guest in its pure state. In the case of m -< f 0 g j 0 , we can eliminate the occupation of cages of type j and consider only the occupation of the remaining cages. In this study, we will not treat this case, such as the C 3 H 8 hydrate, where guest molecules stay only in large cages. Because x j is zero for μ g −f j <0 and is unity for μ g −f j >0 in the limit of T→0 K, the mole fraction, y, is restricted to four values, 0, y l (all large cages occupied), y s (all small cages occupied), and y f (fully occupied) unless there is an accidental coincidence, μ g =f j , which has measure 0.
On the hydrate/guest boundary, the equilibrium condition is given by m m = .
g g 0 For most of the guest species, we have the relation that the free energy is lower than the chemical potential of the guest species in its pure guest solid phase (m -> f 0 g j 0 ) for any cage of type j and thus y=y f on the boundary.
On the water/hydrate boundary, the coexistence occurs at μ i =μ c . We, therefore, examine the chemical potential of water in clathrate hydrate for the following four cases.
(1) For f l >μ g and f s >μ g , no cages are occupied and the chemical potential is The chemical potential of water in ice Ih is always lower than that in the empty clathrate hydrate irrelevant to its structure (either CS-I or CS-II). Hence, a stable hydrate under this condition is precluded.
(2) For f s >μ g >f l , no small cages are occupied while all large cages are occupied, which leads to It is noted that the chemical potential of water increases through either the sequence of (1)→(2)→(4) or (1)→ (3)→(4). The occupancy of the cage is immediately obtained from the magnitude relationship between f j and μ g . The upper and lower bounds of the chemical potential of water in each case are calculated from Equations (17)- (20). Those values of μ c , cage occupancies x, and compositions y are summarized in Table 2 for possible ranges of the free energy of cage occupation, f j , relative to the chemical potential of the pure guest, m g 0 , and that in hydrate, μ g .
As the relationship between f l and f s is dependent on the molar volume for a small guest, we calculate the chemical potential of water for either case prescribed by Equation (18) or (19). The lower chemical potential of water, μ c , under the equilibrium volume at either y l or y s is plotted in Figure 3(a) at  Table 3. The chemical potential of water, m c 0 , at 100 MPa is higher than that at 10 mPa by 2 kJ mol −1 , which arises mainly from the pV term. The water/hydrate boundary of CH 4 as well as Xe is located at y=y f (=0.148) but that of C 2 H 6 along with CF 4 is at y=y l (=0.115) even under extremely different pressures for CS-I clathrate hydrates. It is evident that the boundary is associated intimately with the size of the guest species relative to the cage to accommodate it. A large C 2 H 6 molecule (CF 4 as well) is difficult to fit into the small 12-hedral cage.
A similar calculation is performed for Ne, Ar, Kr, N 2 , and O 2 hydrates, which belong to CS-II. The size of a Kr atom is slightly smaller than CH 4 and it forms CS-II. In those hydrates, the small cages are preferentially occupied rather than the large cages. The chemical potential of water for each hydrate (dotted line in Figure 3 indicating the occupancy of the small cages) shows any hydrate forming CS-II examined here has its water/ hydrate boundary at y=y f (=0.15). The stability of Ne hydrate at 100 MPa is marginal according to our estimation using the LJ parameters in Table 1. Experimental observation shows Ne hydrate is stabilized even around 240 K by the large cages accommodating more than one Ne atom at high pressures, ∼350 MPa (Falenty et al. 2014).

Phase Boundaries for Water/Hydrate and Hydrate/Guest
We examine the equilibria at temperatures mostly lower than the triple point of each guest. Hence, the pure guest is either in the solid or gaseous state except for the high-temperature and high-pressure region. Both the water/hydrate and hydrate/ guest boundaries at 1 kPa are depicted in Figure 4. The two boundaries for each guest are smoothly connected to those at 0 K and are curved to have an intersection (three-phase coexistence). The temperature at the intersection increases with compression. The water/hydrate boundary for the Ar, Xe, and CH 4 hydrate is close to the limiting value of full occupancy. On the other hand, the two boundaries are separated by a large gap for the C 2 H 6 and CF 4 hydrate.

Dissociation Pressures
It is essential to confirm whether our calculation is capable of describing the phase behaviors at cryogenic conditions. To the best of our knowledge, no measurement is available below 80 K except for one report that raises a significant concern as to its reliability (Choukroun et al. 2019;Ghosh et al. 2019). The calculated dissociation pressures for CH 4 hydrate are compared with those from experiments in Figure 5 (Delsemme & Wenger 1970;de Roo et al. 1983;Nakano et al. 1999;Nakamura et al. 2003;Anderson 2004). Although the energy parameter of the LJ potential is so optimized to recover the experimental dissociation pressure only around 273 K (Table 1), overall agreement is clearly seen. The dissociation pressure at the lowest temperature of experimental studies, 82 K, is 0.02 Pa (Delsemme & Wenger 1970). At this condition, CH 4 in equilibrium with hydrate is the gaseous state, practically ideal gas, and further cooling of CH 4 to 40 K along the dissociation curve does not accompany any phase transition of the coexisting CH 4 to the other state. The dissociation pressure according to our calculation down to 80 K seems to exhibit almost perfect agreement with a simple extrapolation of the experimental curve. This corroborates the present method, and the calculated dissociation pressures at lower temperatures are reliable.
Encouraged by the excellent agreement of the temperature dependence of the dissociation pressure for CH 4 hydrate down to 80 K with that from the experiments, we calculate those for various hydrates. Generally, a stronger interaction between a host and a guest molecule results in a lower dissociation pressure as shown in Figure 6. In contrast to the separation of the water/hydrate and hydrate/guest boundaries found for C 2 H 6 and CF 4 at 0 K in Figure 4, no peculiar behavior is observed in the dissociation pressure shown in Figure 6.

Comparison between Spherical and Linear C 2 H 6 Models
So far, all of the guest models are made of a single LJ interaction site. This simplification can be justified by the fact that each type of cage is roughly spherical, and the rotational motion of an anisotropic guest can be described by free rotation inside a cage at high temperatures. A spherical approximation of the CH 4 and CF 4 molecules seems to be legitimate even at low temperatures. However, a C 2 H 6 molecule is rather anisotropic, and the peculiar phase behavior, particularly the large gap between the two boundaries along the horizontal axis in Figure 4, should be examined with a more realistic model.
We adopt for the C 2 H 6 molecule a dumbbell model made of two CH 3 groups, each of which is an interaction site of the LJ potential. The free energy of cage occupation is obtained in the same way as described in Section 3.4. That is, the harmonic frequencies are calculated for the clathrate hydrate composed of water and linear guest molecules according to Equation (14). The energy barrier along one of its rotational coordinates in the 14-hedron cage is very low, around 1 kJ mol −1 (Tanaka 1994), owing to its oblate shape, and it can be approximately a free Note.The chemical potential of the empty state, m c 0 , is calculated at each stage of occupancy, allowing a change in the volume.
rotor. Each C 2 H 6 molecule stays at its optimum orientation at 0 K, and its motion can be described by harmonic oscillations with very low-frequency modes. This can be the highest free energy extreme. The free energy thus calculated provides the upper bound. However, it may not be constrained to that orientation because of the tunneling effect under such a low barrier. For the other extreme, i.e., a rigid free rotor, the lowest energy level is zero, which contributes nothing to the free energy at 0 K. Therefore, the lower bound is obtained possibly    by setting the frequencies of half the number of rotational degrees of freedom to those of the free rotors. The difference is, however, small, less than 0.1 kJ mol −1 . Figure 7 depicts the chemical potentials of water at 0 K in various compositions for the linear and the spherical C 2 H 6 models. The higher chemical potential derived from the harmonic oscillators of all guest motions (the upper bound) is drawn for the linear C 2 H 6 model. The chemical potential at y l =0.115 is lower than the chemical potential of ice Ih, indicating that there is a large gap between the water/hydrate and hydrate/guest boundaries in T-y space for the linear model. The single LJ site approximation for C 2 H 6 is a kind of mean spherical model and works well if an appropriate choice of interaction parameters is made.

Implication of Phase Behavior for Storage and Emission of Volatiles in Satellites
On the surface of Titan, the temperature is around 94 K and the partial pressure of CH 4 is on the order of a kilopascal (Niemann 2005;Tobie et al. 2012). This pressure is much higher than the dissociation pressure of CH 4 hydrate according to our estimation shown in Figure 8(a). The dissociation pressure for O 2 hydrate is also calculated and compared with the conditions of Europa and Ganymede in Figure 8(b) (Hall et al. 1995;Orton et al. 1996;Delitsky & Lane 1998;Hall et al. 1998;Saur et al. 1998). The pressure on the surface of either Europa or Ganymede, whose atmosphere is made of almost 100% of O 2 , is substantially lower than the dissociation pressure of O 2 hydrate at the corresponding temperature.
A further complexity occurs when CH 4 is mixed with N 2 in equilibrium with clathrate hydrate. In fact, the gas mixture of N 2 (95%) and CH 4 (5%) is the main component of the atmosphere on Titan's surface. The dissociation pressure of the mixed hydrate is calculated from where x kj 0 is the occupancy of the guest species k in a cage of type j in equilibrium with the guest gas mixture. It is calculated from the chemical potential of the guest species as where f kj is the free energy of cage occupation by the guest species k in a cage of type j and m gk 0 is the chemical potential of the guest species k in a gaseous mixture. Thus, we can examine the effect of N 2 on the stability of various compositions of CH 4 in the atmosphere of Titan. Plotted in Figure 8(a) are the dissociation pressures for mixed hydrates of CS-II structure with which the gaseous mixtures coexist. The mole fractions of CH 4 are chosen to be 3%, 5%, and 10%. The atmospheric pressure on the Titan is significantly higher than the dissociation pressure.
It has been reported that water (ice) is partially exposed on the surface of Pluto, which is covered mostly with solid N 2 (Schmitt 2017). Provided that water exposed or underneath the surface is in equilibrium with its atmosphere, its form is necessarily clathrate hydrate as shown in Figure 8(a) although the formation of clathrate hydrate may be suppressed by slow kinetics considering the cryogenic condition on the surface of Pluto (Gladstone 2016). A large amount of clathrate hydrate acting as a thermal insulator should serve as a heat reservoir to account for its possession of a subsurface ocean (Kamata et al. 2019). A similar analysis was carried out for Titan to account for the thermal insulation effect (Kalousová & Sotin 2020). It is worth exploring in the future whether the anticipated clathrate hydrate on the subsurface of the ocean is thermodynamically connected to that at the surface without crossing the three-phase  The water/hydrate and hydrate/guest boundaries are drawn in Figure 9 for CH 4 , N 2 , and O 2 hydrates. Only the hydrate/ guest boundary is shown for the CH 4 +N 2 mixed hydrate as the water/hydrate boundary is specified by an arbitrarily chosen composition of CH 4 and N 2 and is out of our current interest. The hydrate/guest boundary is calculated from Equation (22). Because the atmosphere of Titan is nearly 100% N 2 + CH 4 gas mixture at a total pressure of 150 kPa, it is in equilibrium with the hydrate whose composition is given by the intersection with the horizontal line corresponding to the surface temperature (94 K) depicted by the dotted-dashed line in Figure 9. The fraction of CH 4 to the total guest content in the mixed hydrate on that boundary is more than 90% as calculated from Equation (22), which is quite large compared with that in the atmosphere, 5%. On the other hand, either horizontal line corresponding to the surface temperature of Europa (50 K) or Ganymede (70 K) does not intersect with the hydrate/guest boundary for O 2 hydrate at a pressure of 10 −7 Pa. Therefore, O 2 hydrate cannot coexist with the atmospheric air on the surface of either satellite.
These equilibria imply that water on the surface of Titan is all in the form of clathrate hydrate owing to the two-phase coexistence on the hydrate/guest boundary. Contrary to Titan, a stable compound is ice Ih, and a clathrate hydrate cannot exist on the surface of either Europa or Ganymede due to the extremely low O 2 pressure. A similar view on the surface of Titan has been reported with the extrapolation of the threephase coexistence curve to the lower temperature region with empirical parameters (Osegovic & Max 2005). Here, we scrutinize by exploiting a rigorous statistical mechanics theory that the form of water on the surface of each satellite is estimated in terms of two-phase equilibrium, which occurs at thermodynamic conditions much distant from the three-phase equilibrium.
We assume that equilibrium is attained everywhere in the crust and the aqueous solution of Titan. Under this circumstance, we infer a possible phase behavior of CH 4 clathrate hydrate in the crust of Titan and at the subsurface from what we have estimated in the present calculations. Figure 10(a) plots the calculated three-phase coexistence curve, which agrees  (Gladstone 2016), and crosses for Europa and Ganymede (Orton et al. 1996;Hall et al. 1995;Delitsky & Lane 1998;Hall et al. 1998;Saur et al. 1998). The most favorable condition is chosen to form a hydrate in (b). Figure 9. The water/hydrate (solid lines) and the hydrate/guest boundaries (dotted lines) at 1.5×10 5 Pa for CH 4 , N 2 , and CH 4 +N 2 hydrates and at 10 −7 Pa for O 2 hydrate against the mole fraction of the guest species. For the mixed hydrate, only the hydrate/guest boundary is drawn against the sum of the mole fraction of the two guest species. Each dashed-dotted line corresponds to the temperature of the corresponding satellite surface. Two vertical lines for O 2 hydrate overlap.
excellently with the measured ones as shown in Figure 5. This curve tells us that CH 4 clathrate hydrate can form only in the upper region separated by the curve. However, it cannot specify whether the clathrate hydrate alone is a stable phase or it coexists with the water-rich or the guest-rich phase.
The temperature-pressure profile in the crust is assumed tentatively by the line connecting the two positions, the surface (0.15 MPa, 94 K) and the subsurface of the ocean (125 MPa, 262 K) after the bulk composition data (Tobie et al. 2012). This line is also plotted in Figure 10(a). It may deviate more or less from the curve but is unlikely to have any crossing point with the three-phase coexistence curve. Liquid water containing some salts and CH 4 (other hydrophobic molecules as well) below the subsurface can be in equilibrium with clathrate hydrate as the thermodynamic state is located in the upper region and the solubility of CH 4 is low. That is, the rightmost point on the temperature-pressure curve in the crust (cyan filled circle in Figure 10(a)) corresponds to the three-phase coexistence of clathrate hydrate with the aqueous solution and ice Ih. This three-phase equilibrium should be distinguished from the threephase equilibrium of water (either liquid or crystalline), guest-rich phase, and clathrate hydrate. The leftmost point on the curve corresponding to the thermodynamic state at the surface is distant from the three-phase coexistence curve (blue curve), indicating that the atmosphere is necessarily in equilibrium with clathrate hydrate. On the way from the subsurface to the surface, there must be a state (layer) where only clathrate hydrate is stable. Such a state (layer) must intervene between two different coexistence equilibria, clathrate hydrate with the atmosphere at the surface and with ice Ih in the deeper region of the crust, no matter how thin the layer is. As a consequence, a large amount of CH 4 can be stored as clathrate hydrate coexisting with ice Ih even in the deep crust.
Let us consider a simple case where only CH 4 is dissolved in pure liquid water. Salts and ammonia decrease the formation temperature of CH 4 clathrate hydrate. However, this effect will not change the following argument. The temperature-composition diagram is drawn in Figure 10(b) at low temperatures. The CH 4 fraction in the aqueous solution is believed to be small, lower than the solubility limit considering the total content of water and CH 4 in Titan (Tobie et al. 2012). However, clathrate hydrate can form in such an aqueous solution of CH 4 . This is because the solubility of CH 4 in the solution is low in the presence of clathrate hydrate, which cannot exceed 10 −2 around 273 K from our previous calculation, and decreases with cooling contrary to that in the aqueous solution in contact with the CH 4 -rich phase (Tanaka et al. 2018b). Thus, the mean composition of CH 4 is much smaller than the mole fraction of CH 4 in the coexisting clathrate hydrate at the top of the subsurface, which implies that the water/hydrate two-phase equilibrium is achieved at the top of the subsurface. This coexistence condition is maintained even when liquid water is transformed into ice. The content (cage occupancy) of CH 4 in the clathrate hydrate increases with decreasing temperature and reaches the limiting value on approaching the surface along the solid line in Figure 10(b) as long as equilibrium is attained. A small perturbation may make it possible to penetrate the narrow one-phase region enclosed by the solid and dotted lines around 90 K and to arrive at the hydrate/guest boundary.
We also examine the temperature dependence of the density of clathrate hydrate. It is shown in Figure 11 (in g cm −3 ) for the hydrate/water and the hydrate/guest boundaries at 1 kPa for various clathrate hydrates. Only a small difference is observed for Ar, Kr, and CH 4 hydrates between the two boundaries while a large discrepancy is found for C 2 H 6 and CF 4 hydrates. The density is affected by the temperature in two ways: the change in occupancy and the thermal expansion under a fixed composition. It is evident from Figure 4 that the former is much more significant and the difference between the two boundaries arises mainly from the change in occupancy with temperature.
The density relative to that of liquid water is an important property with regard to the distribution of clathrate hydrates in internal (salty) oceans; a higher density clathrate hydrate sinks to the base of the ocean and a lower density one is stored at the top. It is expected that Kr and Xe hydrates above the melting temperature of ice Ih are much denser than liquid water while CH 4 hydrate stays at the subsurface. A detailed investigation on the distribution of hydrates in the internal oceans is out of the scope of the present paper and will be treated in the future.

Concluding Remarks
We investigate the clathrate hydrate stability at cryogenic temperatures down to the 0 K limit where the guest species is a crystalline solid. To this end, we develop a statistical mechanics theory to treat the cryogenic conditions with reliable intermolecular interactions but without introducing the empirical parameters of the thermodynamic properties required for evaluation of the stability. The pressure range covers 10 −12 to 10 8 Pa: from the value at interstellar space to that at the surface of the hydrosphere in icy satellites of the outer planets. A series of noble gases, small hydrocarbon and its derivative, and two common ubiquitous gases N 2 and O 2 are chosen as the guest species. All of the guest species are treated as spherical molecules each having a single interaction site except for C 2 H 6 , which is also treated as a two-site linear molecule.
We first explore the phase behaviors of those hydrates at 0 K. A large guest such as C 2 H 6 and CF 4 has the water/hydrate and hydrate/guest boundaries whose compositions are substantially different. That is, the small cages are left totally unoccupied while the large cages are all filled. This gives rise to an unexpectedly large tolerance of the guest composition in stable hydrate. On the other hand, the two boundaries for the smaller guests Ar, Kr, N 2 , O 2 forming CS-II and CH 4 , and Xe forming CS-I converge to the composition corresponding to the full occupancies of two sorts of cages.
Next, we calculate the two-phase boundaries at temperatures below 200 K. It is shown that the boundaries are extended smoothly from 0 K. It is also found that the dissociation pressure of any clathrate hydrate investigated here decreases significantly on cooling to cryogenic temperatures. We examine a smooth connection of the dissociation pressure from 0 K to that calculated based on the classical partition function above 200 K and compatibility with experiments at low temperatures for CH 4 hydrate. The smooth connection and the excellent agreement suggest that the present method applied to cryogenic conditions is reliable.
For several large guests such as C 2 H 6 and CF 4 , which may be encapsulated in either a large or small cage, the stable state of clathrate hydrate is allowed to exist between two boundaries of the two-phase coexistence, which are separated by a large gap along the y-axis on the y-T plane. At very low temperatures, the condition of its formation, i.e., a water-rich or guest-rich environment, may be recorded in its composition for a fairly long time and then the composition can be adjusted slowly toward the equilibrium composition constrained by the current external environment.
The implication of the phase behaviors of clathrate hydrates for the stable form of water in Titan is that water on the surface exists only as clathrate hydrate in equilibrium with its atmosphere where the partial pressure of CH 4 is sufficiently high to form clathrate hydrate. This situation continues down to a shallow region of the crust, plunging into the one-phase region crossing the hydrate/ guest boundary. However, water ice in the deeper part of the crust still coexists with clathrate hydrate moving beyond the water/ hydrate boundary, so that the two-phase coexistence between ice Ih and clathrate hydrate in the lower part of the crust is continuously connected with the condition at the subsurface. Whereas the surface of Titan is covered by clathrate hydrate, the thin oxygen air coexists exclusively with pure ice at the surfaces of Europa and Ganymede.
Apart from a controversial topic on the rapid formation of clathrate hydrates under interstellar conditions (Choukroun et al. 2019;Ghosh et al. 2019), the equilibrium properties can be obtained from the present theoretical calculations at cryogenic conditions. However, the migration of CH 4 from the subsurface ocean to the surface involves various processes. The temperature dependence of the solubility of guest species in liquid water is opposite that of water in equilibrium with a guest-rich phase, and this may influence the formation of clathrate hydrate in the subsurface ocean (Zatsepina & Buffett 1997;Servio & Englezos 2002;Tanaka et al. 2018b). A further study is required to obtain a conclusive perspective on what extent the stability of CH 4 + N 2 hydrate is involved in the emission of CH 4 inside the ice crust with variable temperature, pressure, and composition, other than by examining the timescale of its formation and dissociation around 90 K.
We are grateful to Prof. G.L. Hashimoto for valuable discussions on CH 4 hydrate on Titan. The present work was supported by JSPS KAKENHI grant No. 17K19106 and the Research Center for Computational Science in providing computational resource.

Appendix
Nomenclature for physical quantities are listed below. The free energies and chemical potentials are all absolute values with reference to the intact molecules separated infinitely.
T: Temperature p: Pressure V : Volume N w : Number of water molecules Figure 11. Densities of Ar, Kr, Xe, CH 4 , C 2 H 6 , and CF 4 clathrate hydrates against temperature at a pressure of 1 kPa on the hydrate/water boundary (solid line) and the hydrate/guest boundary (dotted line).