domain: spi phys math info

text:
R Leroy
email: remy.leroy@onera.fr

P Trouvé-Peloux

F Champagnat

B Le Saux

M Carvalho
email: marcela.carvalho@upciti.com

Pix2Point: Learning Outdoor 3D Using Sparse Point Clouds and Optimal Transport

Good quality reconstruction and comprehension of a scene rely on 3D estimation methods. The 3D information was usually obtained from images by stereophotogrammetry, but deep learning has recently provided us with excellent results for monocular depth estimation. Building up a sufficiently large and rich training dataset to achieve these results requires onerous processing. In this paper, we address the problem of learning outdoor 3D point cloud from monocular data using a sparse ground-truth dataset. We propose Pix2Point, a deep learning-based approach for monocular 3D point cloud prediction, able to deal with complete and challenging outdoor scenes. Our method relies on a 2D-3D hybrid neural network architecture, and a supervised end-to-end minimisation of an optimal transport divergence between point clouds. We show that, when trained on sparse point clouds, our simple promising approach achieves a better coverage of 3D outdoor scenes than efficient monocular depth methods.

Introduction

Geometry estimation is a prior requirement for autonomous agents to comprehend their environment, progress and interact within it. This fundamental task leads to various estimation methods for 3D reconstruction. For long, the parallax between images has been used [START_REF] Faugeras | Three-Dimensional Computer Vision: A Geometric Viewpoint[END_REF][START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], either between two cameras in stereo mode or between multiple acquisitions after displacement. Recently, deep learning techniques have revolutionised 3D estimation from images, allowing even to obtain excellent results with a single view [START_REF] Jahani Amiri | Semi-supervised monocular depth estimation with leftright consistency using deep neural network[END_REF][START_REF] Carvalho | On regression losses for deep depth estimation[END_REF][START_REF] Eigen | Depth map prediction from a single image using a multi-scale deep network[END_REF][START_REF] Fu | Deep ordinal regression network for monocular depth estimation[END_REF][START_REF] Han Lee | From Big to Small: Multi-scale local planar guidance for monocular depth estimation[END_REF]. These impressive results rely upon large and highly accurate databases like KITTI [START_REF] Geiger | Vision meets robotics: The KITTI dataset[END_REF] that involve the simultaneous collection of stereo pairs and LiDAR data, post-processing and temporal integration in order to provide accurate, reliable and dense ground-truth for learning purposes. The overall process requires large scale cooperation and is therefore lengthy. In this paper, we address the question of performing 3D estimation from monocular data with much lighter requirements in terms of training data input. Besides, if monocular point sets prediction has been proposed in the literature, it is only for the reconstruction of a single 3D object [START_REF] Fan | A point set generation network for 3d object reconstruction from a single image[END_REF][START_REF] Mandikal | Dense 3D point cloud reconstruction using a deep pyramid network[END_REF][START_REF] Pumarola | C-flow: Conditional generative flow models for images and 3d point clouds[END_REF].

In contrast, we propose here (i) one of the first approaches to reconstruct a 3D point cloud for an entire outdoor scene given only a single image using a 2D-3D hybrid neural network architecture, illustrated in Fig. 1; (ii) an end-to-end learning scheme of the hybrid model using a sparse point clouds dataset; (iii) a cost function from the optimal transport (OT) that makes possible to obtain good coverage of the scenes on the KITTI dataset [START_REF] Geiger | Vision meets robotics: The KITTI dataset[END_REF] (iv) better performances than state-of-the-art monocular depth map prediction methods when trained with sparse data. This illustrates the advantage of the direct prediction of point clouds, in the case of a low-density ground-truth 3D dataset.

Related Work

The monocular 3D estimation task has been addressed in terms of depth maps prediction with the following works [START_REF] Shariq | Adabins: Depth estimation using adaptive bins[END_REF][START_REF] Carvalho | Deep Depth from Defocus: How can defocus blur improve 3D estimation using dense neural networks?[END_REF][START_REF] Eigen | Depth map prediction from a single image using a multi-scale deep network[END_REF][START_REF] Fu | Deep ordinal regression network for monocular depth estimation[END_REF][START_REF] Han Lee | From Big to Small: Multi-scale local planar guidance for monocular depth estimation[END_REF][START_REF] Zhu | Monocular depth prediction through continuous 3D loss[END_REF][START_REF] Saxena | Learning depth from single monocular images[END_REF] but as already mentioned, they are trained with large and rich datasets. Recently, 3D estimation has also been addressed in terms of 3D point set with PSGN by Fan et al . [START_REF] Fan | A point set generation network for 3d object reconstruction from a single image[END_REF], a method that aims to predict the envelope of an object as an unordered set of points using a single view of it and its location in the image. Mandikal and Babu [START_REF] Mandikal | Dense 3D point cloud reconstruction using a deep pyramid network[END_REF] address the limitations of PSGN regarding the poor number of predicted points with DensePCR, a pyramidal structure allowing to enlarge the number of points. Xia et al . [START_REF] Xia | Realpoint3d: Point cloud generation from a single image with complex background[END_REF], also tackles the generation of a monocular point cloud for objects using prior knowledge over their shapes, making it robust to occlusions and varying poses. A generative flow-based model allowing single view object point cloud prediction has recently been proposed with C-flow [START_REF] Pumarola | C-flow: Conditional generative flow models for images and 3d point clouds[END_REF]. It leverages a back-and-forth prediction loop from image to point cloud, then to image for consistency.

The aforementioned point cloud works only consider the reconstruction of a single 3D object model, i.e., on data that are obtained through demanding procedures, either scanned objects using RGB-D sensors or laser scanners, or handcrafted models. These procedures do not apply to real-life outdoor scenes with various settings and where lies multiple objects. [7] tackles the problem of monocular volumetric reconstruction with occlusion completion for indoor scenes. In addition to the input RGB image, this approach requires a corresponding normal image that is hardly obtainable for outdoor scenes. To our knowledge, we are the first to propose a deep learning method for the reconstruction of complex outdoor scenes in the form of 3D point cloud, solely conditioned by a single RGB image, and learned on sparse point clouds. Our technical contributions with respect to previous works are: a performance study of the 2D-3D hybrid neural network architecture for point cloud reconstruction of scenes "in-the-wild", end-to-end training with OT loss optimisation on sparse point clouds up to ten thousands elements. Finally, we compare our point set prediction approach to state of the art depth map prediction approaches, when they are all trained on the same sparse dataset.

Method

Given a single colour image, our method, namely Pix2Point, predicts a set of 3D point coordinates with an arbitrary number of elements fixed before training. Fig. 1 shows an overview of the proposed method. Like DensePCR, our architecture consists of an encodingdecoding module to predict a first coarse point cloud that will be enriched using a densification module. Our model, unlike DensePCR, is trained to minimise an optimal transport divergence over the point coordinates in an end-to-end fashion.

Encoding The encoding block is a series of convolution, pooling and normalisation layers to extract a feature description of the full RGB input image, which is then processed by a fully connected layer to obtain a first coarse set of 3D point coordinates. We explore and compare several encoding approaches following either VGG, DenseNet and ResNet architectures in section 5. We refer to them as backbones.

Decoding We refer to decoding as the densification of the first coarse 3D point cloud. We duplicate every point n times, and to describe each point we concate-nate: the 3D coordinates, both global point cloud feature vector and the corresponding local feature vector, obtained using dedicated PointNet-like shared multilayered perceptron (MLP) [START_REF] Ruizhongtai | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF][START_REF] Ruizhongtai | PointNet++: Deep hierarchical feature learning on point sets in a metric space[END_REF], and lastly a grid alignment feature vector in order to identify every clone of the same point and to suggest geometric information between every clone. This point description is processed by another shared MLP resulting in 3D coordinates for 1 point.

Optimisation The performance of our approach mainly depends on the training loss function. Unlike the depth map prediction methods which exploit the gridded structure of the image for the evaluation of errors, our method uses distances between unordered point sets. In particular, it requires an additional, computationally expensive step to match points from the predicted and the target point clouds. We expose two usual distances for this task.

