Romain Xu-Darme

Georges Quénot

Zakaria Chihani

Marie-Christine Rousset

PARTICUL: Part Identification with Confidence measure using Unsupervised Learning

Keywords: Part detection, Unsupervised learning, Interpretability, Confidence measure, Fine-grained recognition

In this paper, we present PARTICUL, a novel algorithm for unsupervised learning of part detectors from datasets used in fine-grained recognition. It exploits the macro-similarities of all images in the training set in order to mine for recurring patterns in the feature space of a pre-trained convolutional neural network. We propose new objective functions enforcing the locality and unicity of the detected parts. Additionally, we embed our detectors with a confidence measure based on correlation scores, allowing the system to estimate the visibility of each part. We apply our method on two public fine-grained datasets (Caltech-UCSD Bird 200 and Stanford Cars) and show that our detectors can consistently highlight parts of the object while providing a good measure of the confidence in their prediction. We also demonstrate that these detectors can be directly used to build part-based fine-grained classifiers that provide a good compromise between the transparency of prototype-based approaches and the performance of non-interpretable methods.

Introduction

With the development of deep learning in recent years, convolutional neural networks (CNNs) have quickly become the backbone of all state-of-the-art visual recognition systems. CNNs are highly efficient but also highly complex systems manipulating abstract (and often opaque) representations of the image -also called feature vectors -to achieve high accuracy, at the cost of transparency and interpretability of the decision-making process. In order to overcome these issues, a solution, explored in [1, 3, 6, 7, 10-12, 17, 19, 22, 27, 34, 41, 44], consists in building an intermediate representation associating each input image with a set of semantic attributes. These attributes, usually representing an association between a part of an object (e.g., head of a bird) and a property (e.g., shape, color), can be either used by interpretable methods [START_REF] Chen | This looks like That: Deep learning for interpretable image recognition[END_REF][START_REF] Nauta | Neural Prototype Trees for Interpretable Finegrained Image Recognition[END_REF][START_REF] Rymarczyk | ProtoPShare: Prototype Sharing for Interpretable Image Classification and Similarity Discovery[END_REF] to produce the decision, as post-hoc explanation [START_REF] Hassan | Explaining Visual Classification using Attributes[END_REF][START_REF] Hendricks | Generating Visual Explanations[END_REF] of a particular decision, or as supplementary information in order to discriminate similar categories in the case of fine-grained visual classification (FGVC) [START_REF] Fukui | Attention Branch Network: Learning of Attention Mechanism for Visual Explanation[END_REF][START_REF] Han | Attribute-Aware Attention Model for Finegrained Representation Learning[END_REF][START_REF] Liang | A Unified Multiplicative Framework for Attribute Learning[END_REF][START_REF] Zhao | Recognizing Part Attributes With Insufficient Data[END_REF][START_REF] Zhu | Multi-label Convolutional Neural Network Based Pedestrian Attribute Classification[END_REF]. In practice, attributes are learned through fully supervised algorithms [START_REF] Abdulnabi | Multi-task CNN Model for Attribute Prediction[END_REF][START_REF] Farhadi | Describing Objects by their Attributes[END_REF][START_REF] Fukui | Attention Branch Network: Learning of Attention Mechanism for Visual Explanation[END_REF][START_REF] Han | Attribute-Aware Attention Model for Finegrained Representation Learning[END_REF][START_REF] Liang | A Unified Multiplicative Framework for Attribute Learning[END_REF][START_REF] Wang | Attribute Recognition by Joint Recurrent Learning of Context and Correlation[END_REF][START_REF] Zhao | Recognizing Part Attributes With Insufficient Data[END_REF][START_REF] Zhu | Multi-label Convolutional Neural Network Based Pedestrian Attribute Classification[END_REF], using datasets with hand-crafted annotations. Such datasets are expensive to produce -using expert knowledge or online crowd-sourcing platforms -and prone to errors [START_REF] Jo | Lessons from Archives: Strategies for Collecting Sociocultural Data in Machine Learning[END_REF][START_REF] Northcutt | Pervasive label errors in test sets destabilize machine learning benchmarks[END_REF]. Therefore, a lot of effort [START_REF] Chen | This looks like That: Deep learning for interpretable image recognition[END_REF][START_REF] Ding | Selective Sparse Sampling for Fine-Grained Image Recognition[END_REF][START_REF] Ge | Weakly Supervised Complementary Parts Models for Fine-Grained Image Classification From the Bottom Up[END_REF][START_REF] Han | P-cnn: Part-based convolutional neural networks for fine-grained visual categorization[END_REF][START_REF] Korsch | Classification-Specific Parts for Improving Fine-Grained Visual Categorization[END_REF][START_REF] Li | Attribute Mix: Semantic Data Augmentation for Fine Grained Recognition[END_REF][START_REF] Peng | Object-Part Attention Model for Fine-grained Image Classification[END_REF][START_REF] Sicre | Unsupervised Part Learning for Visual Recognition[END_REF][START_REF] Yang | Learning to Navigate for Fine-grained Classification[END_REF][START_REF] Zhang | Unsupervised Part Mining for Finegrained Image Classification[END_REF][START_REF] Zheng | Learning Multi-attention Convolutional Neural Network for Fine-Grained Image Recognition[END_REF][START_REF] Zhou | Human-Understandable Decision Making for Visual Recognition[END_REF] has been recently put into the development of more scalable techniques using less training information. In particular, the task of localizing different parts of an object, a prerequisite to attribute detection, can be performed in a weakly supervised (using the category of each image) or unsupervised manner. By focusing on general features of the objects rather than discriminative details, part detection represents an easier task than attribute detection. However, it must offer strong evidence of accuracy and reliability in order to constitute a solid basis for a trustworthy decision.

In this paper, we present PARTICUL 3 (Part Identification with Confidence measure using Unsupervised Learning), a plug-in module that uses an unsupervised algorithm in order to learn part detectors from FGVC datasets. It exploits the macro-similarities of all images in the training set in order to mine for recurring patterns in the feature space of a pre-trained CNN. We propose new objective functions enforcing the locality and unicity of the detected parts. Our detectors also provide a confidence measure based on correlation scores, allowing the system to estimate the visibility of each part. We apply our method on two public datasets, Caltech-UCSD Bird 200 (CUB-200) [START_REF] Welinder | [END_REF] and Stanford Cars [START_REF] Yang | A Large-Scale Car Dataset for Fine-Grained Categorization and Verification[END_REF], and show that our detectors can consistently highlight parts of the object while providing a good measure of the confidence in their prediction. Additionally, we provide classification results in order to showcase that classifiers based only on part detection can constitute a compromise between accuracy and transparency.

This paper is organized as follows: Section 2 presents the related work on part detection; Section 3 describes our PARTICUL model and Section 4 presents our results on two FGVC datasets. Finally, Section 5 concludes this paper and proposes several lines of research aiming at improving our approach.

Related work

Part detection (and more generally attribute detection) is a problem which has been extensively studied in recent years, especially for FGVC which is a notoriously hard computer vision task. Learning how to detect object parts in this context can be done either in a fully supervised (using ground-truth part locations), weakly supervised (using image labels only) or unsupervised manner.

Weakly supervised approaches [START_REF] Chen | This looks like That: Deep learning for interpretable image recognition[END_REF][START_REF] Ding | Selective Sparse Sampling for Fine-Grained Image Recognition[END_REF][START_REF] Ge | Weakly Supervised Complementary Parts Models for Fine-Grained Image Classification From the Bottom Up[END_REF][START_REF] Korsch | Classification-Specific Parts for Improving Fine-Grained Visual Categorization[END_REF][START_REF] Kun Duan | Discovering Localized Attributes for Fine-Grained Recognition[END_REF][START_REF] Li | Attribute Mix: Semantic Data Augmentation for Fine Grained Recognition[END_REF][START_REF] Liu | Fine-Grained Recognition with Automatic and Efficient Part Attention[END_REF][START_REF] Peng | Object-Part Attention Model for Fine-grained Image Classification[END_REF][START_REF] Simon | Neural Activation Constellations: Unsupervised Part Model Discovery with Convolutional Networks[END_REF][START_REF] Xiao | The Application of Twolevel Attention Models in Deep Convolutional Neural Network for Fine-grained Image Classification[END_REF][START_REF] Yang | Learning to Navigate for Fine-grained Classification[END_REF][START_REF] Zheng | Learning Multi-attention Convolutional Neural Network for Fine-Grained Image Recognition[END_REF][START_REF] Zhou | Human-Understandable Decision Making for Visual Recognition[END_REF] produce part detectors as by-products of image classification i.e., object parts are jointly learned with categories in an end-to-end manner to help distinguish very similar categories. In particular, the OPAM approach presented in [START_REF] Peng | Object-Part Attention Model for Fine-grained Image Classification[END_REF] obtains parts from candidate image patches that are generated by Selective Search [START_REF] Uijlings | Selective search for object recognition[END_REF], filtered through a dedicated pre-trained network, selected using an object-part spatial constraint model taking into account a coarse segmentation of the object at the image level, and finally semantically realigned by applying spectral clustering on their corresponding convolutional features.

