LOW-TEMPERATURE EQUATION OF STATE OF SOLID METHANE

The theoretical equation of state for solid methane, developed within the framework of perturbation theory, with the crystal consisting of spherical molecules as zero-order approximation, and octupole – octupole interaction of methane molecules as a perturbation, is proposed. Thermodynamic functions are computed on the sublimation line up to the triple point. The contribution of the octupole – octupole interaction to the thermodynamic properties of solid methane is estimated.


I. INTRODUCTION
Methane (CH 4 ) is one of the most common hydrocarbon compounds, which is the main component of natural gas [1].Knowledge of its properties in the condensed state is necessary for a wide range of applications: from the calculation of processes in pipes and refrigeration systems to the description of the structure of the giant planets, consisting largely of condensed (including crystalline) methane [2,3].Solid methane is also used as a neutron moderator [4,5].
In recent years due to the development of the alternative energy sources and the search for the molecular systems on the basis of new carbon materials suitable for storage and easy retrieval of molecular hydrogen as the fuel of the future, there is a growing interest to crystals, composed of XY 4 type molecules (CH 4 , CCl 4 , CF 4 , SiH 4 ) [6].The presence of methane in clathrate systems is also a positive factor, in particular, due to the high content of hydrogen in it.
Properties of gaseous and liquid methane were intensively studied during last decades.Equations of state for the liquid and gas phases are well presented in the literature [7,8], and the accuracy of the description of the thermodynamic properties is satisfactory.
However, in the design of refrigeration and cryogenic facilities it is necessary to estimate the solid phase appearance and stability conditions to avoid unwanted precipitation of methane inside them during their operation.Therefore, when designing cryogenic installations and equipment for separation of atmospheric gases, one requires information on the thermodynamic properties of methane not only in gaseous but also in condensed phase, and above allon the line of sublimation.
The aim of this paper is to develop the equation of state for solid methane, allowing to calculate its thermodynamic properties of a wide range of temperatures and pressures, starting from the sublimation up to the melting line, as well as to predict the properties of solid methane in a wider range of pressures and temperatures.

II. THEORETICAL EQUATION OF STATE FOR A CRYSTAL
Solid methane exists in three different allotropic modifications [9,10], having different crystal lattice structures.One of the phases of crystalline methanea high-temperature -phase is characterized by an almost free rotation of molecules in sites of a face-centered cubic (FCC) lattice.
This phase is stable on the sublimation line from 20.5 K to the triple point temperature (90.7 K), and at higher temperatures on the melting line up to high 81 pressures, where the transition to an orientation-ordered phase occurs.
According to [9,10,11] the pressure of  transition is about 10 kbar at 70 K, and about 15 kbar at 120 K, and continues to increase with further increase of temperature.Orientation ordering at low temperatures and high pressures is associated with the weak anisotropy of the intermolecular interaction of CH 4 molecules.Far from the transition line the noncentral nature of the interaction has little effect on the behavior of crystalline methane.In this case the main contribution of non-spherical forces is octupoleoctupole interaction of CH 4 molecules.
Due to the relative weakness of this interaction, the properties of the high-temperature phase of solid methane are similar to those of solidified inert gases, which allows to apply an equation of state, obtained for the Lennard -Jones FCC crystal [12], to calculate thermodynamic properties of methane.However, in this case the correction for octupoleoctupole interaction of methane molecules should be taken into account.The relative smallness of this correction allows to apply the classical thermodynamic perturbation theory [13], with the FCC-crystal consisting of spherical molecules appears as the zero-order approximation, and octupoleoctupole interaction of methane molecules as a perturbation.Such an equation is proposed below.
The Helmholtz free energy of a crystal can be written as the follows: is the Helmholtz free energy of the FCC-crystal, build from spherical molecules.It can be calculated using the theoretical equation of state, developed in [13].
The general approach to the derivation of the equation of state for strongly anharmonic crystals [13] is based on a generalization of the Mayer's group expansion.
The Helmholtz free energy is written as a sum of a single-particle contribution and some corrections on correlations in displacements of pairs W 2 , triples W 3 ... of molecules from their lattice sites: Here is the Helmholtz free energy within Einstein (single-particle) approximation, which is expressed in terms of so-called "free volume" f v [14], and W 2 , W 3 ,… accounts, respectively, for double, triple, etc. correlations between the displacements of molecules from their lattice sites.
Free volume and correlation corrections W 2 , W 3 in the equation ( 1) are expressed, respectively, by the three-, six-and nine-fold integrals, and special numerical methods are necessary to evaluate them.Such evaluation of free volume, as well as of correlation corrections was carried out numerically in [13] and demonstrated adequate account of rather strong anharmonicity of molecular vibrations in Lennard-Jones solid.
The accuracy of the theoretical equation of state was verified by comparison of the free energy calculations with available Monte Carlo computer simulation data [14] and a good agreement was demonstrated [13].Thus, the adequacy of the description of the strongly anharmonic crystals by theoretical equation of state (1) has been shown.
An analytical approximation of the canonical equation of state for strongly anharmonic crystals proposed recently in [12] allows performing practical calculations easily.The equation has been written in a simple analytical form, similar to what was proposed by Van der Hoef [15]: An attempt to extend the limits of applicability of the Van der Hoef equation [15] was made in [12].The upper limit of the index n was increased and all the coefficients were found by minimizing the deviations between the predictions of the theoretical equation of state and the equation (2).
Using well known thermodynamic relations one can easily calculate on the basis of the canonical equation of state (2): pressure (thermal equation of state), internal energy (caloric equation of state) and entropy, heat capacities C P and C V , temperature expansion α T and isothermal compressibility β T coefficients.Equation (2) provides a good approximation of the computer simulation data at reduced temperatures T* = 0.1…2.0 and reduced densities *= 0.6… 1.39 [12].
Numerical values of n b and nm a coefficients obtained in [12] are listed in Table 1.