Chamfer distance The chamfer distance is the average of squared Euclidean distances to the nearest neighbour from one set to the other.

Optimal transport or OT distance To compute this distance, also known as Earth Mover's distance, one has to find a one-to-one mapping from one set to the other that minimises the sum over each point of the squared distance between them and their corresponding image. This minimised sum is called the OT distance and it informs about the eventual discrepancy between point sets distributions.

The exact computation of an OT distance is time and memory expensive especially for several thousand elements, hence, we consider in our work an approximation of the OT distance obtained by adding a regularising term and solving the Sinkhorn-Knopp algorithm [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Feydy | Interpolating between optimal transport and MMD using sinkhorn divergences[END_REF].

Experimental Settings

This section presents the dataset we are considering for real scene point cloud estimation from a single image and defines several criteria for reconstruction performance evaluation.

Dataset

To assess our method, we operate on RGB image sequences of real urban scenes and corresponding LiDAR point cloud acquisitions from the KITTI depth estimation benchmark dataset [START_REF] Geiger | Vision meets robotics: The KITTI dataset[END_REF]. Every scene point cloud is an accumulation of filtered LiDAR acquisitions over few successive time instants. We use the split defined by Eigen et al . [START_REF] Eigen | Depth map prediction from a single image using a multi-scale deep network[END_REF], that is 22 600 training scenes and 697 testing scenes.

Evaluation Metrics

There is no reference value or method for the KITTI scene reconstruction task regarding 3D point clouds. Therefore, we propose to measure performances using the completeness and accuracy criteria from [START_REF] Tylecek | The second workshop on 3D Reconstruction Meets Semantics: Challenge results discussion[END_REF], and also propose the relative accuracy. All 3 measures are defined as follows:

Completeness is the coverage in per cent of the target point cloud by the predicted points. A target point is covered if a predicted point lies in its surrounding (i.e. fixed radius ball). We evaluate completeness values for radius of 50 cm, 25 cm and 10 cm.

Accuracy is the distance d, in meter, from the r-th percentile of the distances to the nearest neighbour, from the predicted point cloud to the ground-truth point cloud. It measures the greatest distance to the nearest neighbour among the predicted points closest to the ground truth. We choose r < 90% to include most of the points and discard eventual outliers.

Relative accuracy is similar to the accuracy, where every distance to their nearest neighbour is divided by the norm of the corresponding target point. It provides a higher penalty to short-range predictions.

Implementation Details

Experiments were conducted using the Pytorch Framework [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF]. We kept the original image resolution and cropped every picture to 1224 × 370 pixel definition. Due to heavy computational cost for loss back-propagation, parameters were updated after every sample forward, making batch normalisation ineffective. Instance normalisation was applied instead [START_REF] Ulyanov | Instance normalization: The missing ingredient for fast stylization[END_REF]. The number of predicted points was determined to fully load the 8GB GPU during training. Therefore, 2500 elements point clouds are first predicted by the fully connected encoding module, that are then up-scaled by a DensePCR-like module making a point cloud with 10k elements.

Training vs. Testing point clouds: Using an OT loss enables, in principle, the comparison of any ground-truth point cloud to the predicted one, however, in practice, computation and optimisation of an OT loss are computationally much more efficient with point sets of equal cardinality. Therefore, ground-truth point cloud databases are randomly subsampled to 10k points, that is as many points as Pix2Point predicts. When testing we measure performances to the whole ground-truth point cloud.

Experimental Results

In this section we first present the performances of Pix2Point using various encoding backbones and losses, then we compare our method to depth map prediction approaches through evaluation metrics defined in 4.2.

Network parameter study

We trained several models with varying encoding backbones and loss functions. We considered the following configurations: Pix2Point architecture with VGG backbone and training on the minimisation endto-end of either the chamfer or OT distance, and Pix2Point with ResNet backbone and minimising the OT distance. The performances of these models are given in Table 1 respectively as P2P-VGG-C, P2P-VGG-OT and P2P-ResNet-OT. From these figures, we can notice that the minimisation of chamfer distance thrives toward predictions with low local error, and minimising the OT distance grants predictions with higher completeness, hence, better coverage of the scenes. In order to find if these distances could help each other, combinations of both distances have been tested. However, they lead to convergence issues during training and overall worst performances due to opposite objectives of the distances. Changing the backbone from VGG to ResNet has also a slight impact on the completeness and accuracy. The small gain in relative accuracy indicates that far predictions are more accurate. We also provide a comparison to a similar imageto-point-cloud approach, DensePCR [START_REF] Mandikal | Dense 3D point cloud reconstruction using a deep pyramid network[END_REF], initially proposed for 3D graphics models.

Comparison to depth prediction approaches

The current dominant approach to 3D estimation from a single image consists of predicting corresponding depth maps by leveraging the power of image-toimage translation networks. On KITTI, these methods are trained on pseudo-dense depth maps built by accumulating several consecutive LiDAR acquisitions. For comparison in a similar setting, we train two state-of-the-art models for monocular depth estimation from the KITTI challenge1 , namely BTS [START_REF] Han Lee | From Big to Small: Multi-scale local planar guidance for monocular depth estimation[END_REF] and Ad-aBins [START_REF] Shariq | Adabins: Depth estimation using adaptive bins[END_REF], on the same 10k-point-cloud as Pix2Point. At inference time, dense depth maps are projected back to 3D using known camera parameters. Performance comparison with various flavours of Pix2Point is reported in Table 1.

These results show that Pix2Point, with only 15M parameters, trained with Chamfer distance performs better than BTS and Adabins, respectively 45M and 78M parameters. Moreover, when trained with the OT distance, Pix2Point accuracy decreases but it covers three times more points than depth map approaches for the closest neighbourhoods. These observations can be made through Figure 2 where we show point cloud predictions and the coverage error map for each method (for comparison all predictions are visualized with 10k points). This error map displays for each ground-truth point the distance to its nearest predicted point. We display the error from 0 to 50cm using the jet colourmap. While AdaBins and BTS preserve fine features, all Pix2Point variants achieve better coverage of the scene and a lower error, especially for far-away elements and small size features, that are not retrieved by AdaBins and BTS (see for instance the right part of the bottom scene).

Conclusion

We proposed in this work an innovative approach to tackle point clouds reconstruction for complex outdoor scenes from a unique RGB image, using a light 2D-3D hybrid neural network. It recovers properly distributed point clouds by taking advantage of an optimal transport loss. We also provided the first benchmark for this novel task on the KITTI dataset and introduced performance metrics to assess the quality of point cloud reconstruction. We show that our method outperforms state of the art depth map prediction methods, when trained with sparse data, illustrating the interest of direct image to 3D point-cloud translation.

Figure 1 .

 1 Figure 1. Pix2Point: 3D point set prediction for real outdoor scenes from a single image. A first 2D CNN module encodes the RGB image for the following fully connected layer to predict a coarse 3D point cloud. Then a 3D point-wise densification module ensures to lift the number of points of this point cloud using PointNetlike MLP. All models are trained end-to-end to minimise point set distances, e.g., Optimal Transport or chamfer distances.

Figure 2 .

 2 Figure2. For each scene, first row: 3D ground-truth and predictions for the RGB image according to Ad-aBins[START_REF] Shariq | Adabins: Depth estimation using adaptive bins[END_REF], BTS[START_REF] Han Lee | From Big to Small: Multi-scale local planar guidance for monocular depth estimation[END_REF], our Pix2Point VGG-chamfer, VGG-OT and ResNet-OT, all trained on 10k points. We follow the 3D representation of[START_REF] Caccia | Deep generative modeling of LiDAR data[END_REF]. Bird's eye view where the colour encodes the altitude. Second row: the input RGB image and the ground-truth-to-prediction error map for each method. Errors are from 0 (blue) to 50cm (red).

Table 1 .

 1 Comparison of 3D scene reconstructions on KITTI. We report completeness and accuracy. All methods are trained with 10k point clouds.

	Approaches	Completeness ↑(in %)	Accuracy ↓
		50cm	25cm	10cm	in m	rel.
	P2P-ResNet-OT 71.35 48.82 15.12	1.92	0.18
	P2P-VGG-OT	67.4	47.7	14.7	1.79	0.19
	P2P-VGG-C	64.4	36.0	8.0	0.85 0.05
	DensePCR	59.9	23.5	3.5	1.77	0.18
	BTS	67.59	31.29	6.28	1.23	0.06
	AdaBins	65.86	27.52	5.71	1.25	0.06

17th International Conference on Machine Vision Applications (MVA) Fully Online, July 25-27, 2021. © 2021 MVA Organization

http://www.cvlibs.net/datasets/kitti/eval_depth.php? benchmark=depth_prediction


Wilfried B Holzapfel
email: holzapfel@physik.upb.de

Stefan Klotz
email: stefan.klotz@upmc.fr

Coherent thermodynamic model for Ice Ih -A model case for complex behavior

New data on the variation of the thermal expansion of ice Ih with temperature at ambient pressure together with new evaluations of the bulk modulus and earlier data for the heat capacity provide the basis for a coherent thermodynamic modeling of the main thermo-physical properties of ice Ih over its whole range of stability. The quasiharmonic approximation with one Debye and 7 Einstein terms, together with explicit anharmonicity, represents the dominant contributions next to minor "anomalies" from hydrogen ordering and lattice defects. The model fits accurately the main features of all experimental data and provides a basis for the comparison with earlier determinations of the phonon density of states and the Grüneisen parameters.

I. INTRODUCTION

Ice Ih can be considered as a prototype of solids with anomalous thermo-physical properties not only due to the negative slope of its melting curve [START_REF] Bridgman | Water, in the liquid and five solid forms, under pressure[END_REF][START_REF] Feistel | A new equation of state for H 2 O ice Ih[END_REF][START_REF] Choukroun | Thermodynamic model for water and high-pressure ices up to 2.2 GPa and down to the metastable domain[END_REF] and its negative thermal expansion at low temperature [START_REF] Röttger | Lattice constants and thermal expansion of H 2 O and D 2 O ice Ih between 10 and 265 K[END_REF][START_REF] Buckingham | Thermal expansion of singlecrystal H 2 O and D 2 O ice Ih[END_REF][START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] but also due to its anomalous specific heat capacity [START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Leadbetter | The thermodynamic and vibrational properties of H 2 O and D 2 O[END_REF][START_REF] Sugisaki | Calorimetric study of the glassy state. iv. heat capacities of glassy water and cubic ice[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] showing a slow linear increase over a wide range in temperature (at ambient pressure) and a correspondingly strong variation of its thermal Debye temperature 8 , related to a very wide and strongly structured phonon density of states (PDOS) [START_REF] Renker | Lattice dynamics of hexagonal ice[END_REF][13][14][15] . Various first-principles calculations of thermophysical properties of ice Ih [16][17][18][19] have provided a solid background for a better understanding of both the approximations involved in theses calculations and for the basic physics of ice. However the results do not represent the thermophysical properties within the accuracy of the experimental data. Also, the two purely empirical description of the thermodynamics of ice Ih by different series expansions [START_REF] Feistel | A new equation of state for H 2 O ice Ih[END_REF][START_REF] Choukroun | Thermodynamic model for water and high-pressure ices up to 2.2 GPa and down to the metastable domain[END_REF] do not represent the thermophysical properties within the accuracy of the experimental data. In the first case 2 , the negative thermal expansion at low temperature was obviously considered as an experimental artifact and replaced by a positive interpolation down to low temperatures. In the other case [START_REF] Choukroun | Thermodynamic model for water and high-pressure ices up to 2.2 GPa and down to the metastable domain[END_REF] , the series expansion was restricted to temperatures above 180 K. On the other hand, the properties of ice Ih in wide ranges of pressure and temperature are of basic interest to glaciologists as well as for models of icy planets and moons. To the best of our knowledge no comprehensive thermodynamic model for the thermo-physical properties of ice Ih has been reported so far representing all the presently available experimental data within their given accuracy and providing the basis for the prediction of all theses properties over the entire range of existent of ice Ih.

The present study was motivated by the observation that although the use of two Debye and one Einstein term provided a very good description of the recent thermal expansion data [START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] , these three terms are completely insufficient for the representation of the heat capacity: It needs a better description of the PDOS by at least one Debye term and more than 2 Einstein terms, whereby the PDOS derived from the neutron scattering data [START_REF] Renker | Lattice dynamics of hexagonal ice[END_REF][13][14][15] provide good guides for the selection of the representative Einstein temperatures.

The present approach is intended to model the dominant thermo-physical properties of ice Ih. At the same time it helps for the evaluation of the smaller, but significant other contributions detected first in the heat capacity [START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF] . These were attributed to a "proton glass (PG) transition" [START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF] and are clearly visible in more recent volume expansion data [START_REF] Buckingham | Thermal expansion of singlecrystal H 2 O and D 2 O ice Ih[END_REF] where it was interpreted as arising from solitons, i.e. localized librational movements of H 2 O molecules.

Further anomalies were attributed to ferroelectric ordering, to the ice Ih-XI phase transition and to special relaxation effects at low temperatures 20,21 .

II. THE MODEL

Coherent thermodynamic models for all thermo-physical properties of solids are usually based on the formulation of the partition function which provides, being a fundamental function of the system, the Helmholtz free energy F (T, V, n, ...) 22,23 as a function of temperature T , volume V , particle number n, and often other parameters which we do not need in the present case. If we keep the particle number fixed to 1 mole, the state of the system is characterized just by the two state variables T and V , whereby V as well as all the other extensive functions, like the free energy F (T, V ), the internal energy U (T, V ), and the entropy S(T, V ), are represented in the following as molar quantities. The free energy F (T, V ) includes two parts, the energy of the ground state (or "cold part") E zT (V ), and thermal contributions F th (T, V ):

F (T, V ) = E zT + F th (T, V ) (1)
For convenience we include the zero-point energy of the phonons in the ground state energy E zT (V ). Correspondingly, the pressure P and the bulk modulus K are also separated in two parts::

P (T, V ) = ∂F (T, V )/∂V = P zT + P th (T, V ) (2)
K(T, V ) = -V ∂P/∂V = K zT (V ) + K th (T, V ) (3)
For the limited stability range of ice Ih, the ground-state energy and the corresponding pressure and bulk modulus at zero temperature are well represented by the Murnaghan forms 24 :

E zT (V ) = V K 0 K 0 (K 0 -1) ( V 0 V ) K 0 + K 0 -1 (4)
P zT (V ) = K 0 K 0 ( V 0 V ) K 0 -1 (5) K zT (V ) = K 0 ( V 0 V ) K 0 (6)
Thereby V 0 stands for the volume, K 0 for the bulk modulus and K 0 for its pressure derivative at zero temperature and ambient pressure.

The dominant part of the thermal free energy F th (T, V ) and of the thermal pressure P th (T, V ) of insulators comes from phonons and is usually described in the quasi-harmonic phonon approximation as an integral or a sum over contributions from Einstein-oscillators, whereby the volume-dependent frequencies ν(V ) of the oscillators are commonly replaced by corresponding Einstein-temperatures θ(V ) = hν(V )/k B . The free energy of one Einstein term (without its zero-point energy) is

f E (T, θ) = k B T ln(1 -exp(-θ/T )) (7)
whereby we use lower case letters for all the functions related to single modes in contrast to capital letters used for molar quantities. The corresponding formulas for the internal energy u E (T, θ), heat capacity c E (T, θ), thermal pressure p E (T, θ) and thermal contribution to the bulk modulus k E (T, θ) are given in the Appendix.

Convenient relations for the volume dependence of the Einstein temperature θ, its firstorder Grüneisen parameter γ and its second-order Grüneisen parameter Γ are provide by the form:

θ qh (V ) = θ • (V /V 0 ) γ(Γ-1) exp [γΓ • (1 -(V /V 0 ))] (8)
which gives

γ qh (V ) = -∂ln(θ qh (V ))/∂ln(V ) = γ [1 -Γ • (V /V 0 )] (9) Γ qh (V ) = ∂ln(γ qh (V ))/∂ln(V ) = ΓV / [V 0 -Γ • (V -V 0 )] (10)
with the ambient pressure (V = V 0 ) values θ, γ and Γ. For simple, monatomic, nonmetallic, crystalline materials, the Debye model 25 provides a good approximation for the phonon density of state (P DOS) in the form:

P DOS D (T, θ D ) = 3N T 2 /θ 3 D (11)
for N oscillators and a cutoff at the volume dependent Debye temperature θ D . The quasi-harmonic thermal free energy in this Debye model has then the form:

F D (T, θ D ) = θ D 0 P DOS D (τ, θ D )f E (T, τ )dτ ( 12
)
and the corresponding Debye free energy f D (T, θ D ), internal energy u D (T, θ D ) and heat capacity c D (T, θ D ) per mode are given in the Appendix.

In the more general case of crystalline solids with more than one atom per unit cell the Debye approximation remains reasonable for the low energy acoustic phonons, and the remaining phonon contributions are usually well represented by a few "representative"

Einstein terms. In this way one obtains for the quasi-harmonic contribution to the molar free energy:

F qh (T, V ) = N n i w i f i (T, V ) (13)
with N for the number of molecules per mole and the weights w i for the free energies per mode f i (T, V ), whereby the index i refers to the fact that each term with its volume dependence includes the parameters θ i , γ i , and Γ i . In addition, the index i=0 refers to the Debye form (A5) and i >0 to the Einstein form (7). In the same way one obtains the other quasi-harmonic contributions to internal energy, heat capacity, total pressure and bulk modulus as weighted sums over the individual terms.

The question is now, how many Einstein terms are needed to represent the experimental data within their given accuracy? Information on the phonon spectrum for ice Ih from the neutron scattering function 14 motivated Fortes 6 to evaluate his thermal expansion data at ambient pressure with two Debye and one Einstein term. If we use the form (A3) for the contribution of a single Debye or Einstein term to the thermal pressure together with the inverted form (5) for the thermal expansion at ambient pressure as a function of the thermal pressure together with the assumption that the thermal pressure is very small with respect to the bulk modulus and the volume dependence on the right hand side can be neglected, we get as a good first order approximation:

V qh (T ) = V 0 1 + N K 0 V 0 n i γ i w i u i (T, θ i ) (14)
If one would assume that the two transverse and one longitudinal acoustic branches are represented by the two Debye terms and the three degrees of freedom for the librational modes by one Einstein term, the corresponding weights 2:1:3 in the quasi-harmonic approximation result in a very poor description of the heat capacity. If one tries to fit the data for the heat capacity [START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Leadbetter | The thermodynamic and vibrational properties of H 2 O and D 2 O[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] with theses three terms one gets the weights 0.78, 3.10 and 4.23, but still with large systematic deviations indicating that one needs a model with more Einstein terms. However, for the following discussion of different forms for the thermal expansion we may just use these weights to represent the corresponding Grüneisen parameters γ i with the use of the values from Ref. [START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] V 0 =19.304 cm 3 /g and K 0 =10.713 GPa in the form ( 14). In the first line of Table I we reproduce just the corresponding data of Fortes' original fit 6 , whereby we present only the decimals needed for this comparison. The corresponding residues of the thermal expansion data of Fortes 6 from his fitted form are illustrated in fig. 1 together with the deviations of the other optimized forms from his original relation.

For comparison, we have refitted the parameters of the form ( 14) with the algorithm of Mathcad, which we used for all the present calculations, and noticed only a small decrease in the standard deviation σ (second row in Table I marked 2D1E) probably due to a different weighting scheme. The blue dash-dotted line in fig. 1 illustrates the very small difference of this fit. Since the original Debye function with their integral form considerably slows down all least square fits, we applied the simpler pseudo-Debye form [26][27][28][29] . All the refinements were performed with the pseudo-Debye function given in the Appendix.

We tested the pseudo-Debye form with a fit of the form ( 14) to the experimental data (third row in Table I marked 2PD1E and green dotted line in fig. 1) and obtained almost the same results in comparison with the previous fit, justifying the use of this more efficient form in all the following calculations. The form with one pseudo-Debye term and two Einstein terms (row 4 in Table I and thin pink dashed line in fig. 1) deteriorated very slightly the fit but provides still the expected T 4 -variation at low temperatures in contrast to the fit with three Einstein terms (row 5 in Table I and thin dashed-double-dotted curve in fig. 1).

This would be the simplest form for the representation of the thermal expansion coefficient, still with much smaller systematic deviations from the data than the previously proposed form [START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] . The deviation in the thermal expansion calculated by integration of Fortes' previous form is shown in fig. 1 as the grey line. One may notice that the Debye temperatures in the first and second column of Table I are significantly larger than the corresponding Einstein temperatures due to the fact that the Einstein temperatures correspond to the centers of gravity of the respective Debye spectrum.

Before we evaluate in detail all the available data for the volume expansion [START_REF] Röttger | Lattice constants and thermal expansion of H 2 O and D 2 O ice Ih between 10 and 265 K[END_REF][START_REF] Buckingham | Thermal expansion of singlecrystal H 2 O and D 2 O ice Ih[END_REF][START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] the question remains, how many Einstein terms do we need, to represent the heat capacity at ambient pressure within the limits of the experimental accuracy.

III. HEAT CAPACITY AT AMBIENT PRESSURE

In the fit of the experimental C v -data [START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Leadbetter | The thermodynamic and vibrational properties of H 2 O and D 2 O[END_REF][START_REF] Sugisaki | Calorimetric study of the glassy state. iv. heat capacities of glassy water and cubic ice[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] we use only the more accurate data sets [START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] and start with constraining the Debye term with respect to the data at the lowest temperatures [START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] as illustrated in the double logarithmic plot of fig. 2. [START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] By trial and error, the plots in fig. 2 provided good starting values for the refinement of 1 Debye and 7 Einstein terms with the results for the Debye temperature θ 0 and the seven Einstein temperatures θ i as well as for their weights w i given later in Table II together with the other parameters for this 1D7E-model. Thereby, a compromise had to be found in the number of terms, because a larger number of terms improves the quality of the fit but leads to physically unreasonable parameters like negative weight and large uncertainties due to stronger correlations. This anomaly around 100 K had been attributed to molecular disorder 10 and was studied in more detail recently by thermal dilatometry 5 and neutron scattering 6 . It should be noted that the original term "proton glass" 10 is misleading since the molecules are still intact.

What is presumably meant is a glass behaviour in the molecular dipole orientations as the reorientation dynamics becomes extremely slow at these temperatures. It is visible in fig. 3 as a strong peak with significant scattering of the data. Since no theory exists from which an analytical expression for the associated heat capacity could be derived, we simply describe it empirically by:

C P G (T ) = R 2w C P G ∆T C P G T -T C P G ∆T C P G exp - T -T C P G ∆T C P G ) 2 (15)
with the parameters w C P G = 0.214, T C P G = 104 K and ∆T C P G = 22 K and represented in fig. 3 by the light blue line. Since large hysteresis and dependence on heating rate was observed for this transition 10 as illustrated by the green and blue crosses for the two different runs, we have excluded the range marked by the pink solid line in this fit. This transition may be related also to the disorder in the structure of ice Ih observed at first as an anomaly in the entropy 30 and studied in more detail in recent years 16,18,[31][32][33][34] . The corresponding disorder entropy is not included in the present quasi-harmonic approach, but could be added easily later, and it is not clear whether there are other contributions from the disorder to the heat capacity which are included here in the quasi-harmonic terms.

Slightly before melting a small additional contribution shows up as "premelting anharmonicity" (thin red line) related also to the occurrence of defects 35,36 , which we model by:

C a (T, T a ) = w a R T a T (T a -T ) -2 (16)
with the corresponding free energy

F a (T, T a ) = w a R ln(1 - T T a ) ( 17
)
and the internal energy

U a (T, T a ) = w a R T a T T a -T + ln(1 - T T a ) ( 18
)
with weight w a = 8.0•10 -5 and T a = 280 K.

Some problems with the quasi-harmonic approximation are noticed when one uses the experimental data [START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] to calculate the corresponding calorimetric Debye temperature θ C (T ) [START_REF] Leadbetter | The thermodynamic and vibrational properties of H 2 O and D 2 O[END_REF]38 shown in fig. ultrasonic measurements of the elastic coefficients 37 . This difference may be typical for the quantum effects in glasses at low temperatures 39 which lead to a linear contribution to the heat capacity in the milli-Kelvin region but would not effect the acoustic Debye temperature. From this point of view it would be interesting to perform a detailed study of the heat capacity in this very low temperature region, which would than allow to refine the present model further for this region. Due to the strong weight of the many data points below 5 K the fitted curve represents these data very well but deviates significantly in the region around 10 K. This deviation could be modeled by an additional Einstein term. However as long as we have no quantitative information on the possible quantum glass contribution this rather small difference (in absolute values) will be considered as an uncertainty of the present model.

IV. THERMAL EXPANSION AT AMBIENT PRESSURE

In the next step we want to find good starting values for the first-order Grüneisen parameters γ i (Table II) for the present 1D7E-model with its single Debye and 7 Einstein terms from fits of the form ( 14) to various data of the thermal volume expansion. Although the dilatometry data have a much higher precision than the diffraction data at low temperatures, they needed a 6 ppm decrease in the absolute values to match the diffraction data at 0 K. The differences over the full range in temperature are shown in fig. 6.

One may notice that the deviation of the dilatometry data from the later diffraction data is only 6 ppm below 50 K where the diffraction data provide a better reference for V 0 . More significant are the differences at higher temperatures (fig. 6). The dilatometry data show some small wiggles which represent additional deviations from the 1PD7E-model, though in the diffraction data 6 similar wiggles are not observed or resolved. In addition, the large difference between the dilatometry and the diffraction data in the upper temperature range needs some special attention.

On should recall that the dilatometric method probes macroscopic properties (changes in the length along crystallographic axes and hence the volume of the sample), whereas

the diffraction technique provides microscopic information on the crystal structure and the corresponding unit cell. Differences at elevated temperature have been attributed in many cases to the thermal excitation of defects. A difference of 1% just before melting seems not to be too large. However, in general, defects are expected to increase the volume 35,36 . Maybe this is an other anomaly of ice Ih.

With all these uncertainties we use just a reasonable average curve as reference for our model as indicated by the black curve in fig. 6. This reference fits perfectly the dilatometry data in the lower temperature range and gives a reasonable average at higher temperatures.

It corresponds to the 1PD7E-form with just a small change in γ 5 given in Table II One may notice in Table II that the largest differences in the γ j occur for the first Debye term and the last Einstein term. This reflects the fact that, on the one hand, significant differences in the thermal expansion data are observed at very low temperature (fig. 5, giving the variation at low temperatures on double logarithmic scale) and, on the other hand, in the upper temperature range (figure 6, presenting the differences with respect to the best fit of eq. 14 to the dilatometry data 5 ). The less negative values of γ 0 in this case fit also better to the mostly positive pressure dependencies of the elastic coefficients C ij 40,41 and to the dispersion of the mode-Grüneisen parameters measured by neutron scattering under pressure 15 .

V. BULK MODULUS AT AMBIENT PRESSURE

When one compares the data for the temperature dependence of the bulk modulus at ambient pressure [START_REF] Feistel | A new equation of state for H 2 O ice Ih[END_REF][START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF]42,43 in fig. 7, one can notice at first that there is a reasonable agreement in the middle of the given temperature range. In the upper temperature range the differences remain small compared with the uncertainty of the data. However, at low temperatures the differences are more significant. The evaluations [START_REF] Feistel | A new equation of state for H 2 O ice Ih[END_REF][START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] were primarily based on a complete set of elastic constants determined by ultrasonic measurements 42 covering the range from 60 K to 110 K, various data at fixed temperatures in the range 240-260 K 40,[44][45][46] and suitable interpolation schemes between. From this point of view extrapolations [START_REF] Feistel | A new equation of state for H 2 O ice Ih[END_REF][START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] to temperatures below 60 K are not based on experimental data. The most recent study 43 , on the other hand, covers the range from 50 K to 273 K and takes into account all the earlier measurements.

Therefore we give this evaluation twice as much weight when we fit all the data together Our best fit uses a value of K 0 =11.0 GPa and K 0 =5.5 and provides now also the initial values of the second-order Grüneisen parameters Γ i . This includes the inherent volume dependence of the Debye and Einstein temperatures in eq. 19 for the later refinement of the quasiharmonic free energy in the next section.

VI. THE COHERENT THERMODYNAMIC MODEL

For the fundamental function, the total free energy, we now use the form:

F (T, V ) = E zT + F qh (T, V ) + F ah (T, V ) (19)
where the ground-state energy E zT depends on the parameters V 0 , K 0 , K 0 . The quasiharmonic contribution F qh (T, V ) includes now the full volume dependence in the Debye and Einstein temperatures as well as in the Grüneisen parameters ( 8)- (10). The volume dependence in the explicit premelting anharmonicity ( 16)-( 18) is taken into account by the assumption that the characteristic temperature T a in theses relations is replaced by the volume dependent form:

T ah (V ) = T a (V /V 0 ) -γ ah (20)
providing the form:

F ah (T, V ) = F a (T, T ah (V )) ( 21
)
and the volume dependence in other thermodynamic relations (16,18).

The refinement of the 8 weights w i together with the 8 times three parameters θ i , γ i , and

Γ i of

VII. PREDICTIVE POWER OF THE PRESENT MODEL

Since the results of the global fit described in the previous chapter do not show any significant difference with respect to the results of the restricted fit to the heat capacity shown in fig. 3, this figure needs no duplication. Similarly, the differences in the fit of the volume expansion were marginal and we find that the additional constraint in the coherent model leads only to a very minor increase of 0.8% in the standard deviations. Without the contribution from the PG-transition the RMSD would increase by 3%. This indicates again that a small contribution from this transition can be noticed also clearly in the volume expansion derived from the neutron diffraction data 6 . The data for the bulk modulus at ambient pressure 6 are also fitted perfectly with no visible difference to the black fitted curve in fig. 7 and provides in this way the values for the refined second-order Grüneisen parameters given in Table IV.

With the coherent thermodynamic model one can now calculate all the thermo-physical properties of ice Ih from the free energy (19) for the whole range of stability form the initial volume V 0 = 19.30448 cm 3 /mol to the smallest volume V tr = 18.947 cm 3 /mol, both at 0 K. Since the effect of pressure on the heat capacity in this small range of compression is rather small, only the relative change is illustrated here in fig. 8.

One can notice here that this relative change becomes rather large only at low temperature due to the large negative values of the Grüneisen parameters for the low-frequency modes. The small maximum around 10 K can be related to the difference of the Grüneisen parameters for this range.

The effect of pressure on the heat capacity may be seen more clearly in fig. 9, which presents the change of the calorimetric Debye temperature by the red curve with respect to the variation at ambient pressure as blue base-line. In addition, fig. 9 shows also the deviation of the values for ambient pressure calculated from the experimental C v data [START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] with respect to the model curve for comparison with fig. 4, illustrating again the uncertainty below 10 K due to possible effects of disorder and around 100 K due to the PG-transition.

The red curve shows quantitatively at low temperatures the predicted decrease due to the negative Grüneisen parameters and at higher temperatures the increase due to the positive values for the higher phonon energies.

In comparison with the small changes of the heat capacity under pressure, the thermal FIG. 10. Isochores for the thermal pressure P th (V, T ) at V 0 , the volume at ambient pressure and 0 K. The two other values 0.99 V 0 and 0.98 V 0 cover roughly the range of stability of ice Ih (pink, green, blue curves) with the variation along the melting curve (solid red line), corresponding to the free energy ( 19) with the parameters of Table II. The contribution from the PG-transition is too small to be noticed here.

pressure shows much stronger changes as illustrated in fig. 10.

Isochores for the isothermal bulk modulus for the same range up to the melting curve are presented in fig. 11, which does not show any anomalous behavior.

The equation of states (EOS) can be calculated from ( 19)

P (T, V ) = ∂F (T, V )/∂V = P zT + P qh (T, V ) + P ah (T, V ) (22)
either in the form of isochors or by inversion in the form of isobars V (T, P ). Corresponding to the two terms of the thermal contribution to the free energy, the thermal pressure includes also a contribution from the premelting anharmonicity

P ah (T, V ) = γ ah V U ah (T, V ) (23)
and one gets for the isobaric volume the implicit form Isochores for the isothermal bulk modulus K(V, T ) at V m0 , 0.99 V m0 and 0.98 V m0 cover roughly the range of stability of ice Ih (pink, green, blue curves) and the variation along the melting curve (red line). V m0 is the volume at the melting point at zero pressure.

V (T, P ) = V 0 1 - K 0 K 0 (P th (T, V (T, P )) -P ) 1/K 0 -∆V pg (T ) (24)
which can be solved with our reference values K 0 =11 GPa and K 0 =6.6 by iteration to get the three isobars for 0, 100 and for 213 MPa covering the range up to the first high pressure phases as shown in fig. 12. Numerical values are given in the appendix V.

One may notice that the negative thermal expansion is barely visible and the effect of the PG-transition can not be seen on this scale in fig. 12. More interesting are therefore the isobars for the thermal expansion coefficient, α(T, P ), illustrated in fig. 13.

Here one can notice the similarity to the thermal pressure illustrated in fig. 10, with, however, the region of negative values shifted to smaller temperatures, the expected decrease under compression along the melting curve and larger negative values at low temperatures.

More interesting are the isotherms for the isothermal bulk modulus in fig. 14. The isothermal variations look here almost like straight lines, but the two isobars for ambient pressure and the transition pressure of 213 MPa illustrate the anomalous behavior of ice Ih.

On the other hand, the isotherms P (V, T ) (fig. 15) for the total pressure show more directly the anomaly of ice Ih by the fact that the 0 K isotherm overlaps with the 50 K isotherm and the 100 K isotherm appears to be shifted downwards.

The present model allows also to calculate the calorimetric Grüneisen parameter:

γ cal (T, V ) = α(T, V )V K(T, V )/C V (T, V ) (25)
and its isobaric variation. Two isobars are represented in fig. 16, one for ambient pressure and one for the upper boundary of ice Ih, at 210 MPa. contributions from intrinsic anharmonicity not separately modeled in this coherent approach which did not yet take into account the additional information provided by studies on the PDOS to be discussed in the next section.

VIII. COMPARISON WITH PDOS

When we compare the results of the present model with experimental data for the PDOS of ice Ih 13 from neutron scattering, we can see in fig. 19 that the selected Einstein terms reproduce reasonably well the dominant features of the experimental PDOS.

Since the original data for the PDOS of ice Ih from neutron scattering 13 where not available in numerical form, we digitized the corresponding part of the original fig. 2 13 and scaled the area to 3 modes per molecule, whereby we did not correct for the missing part below 30 K. Integration of this PDOS using the standard formula gave the C v -data shown in fig. 20 (blue dots). For this purpose, the background of the PDOS was reduced by 3% with respect to the value of the first peak to obtain a better fit of C v below 100 K. The comparison with the present coherent model based on formula (19) and using only the first 6 terms (0 to 5 in Table IV) (green dashed line) shows a slight deviation above 100 K which is attributed to anharmonic contributions of these translational modes. An accurate estimate of the anharmonicity of the librational modes could only be made by a detailed comparison with not yet available data for the librational PDOS. However, for the . It provides some additional information on the PDOS since it represent its typical features but not the correct relative intensities.

One can see clearly in fig. 21 that the high energy part of the theoretical PDOS 15 shows a too small TO-LO-splitting and a total shift to lower energies in comparison with the corresponding structures in both the experimental PDOS 13 and the scattering function 14 . In the low energy part a significant shift is most clearly visible around 150 K. Theses differences illustrate that the lattice dynamical model [START_REF] Renker | Lattice dynamics of hexagonal ice[END_REF]15 helps very much in understanding the structure of the PDOS, but does not yet provide accurate data for the description of a thermodynamical model covering the whole energy range. In fact this point had already been indicated by the author of the model [START_REF] Renker | Lattice dynamics of hexagonal ice[END_REF] . Also, the more recent theoretical models 17,19 provide primarily a better understanding of the physical background, and the comparison of their results with the experimental data serves as a guide for the choice of the best approximations to be used in theses models.

So far we have compared our results only with the PDOS data for the translational modes (below 500 K, see fig. 22). For higher energies (fig. 22) only the neutron scattering function is available which cannot be converted easily into a PDOS but provides at least the positions of its characteristic features described by the θ 6 and θ 7 terms of the present model.

Since we expect in the strictly harmonic approximation for the librational modes a total weight of 3 per molecule, we see that this value is nearly reached already by the θ 6 term, and the large weight of the θ 7 term must be attributed mostly to anharmonic contributions.

So far we have used what we call a quasi-harmonic approximation which models the anharmonicity just by excess weights of the quasi-harmonic terms. In this way we have limited the number of free parameters of the present approach. However, more correctly we should have used a quasi-harmonic approach with the correctly constrained weights (or number of modes per molecule) together with an implicit anharmonic contribution for each mode. This implicit anharmonicity is commonly related to a temperature dependence of the quasi-harmonic frequencies at constant volume [START_REF] Leadbetter | The thermodynamic and vibrational properties of H 2 O and D 2 O[END_REF]28,29,47 and can be approximated in first order by an anharmonic contribution proportional to the internal energy of the mode i, scaled by a factor β i , i.e. (1

+ β i u i (T, θ i )/(k B θ i ))
. When we applied this scheme to the representation of the heat capacity, we noticed first of all that the θ 7 term had to be shifted into the region of the librational modes together with a smaller shift of the θ 6 term and a dramatic change in the weights of these two terms as indicated by the values given in Table III with the new harmonic weights w i and the additional anharmonicity parameters

β i .
One may notice also in Table III that we need no anharmonic contribution for the first two terms and have to change for the translational modes only the weight of the highest term θ 5 . Since the anharmonicity parameters of the librational modes are strongly correlated with the values of the representative frequencies, we did not attempt to refine the parameters perfectly. In fact, the representative Einstein terms should be determined with respect to a reliable PDOS, which is not yet available. Any extra volume dependence of the anharmonicity parameters has been neglected here and could be handled by an increase in the number of free parameters which could not be constrained by any experimental data at the present time. From this point of view fig. 23 shows now only that this quasi-harmonic approach with the parameters of Table IV fits the experimental data for the specific heat almost as good as the refined unconstrained model. But the distribution of the anharmonicity is physically more reasonable in this constrained approach and allows therefore also for a better estimate of the remaining (systematic) uncertainties.

In any case we can notice that the anharmonic contributions to the translational modes are very moderate and their pressure dependence does not introduce significant uncertainties in any of the thermodynamic relations at moderate temperatures (T < 100). The effect of pressure on the anharmonicity of the librational modes could be positive or negative, and even as a small contribution it represents now the larges uncertainty in the prediction of (black open diamonds), values from the present model (row 6 in Table II) based on the reference data for the bulk modulus (pink diamond with pink connecting line). The other data from Table II, are shown as small green, red, and blue open diamonds respectively, to illustrate the uncertainty related to the different literature data for the thermal expansion [START_REF] Röttger | Lattice constants and thermal expansion of H 2 O and D 2 O ice Ih between 10 and 265 K[END_REF][START_REF] Buckingham | Thermal expansion of singlecrystal H 2 O and D 2 O ice Ih[END_REF][START_REF] Fortes | Accurate and precise lattice parameters of H 2 O and D 2 O ice Ih between 1.6 and 270 K from high-resolution time-of-flight neutron powder diffraction[END_REF] .

allowed us to extract a few values for this energy dependent average Grüneisen parameter from the shifts of some prominent features in theses spectra, which are presented in fig. 24.

At first one may notice reasonably good agreement between the lattice dynamical data 15

and the present data at low energies (below 100 K) and than significantly smaller absolute values with respect to the computer simulation 17 . The values for the librational modes with phonon energies between 800 to 1700 K are obviously not well constrained by the present model due to the influence of effects from anharmonicity and lattice defects.

X. DISCUSSION

Although ice Ih appears at first to be a simple molecular crystal, it shows a very complex behavior which is not yet completely understood. It is clear that the phonon spectrum is separated into 4 distinct regions with the translational modes at the lowest energy (up to 400 K in temperature units), the librational modes from about 700 to 1500 K, the bending modes from 2300 to 2600 K and the stretching modes around 4900 K. However, good values for the corresponding PDOS are only available for the translational modes which leads to the problem that the anharmonic contributions to the thermodynamic properties are well constrained only in the low temperature region (< 150 K) and their unknown pressure dependence limits the accuracy of the present modeling at higher temperatures and elevated pressures. At very low temperatures (< 10 K) possible effects from a quantum glass behavior limit the accuracy of the present thermodynamic model. Furthermore, the lack of a good model for the "proton glass transition" does not yet allow an accurate formulation of the entropy. This effects also the evaluation of the thermal expansion and the related Grüneisen parameters, especially due to the fact, that the small but significant difference between the microscopic and macroscopic determinations of the volume expansion is not yet well understood. Nevertheless, the evaluation of all the experimental data for the present thermodynamic model can be considered as a good basis for the selection of the best approximations to be used in theoretical models which have provided so far a better understanding of the physical phenomena relevant for ice Ih, but still with limited accuracy as far as a comparison with experimental data is concerned. Cold ice remains a hot topic!

FIG. 1 .

 1 FIG.1. Relative deviations of the volume data of Fortes 6 in ppm (red and pink diamonds) with respect to his best fit with 2 Debye and 1 Einstein term as black base line. Additional curves show the small deviations of the present refit with 2 Debye terms and 1 Einstein term as blue dash-dotted line and the fits with 2 pseudo-Debye and 1 Einstein term (green dotted curve), with 1 pseudo-Debye and 2 Einstein terms (dashed pink curve) and with 3 Einstein terms only (orange dash-doubledotted curve), respectively. Integration of Fortes' analytic form for the thermal expansivity (grey solid curve, "int-alfa") shows significantly larger deviations.

  FIG.2. a) Contributions of the Debye term (i = 0) and the 7 Einstein terms to the heat capacity (in units of R) at ambient pressure in log-log scale versus temperature and b) Sum of the contributions to the heat capacity on linear scales. Symbols are experimental data from refs.[START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF]

FIG. 3 .

 3 FIG.3. Experimental heat capacity data7,10,11 (upper) and deviation from the fitted relation for C v (lower) including contributions from the PG-transition (light blue) and explicit premelting anharmonicity (red). The anomalous region (thick pink) has been excluded in the fit.

Fig. 2

 2 Fig.2illustrates that the C v -data deviate already at temperatures below 10 K from the normal Debye behavior and the weight of the Debye term given in TableIIis though very small. The almost linear increase of C v over a wide range in temperature and the large number of Einstein terms indicate here again that the PDOS must deviate strongly from a simple Debye form. The quality of the fit is illustrated in fig.3by the residues with respect to the quasi-harmonic contribution, whereby we excluded in the fit all the data in the region of anomalous behavior around 100 K, the region marked in pink in fig.3.

4 .

 4 FIG. 4. Comparison of the acoustic θ A (blue diamond) 37 and the calorimetric Debye temperature θ C (circles and crosses) below 30 K with the present quasi-harmonic curve (red line) and a rough fit (blue line) illustrating possible quantum glass contributions.

  FIG. 5. Relative volume change in the low temperature region on double logarithmic scale. The dilatometry data 5 are represented by the blue stepped curve with the blue squares as smoothed values, and the solid light blue curve corresponding to the present best 1PD7E-fit. The green curve reproduces the original best fit of earlier diffraction data 4 and the red squares together with their best 1PD7E-fit (red curve) illustrate the later diffraction data 6 . The dashed orange line stands for the expected T 4 -behavior at low temperature.

FIG. 6 .

 6 FIG. 6. Relative volume change (upper) and deviation (lower) from the dilatometry data 5 (blue squares with the blue line being a 1PD7E fit). X-ray diffraction data 4 are represented by the green curve and neutron diffraction results 6 by the red dots and the red curve being the corresponding 1PD7E fit. The black curve represents an estimated average used in this work as reference data.

FIG. 7 .

 7 FIG. 7. Comparison of various reference data for the temperature dependence of the bulk modulus at ambient pressure. The reference data from 2006 2 are represented by orange open diamonds. A later evaluation 6 based on published low-temperature ultrasonic data and Brillouin results at higher temperatures is shown as red line. The most recent evaluation 43 is marked by the blue line with the blue diamonds. Our best fit, discussed in the text, is represented by the black line. The green dashed line illustrates the cold lattice contribution.

FIG. 8 .

 8 FIG. 8. Relative change of the heat capacity by 2% compression corresponding to the complete range of stability of ice Ih.

∆

  FIG.9. Effect of pressure on the calorimetric Debye temperature (red curve) with respect to variation at ambient temperature taken as blue base line. In addition, the deviation of the values calculated from the experimental C v -data for ambient pressure[START_REF] Flubacher | Heat Capacity of Ice at Low Temperatures[END_REF][START_REF] Haida | Calorimetric study of the glassy state X. Enthalpy relaxation at the glass-transition temperature of hexagonal ice[END_REF][START_REF] Smith | Heat capacities and thermodynamic functions of hexagonal ice from T = 0.5 K to T = 38 K[END_REF] (symbols) are illustrated for comparison with fig. 4.

  FIG. 11. Isochores for the isothermal bulk modulus K(V, T ) at V m0 , 0.99 V m0 and 0.98 V m0 cover

FIG. 12 .

 12 FIG.12. Isobars for the molar volume V (T, P ) for P = 0, 100 and 213 MPa (pink, green, blue curves) and the variation along the melting curve (red line).

  FIG. 13. Isobars for the thermal expansion coefficient α(T, P ) for P = 0, 100 and 213 MPa (pink, green, blue curves) and the variation along the melting curve (red line).

  FIG.14. Isobaric variation of the isothermal bulk modulus at ambient pressure (pink curve) and along the upper boundary of ice Ih at 213 MPa (dark green curve) together with the isotherms for 0 K (blue line), 50 K (turquoise dotted), 100 K (green dotted), 150 K (orange dotted), 200 K (pink dotted), 250 K (red dotted) and along the melting line (red line).

Fig. 20 (

 20 Fig. 20 (upper) shows at first that in the range below 100 K the heat capacity (red curve) is reproduced very well by values calculated from the experimental PDOS 13 (blue points).

  FIG. 16. Comparison of the calorimetric Grüneisen parameter from the present model for ambient pressure (pink curve) and for the upper boundary of ice Ih at 210 MPa (blue dashed curve).

  FIG. 17. Volume dependence of the Debye (θ 0 ) and Einstein (θ 1 -θ 7 ) temperatures.

5 FIG 4 and θ 5 .

 545 FIG. 18. Volume dependence of the Grüneisen parameters for the Debye (θ 0 ) and Einstein (θ 1 -θ 7 ) temperatures of the present model.

  FIG. 20. Upper: Heat capacity of ice Ih at ambient pressure (red line) compared with values calculated from the experimental PDOS 13 (blue dots) including only the translational modes as shown in fig. 19 and from the present coherent model (green dashed line) using only the first 6 terms (0-5). An estimate of the anharmonicity in the translational modes results in the green solid line and the blue line illustrates the total harmonic contribution including the librational modes. Lower: Difference between the quasi-harmonic heat capacity obtained by the experimental PDOS 13 (dotted blue line in the upper plot) and values obtained with different distributions of the translational anharmonicity. The green curve corresponds to an equal distribution over all the translational terms (θ 0 to θ 5 ), the pink curve to anharmonicity only on the θ 5 term, and the orange curve to an equal distribution over θ 4 and θ 5 . The gray curve represents the surplus of the heat capacity calculated from the theoretical PDOS 15 .

  FIG. 21. Comparison of the experimental PDOS 13 (green curve) with the theoretical PDOS 15 (blue curve) and the neutron scattering function 14 (red curve with roughly adjusted scale).

  FIG. 22. Comparison of the neutron scattering function 14 (gray curve) with Debye and the Einstein terms of the present model (red and blue crosses represent weights). The weight w 7 = 6.25 is out of scale and only represented by the area of the corresponding Gaussian at θ 7 = 1752 K together with Gaussians for the other terms (red curve).

  FIG. 23. Deviations of the experimental data for the heat capacity from the representation by the "constrained" quasi-harmonic model with explicit anharmonic contributions. The contribution from the PG-transition has been subtracted like in fig.3together with the extra premelting anharmonicity, which is indicated here by the light blue line.

x 2

 2 T, θ, γ, Γ) = (1 + γ E -Γ E ) • p E (T, θ)term, the free energy f D (T, θ D ), internal energy u D (T, θ), heat capacity c D (T, θ), are given by:f D (T, θ D ) = 3k B θ D ln(1 -exp(-x))dx (A5) u D (T, θ D ) = 3k B θ D T θ D

TABLE I .

 I Best fitted parameters for Debye or Einstein temperatures θ i and Grüneisen parameters γ i together with normalized standard deviations σ for the representation of the thermal expansion data of Fortes 6 with various combinations of Debye, pseudo-Debye and Einstein terms labeled with D, PD and E for the terms used in the form (14).

	M odel	θ 1 (K)	θ 2 (K)	θ 3 (K)	γ 1	γ 2	γ 3	σ
	Fortes	87	513	1254	-0.70	1.32	1.33	1.00
	2D1E	87	513	1254	-0.69	1.32	1.33	0.90
	2PD1E	86	513	1239	-0.69	1.32	1.28	0.90
	1PD2E	65	376	1079	-0.43	1.18	1.59	1.08
	0PD3E	41	376	1211	-0.40	1.16	1.37	1.00

  . The corresponding parameters for the preliminary description of the experimental data by theses refinements are presented in Table II as "Our ref.".

TABLE II .

 II Preliminary parameter values for the present 1PD7E-model from independent fits of the data to the heat capacity, thermal expansion and bulk modulus at ambient pressure. γ a,i , γ b,i , and γ c,i correspond to fits to data from refs. 4 , 6 , and 5 , respectively. "Our ref." refers to the parameter set based on the black line in fig. 6.

		0	1	2	3	4	5	6	7
	θ i (K)	55	55	104	227	324	435	955	1693
	w i	0.044	0.142	0.844	0.566	0.626	1.041	2.596	5.94
	γ a,i	-0.71	-1.91	-0.56	0.72	2.51	2.11	0.63	1.12
	γ b,i	-1.41	-1.78	-0.25	0.94	2.22	1.77	0.76	0.65
	γ c,i	-0.61	-1.52	-0.34	0.94	1.94	1.35	1.41	0.09
	Our ref.	-0.61	-1.52	-0.34	0.94	1.94	1.54	1.41	0.09
	Γ i	-18.5	-18.5	-92.9	3.80	8.49	-8.25	-0.93	-0.93

TABLE III .

 III Parameters for the "constrained" quasi-harmonic model with implicite anharmonicity.

		0	1	2	3	4	5	6	7
	θ i (K)	55	55	104	227	340	407	830	1160
	w i	0.044	0.141	0.849	0.576	0.578	0.813	0.300	2.700
	β i	0.0	0.0	0.03	1.14	0.17	0.28	23.00	3.40
	function 14

TABLE V .

 V Calculated values of the specific heat C V at ambient pressure (in units of R) and the relative volume change ∆V /V 0 at ambient pressure, 100 MPa and 213 MPa, with 0 K molar volumes of V 0 =19.30441 cm 3 /mole, 19.13392 cm 3 /mole, and 18.97276 cm 3 /mole, respectively.

	T (K)	CV /R	∆V /V0 • 10 4	∆V /V0 • 10 4	∆V /V0 • 10 4
		0 MPa	0 MPa	100 MPa	213 MPa
	0	0	0	0	0
	5	0.0030	-0.0013	-0.0015	-0.0018
	10	0.0354	-0.0391	-0.0486	-0.0588
	15	0.1174	-0.2065	-0.2528	-0.2991
	20	0.2393	-0.5198	-0.6268	-0.7302
	25	0.3732	-0.9463	-1.1300	-1.3033
	30	0.5018	-1.4486	-1.7224	-1.9763
	35	0.6216	-1.9903	-2.3696	-2.7166
	40	0.7341	-2.5348	-3.0389	-3.4952
	45	0.8420	-3.0441	-3.6962	-4.2819
	50	0.9471	-3.4803	-4.3066	-5.0445
	55	1.0503	-3.8072	-4.8353	-5.7498
	60	1.1519	-3.9925	-5.2501	-6.3663
	65	1.2520	-4.0088	-5.5230	-6.8652
	70	1.3501	-3.8342	-5.6308	-7.2223
	75	1.4460	-3.4526	-5.5555	-7.4182
	80	1.5394	-2.8525	-5.2839	-7.4382
	85	1.6300	-2.0270	-4.8073	-7.2717
	90	1.7179	-0.9722	-4.1202	-6.9117
	95	1.8031	0.3128	-3.2200	-6.3542
	100	1.8858	1.8273	-2.1062	-5.5973
	105	1.9661	3.5692	-0.7798	-4.6411
	110	2.0443	5.5355	0.7573	-3.4867
	115	2.1207	7.7228	2.5023	-2.136
	120	2.1955	10.1274	4.4521	-0.5910
	125	2.2691	12.7459	6.6037	1.1450
	130	2.3415	15.5748	8.9539	3.0693

θ D /T

ACKNOWLEDGMENTS

The authors whish to thank J.J. Neumeier (Montana State University, Bozeman, USA), Th. Strässle (Paul Scherrer Institute, Villigen, Switzerland) and D. Fortes (Rutherford Appleton Laboratory, Chilton, U.K.) for very helpful discussions.

Appendix A: Thermodynamic relations

For one Einstein term, the form of the internal energy, heat capacity, thermal pressure and isothermal bulk modulus are given by:

The formula for the thermal pressure p D (T, θ, γ) and thermal contribution to the bulk modulus k D (T, θ, γ, Γ) correspond to the forms (A3) and ( A4) with E replaced by D.

Appendix B: Pseudo-Debye forms

The free energy of the pseudo-Debye form consists of a leading term f P D 0 for low temperatures:

and two Einstein terms f E in the form of:

with the weights w 0 =0.05652, w 1 =0.90748, w 2 =0.036 and the scaling factors f 1 =0.81 and

.41 for the effective (fixed) Einstein temperatures. Since all the other thermodynamic relations for the pseudo-Debye terms are similar sums, only the expressions involving the leading term f P D 0 are shown here for the internal energy:

and for the heat capacity:

which shows that these terms have the same asymptotic behavior at low and high temperatures as the normal Debye form.

Appendix C: Modelling the "proton glass transition"

When we subtract the contribution of the PG-transition in the volume expansion, we use a Gaussian form:

with the weight w V pg = 0.00014, the position T V pg = 103 K and the width ∆T V pg = 14 K in the fits of the diffraction data which reduces the standard deviation of these by 6%.

But with all the uncertainty related to this sluggish feature, we keep this term as a small correction in the thermal expansion and not as a normal contribution to the equilibrium thermodynamics.

For any application of the present model we give here all the parameters of the final refinement in Table IV.