Recently, unsupervised part detection methods have greatly benefited from the expressiveness and robustness of features extracted from images using deep CNNs. Modern approaches [START_REF] Krause | Fine-Grained Recognition Without Part Annotations[END_REF][START_REF] Sicre | Unsupervised Part Learning for Visual Recognition[END_REF][START_REF] Zhang | Unsupervised Part Mining for Finegrained Image Classification[END_REF][START_REF] Zhang | Picking Deep Filter Responses for Fine-Grained Image Recognition[END_REF][START_REF] Zheng | Learning Multi-attention Convolutional Neural Network for Fine-Grained Image Recognition[END_REF] mainly focus on applying clustering techniques in the feature space and/or identifying convolutional channels with consistent behaviors. More precisely, [START_REF] Krause | Fine-Grained Recognition Without Part Annotations[END_REF] produces part detectors through the sampling and clustering of keypoints within an estimation of the object segmentation produced by a method similar to GrabCut [START_REF] Rother | GrabCut" -Interactive Foreground Extraction using Iterated Graph Cuts[END_REF]. For a given image, [START_REF] Zhang | Unsupervised Part Mining for Finegrained Image Classification[END_REF] uses itemset mining to regroup convolutional channels based on their activation patterns, producing parts per image instead of globally (as in [START_REF] Peng | Object-Part Attention Model for Fine-grained Image Classification[END_REF]) and thus requiring an additional semantic realignment -through clustering -for part-based classification. Rather than working on the full images, [START_REF] Zhang | Picking Deep Filter Responses for Fine-Grained Image Recognition[END_REF] and [START_REF] Sicre | Unsupervised Part Learning for Visual Recognition[END_REF] use convolutional features of region proposals produced by Selective search [START_REF] Van De Sande | Segmentation as Selective Search for Object Recognition[END_REF]. [START_REF] Sicre | Unsupervised Part Learning for Visual Recognition[END_REF] then learns the relation between regions and parts using clustering and soft assignment algorithms, while [START_REF] Zhang | Picking Deep Filter Responses for Fine-Grained Image Recognition[END_REF] builds part detectors from the top-k most activated filters across all regions of the training set. Due to the independent training of each detector, the latter approach produces tens of redundant part detectors which must be filtered out in a weakly unsupervised manner. This issue is partially addressed in the MA-CNN approach [START_REF] Zheng | Learning Multi-attention Convolutional Neural Network for Fine-Grained Image Recognition[END_REF], where part detectors are learned by performing a soft assignment of the convolutional channels of a pretrained CNN into groups that consistently have a peak response in the same neighborhood, then regressing the weight (importance) of each channel in each group through a fully convolutional layer applied on the feature map. Each part detector produces an activation map, normalized using a sigmoid function, corresponding to the probability of presence of the part at each location. By picking the location with the highest activation value for each detector, MA-CNN generates part-level images patches that are fed into a Part-CNN [START_REF] Branson | Bird Species Categorization Using Pose Normalized Deep Convolutional Nets[END_REF] architecture for classification. In practice, part detectors are initialized by channel grouping (using k-means clustering) and pre-trained in an unsupervised manner using dedicated loss functions enforcing the locality and unicity of each detector attention region. Finally, these detectors are fine-tuned during end-to-end learning of the image category. As such, this method falls into the category of unsupervised part learning with weakly supervised fine-tuning. More recently, the P-CNN approach [START_REF] Han | P-cnn: Part-based convolutional neural networks for fine-grained visual categorization[END_REF] implements the part detection learning algorithm of MA-CNN, replacing the fully connected layer in charge of channel grouping by a convolution layer. Its classification pipeline uses Region-of-Interest (ROI) pooling layer around the location of maximum activation to directly extract part features and global features from the part detection backbone. In both cases, the channel grouping initialization requires first to process all images of the training set in order to cluster channels according to the location of their highest response. Moreover, for a given detector, only the area immediately surrounding the location of highest activation is taken into account (either to extract an image patch in [START_REF] Zheng | Learning Multi-attention Convolutional Neural Network for Fine-Grained Image Recognition[END_REF] or a vector through ROI-pooling in [START_REF] Han | P-cnn: Part-based convolutional neural networks for fine-grained visual categorization[END_REF]), thus reducing the ability to detect the same part at different scales.

Contributions Our PARTICUL approach builds part detectors as weighted sums of convolutional channels extracted from a pre-trained CNN backbone. Our detectors are trained globally across the training set, instead of locally, contrary to [START_REF] Zhang | Unsupervised Part Mining for Finegrained Image Classification[END_REF] and [START_REF] Peng | Object-Part Attention Model for Fine-grained Image Classification[END_REF] which require an additional semantic alignment phase. We develop objective functions that are specially crafted to ensure the compactness of the distribution of activation values for each part detector and the diversity of attention locations. Unlike [START_REF] Zheng | Learning Multi-attention Convolutional Neural Network for Fine-Grained Image Recognition[END_REF] and [START_REF] Han | P-cnn: Part-based convolutional neural networks for fine-grained visual categorization[END_REF], which propose similar functions, our approach does not require any fine-tuning of the backbone or an initial channel grouping phase. This induces a fast convergence during training and enables our module to be used with black-box backbones, in a plug-in fashion. Moreover, our detectors use a softmax normalization function, instead of sigmoid function, that simplifies the process of locating the part and enables us to detect parts at different scales. Finally, our detectors supply a measure of confidence in their decision based on the distribution of correlation scores across the training set. Importantly, this measure can also predict the visibility of a given part, a subject which is, to our knowledge, not tackled in any of the related work.

Proposed model

In this section, we present PARTICUL, our proposal for unsupervised part learning from FGVC datasets. In practice, each part detector can be seen as a function highlighting dedicated attention regions on the input image, in accordance with the constraints detailed in this section.

Notations For a tensor T inside a model M, we denote T (x) the value of T given an input x of M. If T has H rows, W columns and C channels, we denote T [h,w] (x) ∈ R C the vector (or single feature) located at the h th row, w th column. We denote * the convolution operation between tensors, σ the softmax normalization function. e define I as the set of all input images of a given dimension H I × W I and I C ⊆ I as the subset of images in I containing an object of the macro-category C. In our experiments, we denote

X train = x i ∈ I C i ∈ [1 . . n]
the training set of size n. Although I C cannot usually be formally specified, we assume that X train is representative of this set.

Part detectors

As illustrated in Fig. 1, we first use a pre-trained CNN F to extract a H ×W ×D feature map. Then, we build each of the p part detectors as a convolution with a 1 × 1 × D kernel k (i) (no bias is added), followed by a 2D spatial softmax layer, each part detector producing a H × W activation map

P (k (i) , x) = σ F (x) * k (i) , ∀x ∈ I (1)
We denote the set of all part detector kernels as K = [k (1) , k (2) , . . . , k (p) ].

Objective functions

Locality: For each part detector i, we force the network to learn a convolutional kernel k (i) which maximizes one region of the activation map for any training image x. In order to allow activations to be localized into a given neighborhood rather than in a single location in the H × W activation map, we first apply a 3 × 3 uniform filter u on P (k (i) , x). Let

G(k (i) , x) = P (k (i) , x) *   1 1 1 1 1 1 1 1 1   u (2)
Learning all convolutional kernels K in order to enforce locality can be translated as the optimization of the objective function