III. OCTUPOLE-OCTUPOLE CONTRIBUTION TO THE SOLID METHANE EQUATION OF STATE
Within the framework of perturbation theory [16], a crystal consisting of spherical molecules corresponds to the zero-order approximation, and octupoleoctupole interaction of methane molecules is a perturbation.
The canonical equation of state includes

), corresponding to the
Lennard-Jones model of crystal with the FCC lattice as a major contribution to the Helmholtz free energy and a perturbation octupoleoctupole correction: Here oct U  stands for the octupoleoctupole interaction of methane molecules, averaged over the orientation and distribution functions of the Lennard-Jones FCC-solid.
Averaging over orientations can be provided in the same manner as for the system of free rotators.In this case, the multipole symmetry of the non-center interaction leads to the fact that the first-order correction oct U  vanishes and the effect of the non-spherical interaction on thermodynamic properties of solid methane provides the mean-squared value of the octupole -octupole interaction: Further calculations are similar to the calculations carried out for the liquid phase [18].Averaging of the octupole -octupole potential according [18] over orientations and crystalline sites gives the following:  (5) where S 14 is so-called lattice sum [18], Using this expression, one can calculate the necessary corrections to all thermodynamic functions.The appropriate kinetic contribution due to the free rotation of the CH 4 molecules should also be included in this contribution.

IV. THERMODYNAMIC PROPERTIES OF SOLID METHANE
To calculate the thermodynamic properties of solid methane using the equations ( 2) -( 5), we need to set values of the Lennard -Jones potential parameters  and  , as well as the quadrupole moment  .
The value of the octupole moment  = 4.510 -34 units of charge/cm 3 was adopted according to [7].In the literature there is a number of parameters of the Lennard-Jones potential , approximating the interaction of CH 4 -CH 4 .But these values effectively take into account the octupole -octupole interaction.Because the octupole-octupole contribution decreases with increasing temperature, the most suitable method of defining of  and  values is use of high-temperature data.We have adopted the high-temperature values of the CH 4 -CH 4 effective interaction potential parameters provided in the monograph [19].At T > 500 K these values became constant: .Specific volumes, internal energy, enthalpy, entropy, thermal expansion and isothermal compressibility coefficients, as well as heat capacities calculated in the range of temperatures (40-90 K) on the sublimation line (at P = 0) are shown in Table 2.The values of the caloric functions are measured from the triple point.

V. DISCUSSION OF THE RESULTS
The proposed equation of state has no specific constants fitted using experimental data on solid methane.Therefore the comparison to the experiment became especially interesting.
Calculated and experimental data [21] are compared on sublimation line in Table 2 and Figs.1-4.Equation of state reproduces specific volume on the sublimation line with an average accuracy of less than 0.65%, and the isothermal compressibility factor within 4.3%.Somewhat more significant deviations from the experiment have been found for the heat capacity C V (8.4%) and temperature expansion coefficient (29%) at lower temperatures.This is due to quantum effects significantly manifested in these properties and not taken into account in this model.Better agreement of calculated heat capacity with the experiment data on the sublimation line was observed near triple point.
The comparison of the calculated volume expansion coefficient with dilatometric and X-ray measurements is presented in Fig. 2. The deviation from experimental data near the triple point (90 K) is traditionally explained by the presence of vacancies in the premelting area.values [21] of isothermal compressibility on the sublimation line of methane.
Our data are in better agreement with X-ray measurements, because in the equation of state the contribution of vacancies was not taken into account.
The calculated value of the isothermal compressibility coefficient (Fig. 2) agrees with the experiment within the measurement error, which can be considered as more than satisfactory for a theoretical equation containing no fitting parameters.
We carried out also assessment of the octupoleoctupole interaction contributions to thermodynamic properties of solid methane.In the considered temperature range this average contribution to the calculated properties of the crystal is: -0.5% (density), 6% (thermal expansion), 2% (compressibility), 7 % (heat capacity).Thus, we can conclude that the theoretical equation of state, presented by equations ( 1) -( 5), having no adjustable parameters, allows quite reliably calculate thermodynamic functions of high-temperature phase of solid methane at least on sublimation line.
free energy, temperature and density referred to the Lennard-Jones potential parameters  and ,  = 3.76 Å were used in calculation of the liquid methane properties[20].The reduced octupole moment of CH 4

Figure 1 -
Figure 1 -Comparison of the calculated (solid line) and experimental values [21] of the specific volume on sublimation line.

Figure 2 -
Figure 2 -Calculated (solid line) and experimentalvalues[21]  of isothermal compressibility on the sublimation line of methane.

Figure 3 - 1 )
Figure 3 -Volume expansion coefficient of solid methane on sublimation line.Comparison of the calculated (solid line) and experimental values[21].

Table 1 .
Coefficients of the Lennard-Jones crystal equation of state(2).

Table 2 .
Specific volume V, enthalpy H, entropy S, temperature expansion  P and isothermal compressibility  T coefficients, and heat capacities C V and C P ) of solid methane calculated on sublimation line.