L l (K) = - 1 p p i=1 1 n x∈Xtrain max h,w G(k (i) , x) [h,w] (3)
Intuitively, solving the optimization problem in Eq. 3 using G(k (i) , x) rather than P (k (i) , x) relaxes the learning constraint and prevents the detectors from focusing on discriminative details between two adjacent feature vectors that would represent the same part.

Unicity In order to prevent the system from learning only a handful of easy parts, we wish to ensure that each feature vector F [h,w] (x) is not simultaneously correlated with multiple convolutional kernels k (i) . Let S(K, x) ∈ R H×W s.t.:

S [h,w] (K, x) = p i=1 P (k (i) , x) [h,w] (4)
Ensuring unicity can be translated as the following objective function, i.e., making sure that no location (h, w) contains a cumulative activation higher than 1:

L u (K) = 1 n x∈Xtrain max max h,w S [h,w] (K, x) -1, 0 (5)
The final objective function for the unsupervised learning of part detectors is a weighted composition of the functions described above:

L(K) = L l (K) + λL u (K) (6)
where λ controls the relative importance of each objective function.

Confidence measure and visibility

Visibility is related to a measure of confidence in the decision provided by each of our detectors and based on the distribution of correlation scores across the training set. After fitting our detectors to Eq. 6, we employ the function

H i : I → R x → max h,w F [h,w] (x) * k (i) (7)
returning the maximum correlation score of detector i for image x (before softmax normalization). The distribution of values taken by H i across I C is modeled as a random variable following a normal distribution N (µ i , σ 2 i ) estimated over X train . We define the confidence measure of part detector i on image x as

C(x, i) = Φ(H i (x), µ i , σ 2 i ) (8)
where Φ(z, µ, σ 2 ) is the cumulative distribution function of N (µ i , σ 2 i ).

Experiments

In order to showcase the effectiveness of our approach, we apply our algorithm on two public FGVC datasets -the Caltech-UCSD Birds 200 (CUB-200) [START_REF] Welinder | [END_REF] dataset containing 11,788 images from 200 bird species (5994 training images, 5794 test images) and the Stanford Cars [START_REF] Yang | A Large-Scale Car Dataset for Fine-Grained Categorization and Verification[END_REF] dataset containing 16,185 images from 196 car models (8144 training images, 8041 test images). In addition to the object subcategory labels (bird species, car model), both datasets also provide additional information in the form of annotations (object bounding box, part locations in [START_REF] Welinder | [END_REF]). However, when learning and calibrating our part detectors, and unless specified otherwise, we do not use any information other than the images themselves and work in a fully unsupervised setting.

Unsupervised learning of part detectors

In this section, we illustrate our unsupervised part detection algorithm using VGG19 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] (with batch normalization) pretrained on the Imagenet dataset [START_REF] Deng | ImageNet: A Large-Scale Hierarchical Image Database[END_REF] as our extractor F (.) (Eq. 1). For part visualization, we use SmoothGrad [START_REF] Smilkov | Smoothgrad: removing noise by adding noise[END_REF], filtering out small gradients using the method described in [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF]. We train our detectors during 30 epochs, using RMSprop with a learning rate of 5×10 -4 and a decay of 10 -5 . These training parameters are chosen using crossvalidation on the training set, with the goal of minimizing the objective function (Eq. 6). Since we do not need to fine-tune the extractor, only p×D convolutional weights are learned during training, which drastically reduces the computation time. Finally, for a given number of parts p, we supervise the importance λ of the unicity constraint (Eq. 6) by measuring the overall attention across the feature map after training. More precisely, we compute the value

E(X train , K) = 1 |X train | x∈Xtrain 1 p h,w max i∈[1. .p] P (k (i) , x) [h,w] (9)
corresponding to the average contribution of each detector. E(X train , K) = 1 corresponds to the ideal case where all detectors focus on different locations in the feature map F (x), while E(X train , K) = 1/p indicates a high redundancy of attention regions among detectors. As illustrated in Fig. 2, for both datasets detectors systematically converge towards the same location inside of the feature map (E(X train , K) ≈ 1/p) corresponding to the head of the bird; in the case of Stanford Cars, where images contain more varied distinctive patterns, the number of unique detectors depends on the random initialization of their convolutional weights, but the average contribution per detector quickly drops with p. For both datasets, we choose λ = 0.2 which maximizes E(X train , K) for all values of p. As a comparison, we also train our detectors after fine-tuning the extractor on each dataset (in this case, the learning process can be considered weakly supervised), leading to only marginally better results. This supports our claim that, in practice, our detectors do not require a fine-tuned extractor and can work as a plug-in module to a black-box extractor. As illustrated in Fig. 3, our detectors consistently highlight recurring parts of the objects and are relatively insensitive to the scale of the part. Although we notice some apparent redundancy of the detected parts (e.g., around the head of the bird), the relatively high corresponding value for E(X train , K) (0.87 for λ = 0.2 and p = 6, see Fig. 2a) indicates that each detector actually focus on a different location of the feature map, leading to a richer representation of the object. Note that when training our detectors, we do not know beforehand towards which part each detector will converge. However, using part visualization, we can re-attach a semantic value to each detector after training.

Confidence measure

After training our detectors, we perform a calibration of their confidence measure using the method presented in Sec. 3.3. In order to illustrate the soundness of our approach, we exploit the annotations provided by the CUB-200 dataset to extract a subset of images where the legs of the bird are non-visible (2080 images from both the training and test set, where the annotations indicate that the color of the legs is not visible). As shown on Fig. 4a, there is a clear difference in the distributions of maximum correlation scores between images with and without visible legs. This also confirms that images not containing the part tend to produce lower correlation scores, and by consequence our proposal for a confidence measure based on a cumulative distribution function (Eq. 8). In practice (Fig. 4b), a calibrated detector can detect the same part at different scales while ignoring images where the confidence measure is below a given threshold (2% in the example). It is interesting to note that for all images in the middle row, the annotations actually indicate that the legs are not visible, i.e., our detectors can also be employed as a fast tool to verify manual annotations.

Classification

In all related works closest to ours, the performance of part detectors is never evaluated independently but rather w.r.t. to the classification task, a method which is highly dependent on the architecture chosen for the decision-making process: e.g., SVMs in [START_REF] Krause | Fine-Grained Recognition Without Part Annotations[END_REF] and [START_REF] Sicre | Unsupervised Part Learning for Visual Recognition[END_REF], Part-CNN [START_REF] Branson | Bird Species Categorization Using Pose Normalized Deep Convolutional Nets[END_REF] in [START_REF] Zheng | Learning Multi-attention Convolutional Neural Network for Fine-Grained Image Recognition[END_REF], Spatially Weighted Fisher Vector CNN (SWFV) in [START_REF] Zhang | Picking Deep Filter Responses for Fine-Grained Image Recognition[END_REF], multi-stream architecture in [START_REF] Zhang | Unsupervised Part Mining for Finegrained Image Classification[END_REF]. Nevertheless, in order to provide a basis of comparison with state-of-art techniques, we also provide classification results based on the extraction of part feature vectors. As illustrated in Fig. 5, we use a method similar to [START_REF] Zhao | Recognizing Part Attributes With Insufficient Data[END_REF], where the feature Fig. 5: Our classification model for FGVC datasets. The feature map produced by the extractor is masked out using the activation map of each detector and processed through a pooling layer, resulting in a set of part vector features. These vectors are processed independently through a set of fully connected layers to generate part-based logits that are summed up to produce the final prediction. map F (X) is multiplied element-wise with the activation map produced by each detector to compute a set of part feature vectors that are used as inputs of p independent classifiers (each containing 2 intermediate layers -with 4096 neurons, ReLU activation and dropout -before the classification layer). As such, our detectors operate a form of semantic realignment and extraction of relevant feature vectors from the image. Note that, contrary to the fixed size ROI pooling used in [START_REF] Han | P-cnn: Part-based convolutional neural networks for fine-grained visual categorization[END_REF], this method has the advantage of dynamically adjusting the number of feature vectors taken into account depending on the scale of the part. For the final decision, we sum up the logits of all part-based classifiers to produce prediction scores for each category. Therefore, the final decision can be directly traced back to the individual result of each part-based classifier. This approach does not provide the transparency [START_REF] Lipton | The mythos of model interpretability[END_REF] of prototype-based methods [START_REF] Chen | This looks like That: Deep learning for interpretable image recognition[END_REF][START_REF] Nauta | Neural Prototype Trees for Interpretable Finegrained Image Recognition[END_REF], but constitutes a good compromise to non interpretable approaches using global features. 1 summarizes the results obtained when using our part-based classifier, along with the results obtained on a fine-tuned VGG-19 (acting as a baseline) and other state-of-the-art methods. With 8 detectors, we outperform the prototype-based approach [START_REF] Nauta | Neural Prototype Trees for Interpretable Finegrained Image Recognition[END_REF] -which uses a more efficient extractor (Resnet50) -on Stanford Cars and obtain similar results on CUB-200. When compared with other methods using only part-level features, our PARTICUL model outperforms the OPAM [START_REF] Peng | Object-Part Attention Model for Fine-grained Image Classification[END_REF] approach on both datasets. We also obtain comparable results to other methods requiring either the object bounding box [START_REF] Krause | Fine-Grained Recognition Without Part Annotations[END_REF] or to preselect detectors (in a weakly supervised manner) based on their classification accuracy [START_REF] Zhang | Picking Deep Filter Responses for Fine-Grained Image Recognition[END_REF]. When compared with less transparent methods using image-level (global) features in addition to part-level features -a method which usually has a significant impact on accuracy (e.g., from 80.7% to 85.8% on CUB-200 using OPAM [START_REF] Peng | Object-Part Attention Model for Fine-grained Image Classification[END_REF]) -again we achieve comparable results and even outperform the UPM [START_REF] Zhang | Unsupervised Part Mining for Finegrained Image Classification[END_REF] approach on CUB-200. Finally, it is also interesting to note that for both datasets, our classification results obtained by picking 10 feature vectors out of the 14×14 = 196 possible vectors of the extractor feature map correspond to a drop of less then 1% in accuracy when compared with the baseline, indicating that we are indeed selecting the most relevant vectors for the classification. Moreover, in 80% (resp. 77%) of these cases where only the baseline provides a correct prediction on the CUB-200 (resp. Stanford Cars) dataset, at least one of our individual part-based classifier does provide a correct prediction (see Fig. 6). Thus, our proposed model could be further improved by fine-tuning the relative importance of part logits (e.g., using the confidence measure).

Conclusion and future work

In this paper, we presented our algorithm for unsupervised part learning using datasets for FGVC. We showed that our detectors can consistently highlight parts of an object while providing a confidence measure associated with the detection. To our knowledge, our method is the first to take the visibility of parts into account, paving the road for a solid attribute learning and ultimately for interpretable visual recognition. In the particular context of FGVC, our detectors can be integrated in a part-based classification architecture which constitutes a good compromise between the transparency of prototype-based approaches and the performance of non-interpretable methods. As a future work, we will study the integration of our detectors into a prototype-based architecture, learning prototypes from part feature vectors rather than from the entire image feature map. We will also study the impact of weighting part logits by the confidence score associated with the detected part on the overall accuracy of the system.

Fig. 1 :

 1 Fig. 1: Architecture of our PARTICUL model. (1) Convolutional features F (x) are extracted from the image x. (2) Each part detector produces an activation map P (k (i) , x). (3) We apply a uniform kernel u to each part activation map before computing the Locality loss which ensures the compactness of activations. (4) All part activation maps are summed in S(K, x). (5) Unicity loss is applied to ensure the diversity of part detectors. Here, part detectors 2 and p are very similar, leading to a high peak in bright red in S(K, x). Best viewed in color.

Fig. 2 :

 2 Fig. 2: E(X train , K) v. number p of detectors for various values of λ, with pretrained or fine-tuned extractor.

Fig. 3 :

 3 Fig. 3: Part visualization after training 6 detectors using the VGG19 extractor on Stanford cars (top) or CUB-200 (bottom). Using this visualizations, we can manually re-attach a semantic value to each detector, e.g., the second part detector trained on CUB-200 is probably a "leg" detector. Best viewed in color.

  (a) Distribution of maximum correlation scores on the CUB-200 training set (in blue) and on a subset containing only images with non-visible legs (red). (b) Confidence scores and part visualizations on images with non visible legs (top-row) and with visible legs (bottom rows).

Fig. 4 :

 4 Fig. 4: Confidence measure applied on a bird leg detector trained and calibrated on CUB-200 dataset. Best viewed in color.

Fig. 6 :

 6 Fig. 6: Examples of images correctly classified by the baseline and incorrectly classified by our model (only 5 parts are shown for clarity). In both cases, at least one part-based classifier provides a correct prediction. Best viewed in color.

Table 1 :

 1 Classification accuracy on CUB-200 and Stanford Cars and comparison with related works, from less transparent to most transparent. Table

		Method	Train anno. CUB-200 Accuracy (%) Stanford Cars
		Global features only
		Baseline (VGG-19)	83.3	89.3
		Part + global features
		UPM [39] (VGG19)	81.9	89.2
	Transparency [20] ←-	OPAM [25] (VGG-16) MA-CNN [42] (VGG-19) P-CNN [9] (VGG-19) Part-based Ours (VGG-19) 2 parts 4 parts 6 parts 8 parts	85.8 86.5 87.3 70.1 79.2 81.5 82.3	92.2 92.8 93.3 76.0 84.2 87.5 88.3
	←-	10 parts OPAM [25] (parts-only, VGG-16)	82.3 80.7	88.6 84.3
		No parts [15] (VGG19)	BBox 82.0	92.6
		+ Test Bbox	BBox 82.8	92.8
		PDFS [40] (VGG-19)	84.5	n/a
		Prototypes
		ProtoPNet [3] (VGG-19)	BBox 78.0 ± 0.2 85.9 ± 0.2
		ProtoTree [22] (Resnet-50)	82.2 ± 0.7 86.6 ± 0.2

Acknowledgements Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations. This work has been partially supported by MIAI@Grenoble Alpes, (ANR-19-P3IA-0003) and TAILOR, a project funded by EU Horizon 2020 research and innovation programme under GA No 952215. This work was financially supported by European commission through the SAFAIR subproject of the project SPARTA which has received funding from the European Union's Horizon 2020 research and innovation programme under GA No 830892, as well as through the CPS4EU project that has received funding from the ECSEL Joint Undertaking (JU) under GA No 826276.


J Bouchet

F Bottin

D Antonangeli

G Morard

properties of hcp iron at high pressure and temperature

Sound velocities and thermodynamical properties of hcp iron have been computed using ab initio calculations over an extended density and temperature range, encompassing the conditions directly relevant for the Earth's inner core. At room temperature, and up to 350 GPa, an excellent agreement is obtained between present results and experimental data for many thermodynamical quantities: phonon density of states, vibrational entropy, heat capacity, Grüneisen parameter and thermal expansion. With increasing temperature, along an isochore, we observe a strong decrease of the phonon frequencies, demonstrating that intrinsic anharmonic effects cannot be neglected. We also carefully compare previous theoretical data for the sound velocities and try to explain the discrepancies observed with experiments. Finally, we propose a temperature dependant Birch's law that we compare with previous experimental work.

Introduction

Due to its technological and geophysical importance, iron is one of the most studied element of the periodic table. A large experimental effort has been devoted to the measure of its thermodynamical properties at pressures and temperatures relevant to Earth's interior. Among these properties the sound velocities [1,2,3,4,5] are key quantities since they can be directly related to seismological observations [6,7]. Probing matter at such thermodynamical conditions (∼ 330 GPa, ∼ 6000 K) is extremely challenging and in fact constraints are obtained by extrapolationof data at lower pressures and temperatures. In parallel to these experimental works and to circumvent the experimental limits, ab initio calculations have been performed to obtain sound velocities [START_REF] Price | geophysical and Geochemical Evolution of the Deep Earth[END_REF][START_REF][END_REF]10,[START_REF] Tsuchiya | Physics of the Earth and Planetary Interiors[END_REF][START_REF] Sha | [END_REF]. In addition, other thermodynamical properties as thermal expansion or Grüneisen parameter have been calculated [13,14,15] to build reliable equation of state (EOS) and compare to experimental data [16,5,17].

Calculations at high pressure with density functional theory (DFT) can be performed simply by changing the value of the volume. Including temperature is a longstanding task, much more complex. Many approximations have been made, mostly based on the harmonic approximation, to indirectly simulate the temperature and describe the thermal behaviour of physical quantities. One of the best known is the quasi-harmonic approximation (QHA) [18] which introduces an implicit temperature dependence of the phonon frequencies and leads to a non-zero thermal expansion. The QHA is easy to implement and only requires 0 K calculations [14]. However, its temperature range of validity is difficult to assess and can strongly depend on the materials under consideration, since it is now well established that anharmonic effects, beyond the QHA, can be important far below the melting temperature [19,20,21].

In the past ten years, strong efforts have been made to take into account explicit temperature effects and go beyond the QHA. New methods capturing the thermal properties of solids at non-zero temperature are now available and can be applied in ab initio calculations. These approaches combine ideas including finite large displacements, molecular dynamics sampling, self consistent harmonic theories, and different force fitting schemes. The most widely used methods include Self-Consistent Ab Initio Lattice Dynamics (SCAILD) [22], Stochastic Self-Consistent Harmonic Approximation (SSCHA) [23,24], Temperature Dependent Effective Potential (TDEP) [25,26,27], Anharmonic LAttice MODEl (ALAMODE) [28] and Compressive Sensing Lattice Dynamics [29]. Other methods obtain anharmonic contributions via thermodynamic integration [30], or a series expansion of the interatomic forces constants [31]. A large number of new phenomena, intrinsically temperature dependent, can now be captured: the modification of the phonon density of states (vDOS) and free energy, the (T , P ) phase transition boundaries [32], the evolution of elastic constants [33] or Grüneisen coefficients, the phonon lifetimes [34], mechanisms behind superconductivity in pressurestabilized hydrides [35] or the thermal conductivity [36].

In this paper we perform ab initio molecular dynamics (AIMD) simulations and A c c e p t e d M a n u s c r i p t use the TDEP method to calculate the vibrational properties of iron up to pressure and temperature relevant for the Earth's core. In the next two sections we give the details of the calculations and the different equations used to calculate various physical quantities. Then we present results for the thermodynamical properties in section 4 and for the sound velocities in section 5 with careful comparison with previous theoretical and experimental data.

Computational Details

Ab initio simulations were performed using the ABINIT package [START_REF]The ABINIT code is a common project of the Catholic University of Louvain (Belgium), Corning Incorporated, CEA (France) and other collaborators[END_REF][START_REF] Gonze | [END_REF]39] in the framework of the Projector Augmented Wave (PAW) method [40,41]. We employ the generalized gradient approximation (GGA) according to the parametrization of Perdew, Burke and Ernzerhof (PBE) for the exchange-correlation energy and potential [42].

Using ATOMPAW [43,[START_REF]ATOMPAW is a general license public code developed at Wake Forest University. Some of its capabilities have been developed at the Commissariat à l'énergie atomique[END_REF][START_REF] Dewaele | [END_REF], we generate a small core PAW atomic data with a radius r PAW equal to 2.0 u.a, with 3s, 3p, 3d, and 4s states as valence electrons and a cutoff energy equal to 350 eV for the plane wave basis set. The simulation box includes 180 atoms and corresponds to a 5 × 3 × 3 supercell of a 4-atoms C-centered orthorhombic cell (hcp structure with orthogonal axes). A good integration of electronic quantities is fulfilled using a 2 × 2 × 2 Monkhorst-Pack mesh leading to the inclusion of 4 special k-points in the calculations. The size of the supercell and the k-points mesh used in present calculations give an uncertainty of less than 1% on the velocity, which cannot be obtained using smaller supercells (128 and 144 atoms) whatever the k-points mesh used (Γ point or 3 × 3 × 3).

Tens of AIMD simulations are performed in the N V T ensemble (constant number of particles, constant volume and temperature) and are run for about 3-5 ps using a time step (τ ) of 1.13 fs. In particular, AIMD simulations are carried out for five densities ρ equal to 10, 10.8, 12, 13.15 and 14 g.cm -3 corresponding to pressures between 50 and 350 GPa, and several temperatures from 300 up to 7000 K depending on the density considered. The brodening of the occupation numbers was beased on the Fermi-Dirac statistics with the electronic temperature equal to the ionic one. Taking benefit of an efficient scheme of parallelization [46] and using hundreds to thousands of processors, the recovery time is a few months.

Calculations of thermodynamic properties and sound velocities

To obtain dynamic, elastic and thermodynamic properties as a function of the temperature we used the TDEP approach developed by Hellman and coworkers [25,26] and implemented in ABINIT [47,27] as aTDEP. In this method, a model Hamiltonian expanded as a function of the atomic displacements around equilibrium and truncated at the second or third order is adjusted to fit the potential energy surface obtained using molecular dynamic simulations at finite temperature. As a results, we obtain the effective interatomic force constants (IFC) at the second Φ αβ ij (T ) and third Ψ αβγ ijk (T )

A c c e p t e d M a n u s c r i p t

Sound velocities and thermodynamical properties of hcp iron at high pressure and temperature4 order, between i, j, k atoms and for any α, β, γ directions. These effective IFC are no longer constant as a function of the temperature but now depend on it. By performing a Fourier transform of the effective 2 nd order IFC, we can get the dynamical matrix, the eigenvectors X α is (q, T ) and eigenvalues (phonon frequencies) ω s (q, T ) at any q-point of the Brillouin zone and for each phonon mode s. In addition, aTDEP provides the vDOS g(ω, T ) and the associated thermodynamic properties: the specific heat C V (T ), the vibrational entropy S vib (T ), the vibrational energy U vib (T ) and the vibrational free energy F vib (T ).

Another thermodynamic quantity which can be derived from the IFCs is the Grüneisen parameter

γ = V ∂P ∂U V
, where U is the internal energy. For iron at Earth's inner core conditions, its determination is essential because it provides the link between the thermal pressure P vib and the temperature (at high temperature, in the classical limit, U vib (T ) = 3N k B T , with k B the Boltzmann constant). Consequently, this parameter is a key quantity for the Earth's interior where P vib as a function of depth is well constrained, thanks to seismological data, but T is poorly known. The evolution of γ(ρ, T ) with density ρ, but also with temperature is therefore essential to build reliable equation of states (EOS) for planetary modeling [48,49].

The mode Grüneisen parameters are defined by the volume derivatives of the phonon-mode frequencies:

γ s (q, T ) = - V ω s (q, T ) ∂ω s (q, T ) ∂V T (1)
These later can be obtained [27] from the effective 3 rd order IFC Ψ αβγ ijk (T ) as follows:

γ s (q, T ) = - 1 6ω 2 s (q, T ) ijkαβγ Ψ αβγ ijk (T ) × X α is (q, T )X β js (q, T ) M i M j r γ k exp [iq.R j ] (2)
where M i is the mass, r α i the vector position of atom i along α and R i the lattice vector of the unitcell of atom i. From the mode Grüneisen parameters γ s (q, T ), the thermodynamical Grüneisen parameter writes:

γ(T ) = 3Na s=1 q∈BZ γ s (q, T )C V,s (q, T ) C V (T ) (3)
where the C V,s (q, T ) are the mode heat capacities such as C V (T ) = 3Na s=1 q∈BZ C V,s (q, T ) and BZ is the Brillouin zone.

As previously performed for the Grüneisen parameter, we can also obtain the temperature dependence of the isothermal elastic constants by using the effective 2 nd IFC and the formulation proposed in Refs. [START_REF] Leibfried | Theory of anharmonic effects in crystals[END_REF][START_REF] Wallace | Thermodynamics of Crystals Dover books on physics[END_REF]33]:

C αβγδ (T ) = A αγβδ (T ) + A βγαδ (T ) -A αβγδ (T ) (4)
A c c e p t e d M a n u s c r i p t

Sound velocities and thermodynamical properties of hcp iron at high pressure and temperature5

with,

A αβγδ (T ) = 1 2V ij Φ αβ ij (T )d γ ij d δ ij (5)
where d α ij is the distance between the atom i and j along α. The full isothermal elastic tensor C αβγδ (T ) can be reduced following the symmetries and the Voigt notation. The temperature effects are then directly introduced in the elastic constants by the temperature variation of the IFCs.

Once the C αβ (T ) are obtained, the isothermal incompressibility K T (T ) and the shear modulus G(T ) are given using the Voigt average for the hcp structure by:

K T (T ) = 2(C 11 + C 12 ) + 4C 13 + C 33 9 (6)
G(T ) = 7C 11 -5C 12 + 2C 33 -4C 13 + 12C 44 30 (7)
We have also calculated K T and G using the Reuss and Hill formulas, the differences were less than 0.3% compared to the Voigt average. From the Grüneisen parameter, the isothermal compressibility β T (T ) = 1 K T (T ) and the heat capacity, it is also possible to compute the thermal expansion of the system, with all the intrinsic and extrinsic thermal effects:

α(T ) = γ(T )C V (T ) β(T )V (8)
Finally, the adiabatic incompressibility K S (T ), the compressional V P (T ) and shear V S (T ) velocities are obtained using all the previous results coming from the 2 nd and 3 rd orders, as follows:

K S (T ) = K T (T ) 1 + α(T )γ(T )T (9)
V P (T ) = K S (T ) + 4/3G(T ) ρ (10) V S (T ) = G(T ) ρ (11)
This procedure enables us to take into account the anharmonic effects not included in the QHA where only implicit effects coming from the thermal dilation is considered. This missing explicit anharmonic part, also called intrinsic, can be crucial for finite temperature elasticity, as shown on TiAlN alloys [33].

For comparison, the procedure to obtain the C αβ with the QHA [START_REF] Steinle-Neumann | [END_REF]14] is very tedious since the phonon spectra has to be calculated not only for the parent structure but also for all the strained structures with different values for the applied strain. So between 20 to 30 first principles linear-response calculations are necessary at each volume. Note that if the internal degree of freedom are relaxed for the strain structures at 0 K, their temperature dependency is not taken into account with this procedure.

A c c e p t e d M a n u s c r i p t

Similar criticisms can be made for the method where the strains are applied to supercells [START_REF][END_REF]53,10]. At least 10 AIMD simulations are necessary with running time around 10 ps [10], resulting in an overall simulation time of 100 ps to obtain the set of elastic constants of the hcp structure at one thermodynamic point. By comparison, using TDEP, we obtain the elastic constants with a single simulation of about 5 ps. The main drawback of this method is that since the elastic constants are related to phonon frequencies in the long wavelenght limit, it requires large supercells to converge the sum in Eq. ( 5). To circumvent this difficulty, we only retain the temperature dependence of Eq. ( 4) that we apply to the static calculations at 0 K [33,32].

Thermodynamic properties of -Fe

Phonon spectrum, phonon density of state and related quantities

As a first illustration, we compare in Figure 1 the vDOS obtained in present AIMD calculations at room temperature to the results obtained in nuclear resonant inelastic X-ray scattering (NRIXS) measurements on isotopically enriched 57 Fe [54](we consider natural 56 Fe in present calculations). The excellent agreement obtained attests the validity of vDOS computed as a function of the pressure (between 50 and 150 GPa) at low temperature and assess the quality of present calculations at higher pressure and temperature where no experimental data are available.

The phonon dispersion relations and the associated vDOS of hcp Fe at 13.06 g.cm -3 as a function of temperature are presented in Fig. 2. First, we do not observe any dynamical instabilities for the hcp structure even at 7000 K, a temperature close to the melting point (ab initio calculations evaluate the melting temperature around 6300 K at this density [56,57] but the solid state can be maintained at higher temperature due to the overheating effect [15]). Secondly, we observe that the frequencies soften with temperature but not with the same amplitude. Up to 5000 K the temperature effect is rather small for the longitudinal branches (acoustic and optic) while we observe a stronger effect on the transversal branches. Between 5000 and 7000 K, the longitudinal branches decreases considerably. This is also clearly seen on the vDOS where the peak at high frequencies is abruptly shifted between 5000 and 7000 K, while the peak at low frequencies softens monotonously with temperature. The initial slope at low energy, from which the Debye sound velocity can be derived, strongly increases with temperature, an effect already noticed experimentally [1] but at lower temperature (from 300 to 1100 K).

As already noticed for bcc Fe [58] and for hcp Fe [13], Figure 2 demonstrates the importance of the explicit anharmonic effects in iron at high temperature. All the temperature effects highlighted in Fig. 2 cannot be reproduced using QHA with fixed phonon spectrum at constant volume.

From the vDOS computed at each thermodynamic point we extract the vibrational entropy S vib and the vibrational heat capacity C V , and compare them to experimental A c c e p t e d M a n u s c r i p t values (see Fig. 3). Once again the agreement is really encouraging and proves the reliability of the present method. Our results for the vibrational entropy at 300 K are also in good agreement with the calculations of J. Zhuang et al. [59]. As shown in this work, above 1000 K, the electronic contribution to the specific heat is no longer negligible and it starts to deviate from the Dulong-Petit law.

Grüneisen parameter

The Grüneisen parameter is an important quantity for thermal EOS and to extrapolate thermophysical properties to high pressures and temperatures. It can be deduced from several measured quantities, as the vDOS [55], the atomic mean square displacements [61], Raman spectroscopy [62] or by comparing the Hugoniot and measured isentrope, as done recently at the National Ignition Facility (NIF), up to 1.4 TPa [63].

Concerning the building of equation of state (EOS), a commonly used expression for the volume dependence of the Grüneisen parameter is:

γ = γ 0 V V 0 q ( 12
)
where γ 0 and V 0 are the Grüneisen parameter and the volume at ambient conditions, Table 1. Parameters of Eq.( 12) for the variation of the Grüneisen parameter.

Reference γ 0 q

Dubrovinsky [61] 1.78 0.69 Murphy [55] 1.88 0.8 Merkel [62] 1.68 0.7 Miozzi [17] 1.11 0.3 Present study 1.73 0.68 and q is a fitted parameter. The result of the fitting of present data with Eq. ( 12) is given in Table 1 in comparison with previous studies.

Present raw data at 300 K and the corresponding fitting curve are reported in Fig. 4 and compared to various experimental data and fits. Present room temperature values are in close agreement with the experimental points, specially those obtained by X-ray diffraction [61]. At high pressure, the fit also confirms the experimental constraint on the Grüneisen parameter obtained at the NIF.

The behaviour in temperature of the Grüneisen parameter is shown in Fig. 5. At the highest densities here considered (13 to 14 g.cm -3 ), it slightly increases at low temperature then shows a strong decrease due to anharmonic effects, following the decrease observed on the phonon frequencies (see Fig. 2). At high temperatures and densities, the Grüneisen parameter is almost constant, as already observed by Alfè et al. [13]. At 6000 K, they obtain a constant value of 1.45, while at 7000 K we have around 1.0. Since here we only calculate the vibrational contribution, we can estimate the electronic part around 0.3-0.4 at temperatures and densities relevant for the Earth's core.

Thermal expansion

As mentioned in section 3, the thermal expansion matrix can be calculated using the third order effective IFCs. For an hcp crystal, there is two coefficients of thermal expansion (CTE), one for the basal plane (α a ) and a second one for the perpendicular direction along the c axis (α c ), while the volumetric CTE is given by α V = 2α a + α c . Present results at 300 K are presented in Fig. 6 and compared to the volumetric CTE deduced from the experimental vDOS obtained with NRIXS [60]. As for the vDOS and the vibrational entropy, the agreement is excellent. We observe an anisotropy between the thermal expansion along the a and c axis, α c being larger than α a . Therefore the c/a ratio increases with temperature as found in previous theoretical [65,66] and experimental studies [67,68]. This anisotropy decreases at core density where α a and α c are close to each other at 300 K [69]. The black diamonds are the results of present calculations. The green circles, the red triangles and the blue squares are experimental data from Ref. [55,64,63] respectively. The black, orange, violet and green lines are fits following Eq. ( 12) for present data and Ref. [61,62,55].

A c c e p t e d M a n u s c r i p t

The effect of temperature on the volutmetric CTE is presented in Fig. 7. It increases between 300 and 1000 K due to the filling of the phonon energy levels, then decreases at higher temperatures due to anharmonic effects (a trend similar to the one observed for the Grüneisen parameter, according to Eq. ( 8)). At 5000 K, we slightly underestimate the experimental value of Duffy and Ahrens [69], possibly because of the electronic contribution. In present calculations we also observe that the anisotropy between the basal plane and the c axis decreases with temperature so the c/a ratio stays below the ideal value of 1.633 in agreement with Gannarelli et al. [66].

Elastic properties and sound velocities

Comparison with previous ab initio calculations at 0 K

In Fig. 8, we compare present compressional and shear velocities for hcp-Fe at 0 K to previous ab initio calculations with different methods: PAW [START_REF] Price | geophysical and Geochemical Evolution of the Deep Earth[END_REF]53,10], ultrasoft pseudopotentials [START_REF] Tsuchiya | Physics of the Earth and Planetary Interiors[END_REF] or all-electrons methods [START_REF] Steinle-Neumann | [END_REF][START_REF] Sha | [END_REF]70]. All these calculations use the GGA functional to describe the exchange-correlation energy and potential. As expected, the agreement between all these results is really good. The discrepancies observed with the PAW calculations of Vocadlo et al. [START_REF] Price | geophysical and Geochemical Evolution of the Deep Earth[END_REF] (filled red squares in Fig. 8) are certainly due to a less dense k-points grid or a smaller number of valence electrons. Their more recent results [10,53] (open red squares in Fig. 8) are in better agreement with other calculations. We also observe large discrepancies with the data of Sha and Cohen [START_REF] Sha | [END_REF] (filled orange down triangles in Fig. 8) obtained using full-potential linearmuffin-tin-orbital (FPLMTO) calculations. This is unexpected and seems incoherent with calculations using full charge density exact muffin-tin orbital (FCD-EMTO) method [71,70]. The difference increasing with density, it is possible that a too large muffin-tin radius was used in the work of Sha and Cohen, resulting in an overlap between the spheres around each ion and a smaller pressure than expected. This would explain why their compression curve strongly diverges from the experimental one at low volumes (see the figure 4.e of Ref. [START_REF] Sakai | Physics of the Earth and Planetary Interiors[END_REF]). If we put aside these two works [START_REF] Price | geophysical and Geochemical Evolution of the Deep Earth[END_REF][START_REF] Sha | [END_REF], we observe a linear dependence on density for the sound velocities, with differences below 2% between the calculations, confirming the validity of the Birch's law for athermal calculations (i.e. the validity of the quasi-harmonic approximation).

Comparison with experimental data at 300 K

We compare present theoretical velocities to the experimental data on Fig. 9. Two main techniques are used to obtain the velocities: the NRIXS and the inelastic X-Ray Theoretical compressional V P and shear V S velocities for hcp-Fe. Black circles, PAW, this study; red squares, PAW, Ref. [START_REF] Price | geophysical and Geochemical Evolution of the Deep Earth[END_REF]; open red squares, PAW, Ref. [53,10]; green triangles, ultrasoft pseudopotentials, Ref. [START_REF] Tsuchiya | Physics of the Earth and Planetary Interiors[END_REF]; violet diamonds, FPLAW, Ref. [START_REF] Steinle-Neumann | [END_REF]; orange down triangles, FPLMTO, Ref. [START_REF] Sha | [END_REF]; brown left triangles, FCD-EMTO Ref. [70]; blue right triangles, FPLMTO, Ref. [71].

scattering (IXS), with different drawbacks and advantages. Basically NRIXS is more reliable to provide V S than V P while for IXS it is the opposite. For a complete discussion and a comparison of these techniques, see Ref. [2], which also provides a fit to combined data sets available in literature, aslo including results by impulsive stimulated light scattering [START_REF] Crowhurst | [END_REF] and picosecond acoustics [74].

Clearly, the theoretical results overestimate the experimental velocities for the low densities, but the difference diminishes at higher densities. For V P , at 10 g.cm -3 , we have a difference of about 8 % with the fit of Ref. [2], which drops to 2 % at 13 g.cm -3 if we extrapolate the experimental data. The comparison is worse with the sub-linear relation proposed by Mao et al. [75] at core densities with differences up to 8%. The difference with the NRIXS measurements is larger, about 12 % from 10 to 12.5 g.cm -3 . For V S , the situation is even worse, with differences going from 18 to 12% respectively at 10 and 13 g.cm -3 .

Can we understand these discrepancies and the reasons behind, and more, can we propose a way to correct them? To answer these questions, we need to turn to the comparison between the experimental and ab initio EOS of iron. Fortunately, the EOS obtained with the PAW atomic datasets that we use here have been extensively compared with the more recent EOS derived from diamond anvil cell data [77,[START_REF] Dewaele | [END_REF][START_REF] Sakai | Physics of the Earth and Planetary Interiors[END_REF]16],

A c c e p t e d M a n u s c r i p t Figure 9. Comparison with experimental compressional V P and shear V S velocities for hcp-Fe. Theory: black circles, PAW, this study (same as in Fig. 8); open circles, PAW with corrected pressure. Experiments: red line, linear fit of combined datasets, Ref. [2]; green line, linear fit of IXS data, Ref. [3]; orange line, power law fit of IXS, Ref [75]; red circles, blue diamonds, violet triangles and brown down triangles, NRIXS measurements of Ref. [76], [60], [54], and [1] respectively. PREM is shown as crosses. For the fits, the lines cover the data points while the dashed ones are extrapolations.

either using Vinet or Birch-Murnaghan EOS. These comparisons bring out a number of important points. Below 200 GPa, the ab initio EOS underestimates the experimental pressure, while above the pressure is overestimated. Dewaele et al. [START_REF] Dewaele | [END_REF] also emphasize an overestimation of the bulk modulus with the filling of the d-orbitals due to the GGA approximation itself. In conclusion they propose a procedure to correct the ab initio EOS (i) by using the experimental equilibrium volume instead of the theoretical one, and (ii) by using the K T and K T recalculated at this volume.

In the same spirit we can try to correct present theoretical sound velocities. First, using the theoretical and experimental parameters of the Vinet EOS, we have corrected present densities to obtain the experimental pressure (see the open circles in Fig. 9). As direct consequence of this correction, velocities are decreased for densities lower than 12.5 g.cm -3 (the density where both EOS give the same pressure) and increased above. The effect is significant for V P and improve the agreement with experiments but is almost negligible for V S due to a weaker slope. Interestingly this correction gives a density dependence of the compressional velocity in close agreement with the fit proposed by Antonageli et al [2] and Sakamaki et al [3] with a constant overestimation A c c e p t e d M a n u s c r i p t Sound velocities and thermodynamical properties of hcp iron at high pressure and temperature15 of the compressional velocities around 4% between 10 and 14 g.cm -3 .

Once the pressure correction performed, the discrepancies that we still observe for the shear velocity can only be explained by an overestimation of the shear modulus G in the calculations. If we compare the K T obtained with present athermal calculations and the values taken from the EOS of Sakai et al. the difference is around 9% for the range of densities considered here. An increase of G around 20 % would be needed to find a good agreement with experiments for both velocities.

The second correction that can be made to the theoretical results comes from the temperature difference between calculations and experiments, (0 K vs. 300 K). In the results presented in the next section we have included the velocities that we obtain at room temperature. This correction lowers the velocities, with a larger effect for V S , but is still not sufficient to explain the discrepancies. 8) and (3) (see sections 4.3 and 4.2). Both velocities decrease with temperature but V S much more rapidly than V P . This is simply due to the fact that C 11 and C 33 decrease with temperature while C 12 and C 13 increase. So, according to equation (6), these effects compensate each other for K T while they amplify each other for G. The temperature effects on velocities weaken as the denity increase.

We fit these data using the formula proposed by Sakamaki et al. [3]:

V (ρ, T ) = M ρ + B + A(ρ -ρ * )T (13)
where M and B are the coefficient of Birch's law at 0 K while A and ρ * give the temperature dependence of the velocities. In equation ( 13), T is in K, density in g.cm -3 and V in kms -1 . In the present fit, we discard the data at ρ = 10.22 g.cm -3 and 3000 K since at these thermodynamic conditions iron is either liquid or fcc [78]. For V P we found M = 1.20, B = -3.66, A = 4.55 × 10 -5 and ρ * = 15.68 g.cm -3 . For V S we found M = 0.50, B = -0.36, A = 5.24 × 10 -5 ρ * = 19.6 g.cm -3 . We plot in Fig. 11 a comparison between the results of the present fit and the direct calculation of the velocities. The differences do not exceed 2 % and are larger for V S than for V P . We therefore conclude that for the range of densities considered here, and up to 7000 K, velocities along an isotherm linearly depend on density and that it is not necessary to introduce a power law function for V P as proposed in Ref. [75]. As expected the temperature effects are more important at low densities and tends to reduce at the inner core boundary (ICB) but since the melting temperature is much higher they cannot be neglected.

5.3.2.

Comparison with previous theoretical work At 12.52 g.cm -3 , when the temperature increases, Sha and Cohen [START_REF] Sha | [END_REF] obtain a reduction of the compressional A c c e p t e d M a n u s c r i p t velocity equal to 0.13 m.s -1 .K -1 , which is comparable to the present value of 0.15 m.s -1 .K -1 . The difference in the case of V S is larger, they found 0.35 m.s -1 .K -1 to be compared with present value of 0.39 m.s -1 K -1 . At 10.43 g.cm -3 , they obtain a value around 0.09 m.s -1 .K -1 for V P , much lower than present result 0.27 m.s -1 .K -1 . This difference is certainly due to anharmonic effects beyond the thermal expansion and not taken into account by the quasi harmonic approximation. On an isotherm, these anharmonic terms are supposed to be more important at low densities and this would explain these discrepancies.

We use the (ρ, T ) values of Refs. [53,10] in the present fit (see Eq. ( 13)) to compare their sound velocities values with present results (see Fig. 12). The agreement is really good, all the more so when considering that the correction on present densities (coming from the discrepancy between the theoretical and experimental EOS) is not taken into account in these previous studies.

However we observe larger differences at high temperature. Elastic properties are difficult to obtain at temperatures close to the melting point. For example at 13 g.cm -3 and 7200 K, we observe a distorsion in the hcp supercell. A plane of atoms slides along the [010] direction and the hcp structure is no longer stable. This is expected to have a strong impact on the values of the elastic constants which can considerably 13)). The colors of the temperatures are the same as in Figure 10.

drop. Unfortunately these phenomena are related to the size of the supercell, and the conclusions from those simulations are limited. It could be premelting effect, but also a transition towards another solid structure as bcc which could be more stable than the hcp structure at temperatures close to the melting point [79]. Martorell et al. [10] observed the same phenomena but on smaller supercells (64 atoms compared to 180 atoms in present simulations). A deformation of the supercell as used in this work to calculate the elastic constants could also induce sliding planes. We believe that these results should be treated with caution.

Comparison with experimental data

If temperature effects on the sound velocities of iron have been calculated for thermodynamic conditions corresponding to those expected at the inner-core, outer-core boundary (ICB) [START_REF] Sha | [END_REF]53,10], results of calculations have never been directly compared to experimental measurements at the actual density and temperature where data have been collected (i.e. at lower densities and temperature). V P has been much more studied and we are aware of only one NRIXS experiment to measure V S up to 1700 K [1]. The temperature dependence deduced from the experiments for V P are in open disagreement. Lin et al. [1] using NRIXS estimated a temperature decrease of about 0.35 m.s -1 .K -1 at a constant density value of 10.25 g.cm -3 . For the same density and temperature from 300 to 700 K, Mao et al. [75] A c c e p t e d M a n u s c r i p t using IXS give a value almost two times larger: 0.67 m.s -1 .K -1 . From the fit proposed by Sakamaki et al. [3] based on their IXS measurements, we obtain a value of 0.28 m.s -1 .K -1 at this density. In contrast, Antonangeli et al. [80] and Ohtani et al. [81], starting with IXS reported the temperature effects to be negligible up to 1000 K.

From Eq. ( 13) we find a value of 0.28 m.s -1 .K -1 , in agreement with Sakamaki et al. [3] and we do not support the large value found by Mao et al. [75]. The agreement remains very good at high pressure and temperature. They predict a value of 11.85 km.s -1 at 330 GPa, 5500 K) while we obtain 11.93 km.s -1 . For V S , Lin et al. [1] using NRIXS estimated a temperature decrease of about 0.46 m.s -1 . For the same density, Eq. ( 13) gives a similar value.

Conclusion

Using a new method to calculate the vibrational contribution to the free energy we have studied the temperature dependence of physical quantities of hcp iron up to the Earth's core conditions. To assess the validity of the present approach we first compare the results obtained to recent experimental data in the (P , T ) domain accessible to experiments. We obtain an excellent agreement between present calculations of the phonon density of states at room temperature and NRIXS measurements which gives

Figure 1 .

 1 Figure 1. Room temperature phonon density of states of hcp Fe at different pressures. Straight lines: present work at 51 (black), 87 (red), 158 (blue), 245 (green) and 323 (orange) GPa. Open circles: NRIXS measurements [54, 55] at 51 (black), 85 (red) and 151 (blue) GPa.

A c c e p t e d M a n u s c r i p tFigure 2 .

 2 Figure 2. Temperature-dependent phonon dispersion and vDOS of hcp Fe at 300, 1000, 3000, 5000 and 7000 K at constant volume corresponding to a density of 13.06 g.cm -3 .

Figure 3 .

 3 Figure 3. Vibrational entropy (squares) and vibrational heat capacity (circles) as a function of density for hcp Fe on isotherms T = 300 K (filled symbols) and 1000 K (opened symbols). Present work and NRIXS measurements [60, 54, 1]] are respectively in black, green, red and blue.

Figure 4 .

 4 Figure 4. Grüneisen parameter of hcp Fe at T = 300 K as a function of density.The black diamonds are the results of present calculations. The green circles, the red triangles and the blue squares are experimental data from Ref.[55,64,63] respectively. The black, orange, violet and green lines are fits following Eq. (12) for present data and Ref.[61,62,55].

Figure 5 .

 5 Figure 5. Grüneisen parameter of hcp Fe as a function of density and temperature.

Figure 6 .

 6 Figure 6. CTE of hcp Fe at T = 300 K as a function of density. α V is the volumetric thermal expansion while α a and α c are the linear coefficient in the x (y) and z direction respectively. The black symbols are the results of the calculation while the red squares are the data extracted from the NRIXS measurements[60].

Figure 7 .

 7 Figure 7. Temperature dependence of the volumetric CTE as a function of density. The orange diamonds with error bar is the experimental value of Duffy and Ahrens [69] at T = 5200 ± 500 K.

  A c c e p t e d M a n u s c r i p t

Figure 8 .

 8 Figure 8.Theoretical compressional V P and shear V S velocities for hcp-Fe. Black circles, PAW, this study; red squares, PAW, Ref.[START_REF] Price | geophysical and Geochemical Evolution of the Deep Earth[END_REF]; open red squares, PAW, Ref.[53,10]; green triangles, ultrasoft pseudopotentials, Ref.[START_REF] Tsuchiya | Physics of the Earth and Planetary Interiors[END_REF]; violet diamonds, FPLAW, Ref.[START_REF] Steinle-Neumann | [END_REF]; orange down triangles, FPLMTO, Ref.[START_REF] Sha | [END_REF]; brown left triangles, FCD-EMTO Ref.[70]; blue right triangles, FPLMTO, Ref.[71].

5. 3 .

 3 Temperature dependence of the sound velocities 5.3.1. Results and Birch's law We present the temperature dependence of the compressional and shear velocities of hcp iron on Fig. 10. To calculate K S we use the thermal expansion α and the Grüneisen parameter γ obtained using Eq. (

Figure 10 .

 10 Figure 10. Compressional V P and shear V S velocities for hcp-Fe at 0 (black), 300 (red), 1000 (blue), 3000 (green), 5000 (orange) and 7000 K (violet) as a function of density. The compressional and shear velocity of the PREM[6] inner core are shown as large crosses.

  A c c e p t e d M a n u s c r i p t

Figure 11 .

 11 Figure 11. Comparison between the compressional (filled circles) and shear (open circles) velocities obtained form the AIMD and from the fit (see Eq. (13)). The colors of the temperatures are the same as in Figure10.

Figure 12 .

 12 Figure 12. Comparison between the sound velocities obtained with the present fit (filled symbols) and the values obtained in Refs.[53] (open red circles) and Refs.[10] (open black squares).

  A c c e p t e d M a n u s c r i p t A c c e p t e d M a n u s c r i p t

Acknowledgements

We would like to thank M. Torrent for helpful discussions. This research was supported by the PlanetLab program of the Agence Nationale de la Recherche (ANR) grant ANR-12-BS04-0015-04.