A Pisch

A Pasturel

G Deffrennes

O Dezellus

P Benigni

G Mikaelian

Investigation of the thermodynamic properties of Al 4 C 3 : a combined DFT and DSC study

Keywords: Aluminum Carbide, Thermodynamic properties, DFT, SCAN+vdW functional

The question of material stability through the determination of the thermodynamic properties is of fundamental and technological importance to any analysis of system properties in many applications.

For Al4C3, its experimental heat of formation varies widely, from -0.187 to -0.363 eV/atom, which makes it difficult to use such experimental information for any reactivity assessment. Here, we demonstrate that density functional theory (DFT), with the recently developed strongly constrained and appropriately normed (SCAN) functional, is especially powerful in critically assessing these experimental data. In order to have a more complete description of thermodynamic properties of Al4C3, we also determine the temperature dependence of its heat capacity using both the harmonic and quasi-harmonic approximation. In addition, and to select the most efficient approximation, we come back to the experimental background by using differential scanning calorimetry (DSC) to measure this quantity at constant pressure up to very high temperatures, namely 850K.

Introduction

Aluminum Carbide is a refractory material with applications in various technical fields. It is a potential stable phase in the carbothermic reduction of alumina [START_REF] Worrell | Carbothermic Reduction of Alumina: A Thermodynamic Analysis[END_REF]. It may be a product in the reaction of aluminum with silicon carbide [START_REF] Oden | Phase equilibria in the Al-Si-C system[END_REF]. Al4C3 is also present in selected composite ceramics in the Y-Al-Si-C-O system [START_REF] Gröbner | Thermodynamische Berechnungen im System Y[END_REF]. In addition, Al4C3 is an important phase for Al-based metal matrix composites.

Recently, some authors reported its use as reinforcement [START_REF] Besterci | Preparation, microstructure and properties of Al-Al4C3 system produced by mechanical alloying[END_REF][START_REF] Birol | Response to thermal exposure of the mechanically alloyed Al/C powder blends[END_REF][START_REF] Nagendra Naik | Exploration of Mechanical behaviour and Wear Behaviour of Al 4 C 3 Reinforced Aluminium Metal Matrix Composites[END_REF]. However, the formation of Al4C3 is most widely associated with a severe degradation of properties of composites [START_REF] Rodríguez-Reyes | Limiting the development of Al4C3 to prevent degradation of Al/SiCp composites processed by pressureless infiltration[END_REF][START_REF] Etter | Aluminium carbide formation in interpenetrating graphite/aluminium composites[END_REF]. For instance, Al4C3 might form at the matrix/reinforcement interface during processing, dissolving the reinforcements and causing embrittlement. Last but not least, Al4C3 is a semiconducting material and can be used in electronic components such as diodes [START_REF] Kim | Electrical properties of Al/Al 4 C 3 /4H-SiC diodes[END_REF].

In order to control the stability of the phase and the synthesis conditions, a thorough knowledge of the thermodynamic properties of solid Al4C3 is mandatory. Its heat of formation, heat capacity at constant pressure and Gibbs energy were then studied intensively over the last 120 years but surprisingly the reported experimental values display large discrepancies. Indeed the standard heat of formation obtained from different experimental techniques involving variety of calorimetric methods based on combustion, acid-solution and direct-reaction and 2 nd / 3 rd law analysis of vapor pressure measurements [START_REF] Satoh | The Heat of Formation and Specific Heat of Aluminium Carbide[END_REF][START_REF] Meichsner | Beiträge zur Thermochemie des Aluminiums[END_REF][START_REF] Roth | Die Bildungswärmen von Calciumaluminaten[END_REF][START_REF] Meschi | The Dissociation Pressure of Aluminum Carbide[END_REF][START_REF] King | Heat of combustion and heat of formation of aluminum carbide[END_REF][15][16][START_REF] Plante | Dissociation pressure of aluminum carbide using a rotating Knudsen cell[END_REF][18][START_REF] Blachnik | Enthalpies of formation of the carbides of aluminium and beryllium[END_REF][START_REF] Rinehart | Vaporization thermodynamics of aluminum carbide[END_REF][START_REF] Meschel | Standard enthalpies of formation of AlB12 and Al4C3 by high temperature direct synthesis calorimetry[END_REF] display a large range of values, varying from -0.187 to -0.363 eV/atom. More recently, DFT is applied as a computationally efficient component of a methodology capable of accurately determining structural properties and energies of formation of many compounds. Furthermore, this methodology should also be used where experimental data are lacking, i.e. in the prediction of new materials or in the case of experimental discrepancies or controversies. DFT with the PBE [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] generalized gradient approximation (GGA) to its exchange-correlation energy is the most common approach due to its relatively cheap computational cost and reasonable accuracy. However, the theoretical determination at 0K of the heat of formation of Al4C3 using this approximation [START_REF] Jain | Commentary: The Materials Project: A materials genome approach to accelerating materials innovation[END_REF] is considerably less negative than all experimental values, namely -0.013 eV/atom, preventing any use of theory to critically assess the experimental data. Note that corrections to include thermal contributions at T=298K [START_REF] Furukawa | Heat Capacity and Enthalpy Measurements on Aluminum Carbide (A14C3) From 15 to 1173 OK. Thermodynamic Properties From 0 to 2000[END_REF][START_REF] Chase | National Institute of Standards and Technology[END_REF] yields a value of -0.015 eV/atom but do not change the previous conclusion. Such a discrepancy can be attributed to the use of PBE, a common deficiency found in DFT calculations to capture the cohesive energy of graphite accurately [START_REF] Buda | Characterization of Thin Film Materials using SCAN meta-GGA, an Accurate Nonempirical Density Functional[END_REF]. The absence of van der Waals (vdW) interactions in PBE is the major source of error to describe the layered structure of graphite. We also emphasize that Al4C3 has an unusual crystal structure that consists of alternating layers of Al2C and Al2C2 (see Fig. 1) and PBE seems not to be accurate enough to describe its structural and cohesive properties. In this work, we investigate the structural and thermodynamic properties of Al4C3 using the recent SCAN many-body interaction functional [START_REF] Sun | Strongly Constrained and Appropriately Normed Semilocal Density Functional[END_REF]. First, we find that SCAN significantly improves the agreement between the experimental and calculated ground state properties of C-graphite and those of Al4C3. As a consequence, the energy of formation of Al4C3 at 0K is -0.236 eV/atom, a value that is now located inside the range of experimental values. Moreover, we complete our theoretical analysis at T = 0K by including temperature effects through the calculation of the temperature dependence of the heat of capacity of Al4C3 as well as those of Al-fcc and C-graphite. Both harmonic and quasi-harmonic approximations are tested to determine these temperature dependences. Using the quasi-harmonic approximation, we obtain a standard heat of formation equal to -0.250 eV/atom.

As experimental heat capacity data of Al4C3 are only available at temperatures below 350K, we also determine its temperature dependence up to 850K by using the "discontinuous method" as described by Höhne et al. [START_REF] Höhne | Differential scanning calorimetry[END_REF]. This experimental study allows to identify that the quasi-harmonic approximation provides theoretical values in close agreement with experimental data up to 500K.

Theoretical and experimental details

DFT calculations

The ground state properties of pure carbon (graphite, diamond), aluminum (fcc) and Al4C3 were calculated using density functional theory DFT [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF][START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF] and the VASP software package [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] in its most recent version (5.4.4).

The Strongly Conditioned and Appropriately Normed (SCAN) semi-local density functional [START_REF] Sun | Strongly Constrained and Appropriately Normed Semilocal Density Functional[END_REF] was chosen for the calculations. In addition, we include van der Waals interactions to take into account the long range order phenomena in these carbon based materials. The numerical routines from Klimes et al. [START_REF] Klimeš | Van der Waals density functionals applied to solids[END_REF][START_REF] Klimeš | Chemical accuracy for the van der Waals density functional[END_REF] as implemented in the VASP code were used. For Al, the 3p and 3s orbitals and for C, the 2p and 2s orbitals are considered as valence states in the calculations. The energy cutoff for the projector augmented plane-wave bases was set to 800 eV. An automatically generated, gamma centered grid of k-points in the irreducible part of the Brillouin zone was used following the Monkhorst-Pack scheme [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF].

The lattice parameters of all compounds as well as their internal atomic coordinates were fully relaxed. The linear tetrahedron method with Blöchl corrections [START_REF] Blöchl | Improved tetrahedron method for Brillouin-zone integrations[END_REF] was used to calculate the electronic density of states (DOS) for C-graphite and diamond as well as for Al4C3. For Al-fcc, the Methfessel-Paxton method of 2 nd order was used for the relaxation calculation. All relaxations were performed with a convergence criterion of 10 -8 eV/Å for the total energy. Zero Point Energy (ZPE) and heat capacities at constant volume / pressure to access finite temperature properties can be obtained from lattice dynamics theory. The phonon spectra of all compounds were calculated with the frozen phonon (supercell) method using the phonopy code [START_REF] Togo | First principles phonon calculations in materials science[END_REF] coupled to VASP. The generated supercells for the different compounds were: 5x5x2 for C-graphite, 3x3x3 for C-diamond and Al-fcc and 2x2x1 for Al4C3. The convergence criteria for the Hellman-Feynman forces was set to 10 -6 eV/Å to avoid any kind of residual strain in the lattice.

Experimental study

Al4C3 powder was purchased from the Alfa Aesar company. The chemical analysis obtained from the manufacturer is shown in Table 1. As this product is moisture sensitive, the as received bottle containing the carbide powder was opened and further handled in a glove box maintained under a protective atmosphere of purified argon (O2 < 10 ppm, H2O < 5 ppm). The as received carbide powder was first analyzed by XRD. To avoid any potential hydration of the product during the XRD analysis, the aluminum carbide powder was loaded within the glovebox, in an airtight specimen holder equipped with a dome like X-ray transparent cap. During the specimen holder filling, it was not possible to perfectly flatten the surface of the powder because of the risk to trap powder grains in the groove of the O-ring of the holder, which could jeopardize the airtightness. The specimen holder was then transferred to the diffractometer equipped with a copper anode and analyzed in the The low angles (2 < 30°) broad peak is due to the x-ray scattering by the glassy specimen holder dome. There is a good match between the diffractogram of the tested powder and the reference pattern both in term of peak positions and intensities. Moreover, even if some few and low intensity additional peak can be observed, no secondary phases can be detected on the diffractogram.

Therefore, the amount of secondary phases can be considered to be lower than 1 wt.%.

The heat capacity was measured by the so-called "small temperature steps procedure" or "discontinuous method" using the classical three steps method as described by Höhne et al. [START_REF] Höhne | Differential scanning calorimetry[END_REF].

with a DSC 111 calorimeter from the Setaram Company. Further details on the calorimeter and on the experimental protocol can be found in Benigni et al. [START_REF] Benigni | Measurement of the heat capacity of ZnSb by DSC between 300 and 673 K[END_REF]. Four stainless steel crucibles provided by Setaram were used in the experiments. The crucibles were sealed with a crimping tool inside the glovebox. The first crucible was filled with 132.98 mg of Al4C3, the second with 133.571 mg of Standard Reference Material 720 -Al2O3 [START_REF] Uriano | Standard Reference Material 720 Synthetic Sapphire (α-Al2O3)[END_REF], the other two were sealed empty. The temperature program and the processing of the recorded thermograms were performed with the SETSOFT 2000 software provided by Setaram. The parameters in the temperature program were set as follows: 2.5 K temperature step, 1.5 K/min heating rate between each step, and 800 s stabilization time after each step. The heat capacity was measured between 300 and 873 K. A point by point correction factor calculated from the ratio between the reference [START_REF] Uriano | Standard Reference Material 720 Synthetic Sapphire (α-Al2O3)[END_REF] and the apparent heat capacity of α-Al2O3 was applied. The uncertainty of the heat capacity measurements is estimated to be around ±2.5% [START_REF] Benigni | Measurement of the heat capacity of ZnSb by DSC between 300 and 673 K[END_REF].

3. Results and Discussion:

Structural and cohesive properties of C-graphite and C-diamond

The crystal structure of hexagonal C-graphite was determined by Nelson & Riley [START_REF] Nelson | The thermal expansion of graphite from 15C to 800C: Part I. Experimental[END_REF], Baskin & Meyer [START_REF] Baskin | Lattice Constants of Graphite at Low Temperatures[END_REF], Trucano & Chen [START_REF] Trucano | Structure of graphite by neutron diffraction[END_REF], Zhao & Spain [START_REF] Zhao | X-ray diffraction data for graphite to 20 GPa[END_REF] and Howe et al. [START_REF] Howe | Improved crystallographic data for graphite[END_REF]. It crystallizes in the P63mmc space group with two fixed atomic positions: C1 (0 0 1/4) and C2 (1/3 2/3 1/4). C-diamond crystallizes in the cubic Fd-3m structure. Lattice parameter studies were performed by Riley [START_REF] Riley | Lattice Constant of Diamond and the C-C Single Bond[END_REF], Straumanis & Aka [START_REF] Straumanis | Precision Determination of Lattice Parameter, Coefficient of Thermal Expansion and Atomic Weight of Carbon in Diamond 1[END_REF] and Hom et al. [START_REF] Hom | Accurate lattice constants from multiple reflection measurements. II. Lattice constants of germanium silicon, and diamond[END_REF]. DFT calculations of the two allotropes of carbon were also performed using PBE [START_REF] Jain | Commentary: The Materials Project: A materials genome approach to accelerating materials innovation[END_REF]. While the agreement for C-diamond is good, the agreement for C-graphite is poor, especially the calculated c-axis is much longer than the experimental one. As already discussed in the introduction, the absence of vdW interactions in PBE does not allow to reproduce the interlayer coupling in the C-graphite structure. Note this remark holds also for other layered materials [START_REF] Buda | Characterization of Thin Film Materials using SCAN meta-GGA, an Accurate Nonempirical Density Functional[END_REF].

The available experimental [START_REF] Nelson | The thermal expansion of graphite from 15C to 800C: Part I. Experimental[END_REF][START_REF] Baskin | Lattice Constants of Graphite at Low Temperatures[END_REF][START_REF] Trucano | Structure of graphite by neutron diffraction[END_REF][START_REF] Zhao | X-ray diffraction data for graphite to 20 GPa[END_REF][START_REF] Howe | Improved crystallographic data for graphite[END_REF][START_REF] Riley | Lattice Constant of Diamond and the C-C Single Bond[END_REF][START_REF] Straumanis | Precision Determination of Lattice Parameter, Coefficient of Thermal Expansion and Atomic Weight of Carbon in Diamond 1[END_REF][START_REF] Hom | Accurate lattice constants from multiple reflection measurements. II. Lattice constants of germanium silicon, and diamond[END_REF] and theoretical [START_REF] Jain | Commentary: The Materials Project: A materials genome approach to accelerating materials innovation[END_REF][START_REF] Zunger | Self-consistent LCAO calculation of the electronic properties of graphite. I. The regular graphite lattice[END_REF][START_REF] Yin | Ground-state properties of diamond[END_REF] structural data for C-graphite and Cdiamond are reported in Table 2 together with our values of the lattice constants obtained from SCAN calculations including vdW interactions. In the following we will use the term SCAN when referring to our findings. For C-graphite, the agreement with the low temperature data from Baskin & Meyer [START_REF] Baskin | Lattice Constants of Graphite at Low Temperatures[END_REF] is excellent. The c-axis elongation, as observed using PBE, disappears with SCAN. This confirms that the SCAN functional is particularly adapted describe the electronic properties of layered materials such as C-graphite. For C-diamond, an excellent agreement between the calculated and measured lattice parameters is also observed whatever the functional. In Table 3, we compare the cohesive energy of C-graphite and C-diamond computed at T=0K using both PBE and SCAN functionals with the experimental data [START_REF] Kittel | Introduction to solid state physics[END_REF].

Structural properties of Al4C3:

The first determination of the crystal structure of Al4C3 was performed by v. Stackelberg and Schnorrenberg in 1934 [START_REF] Stackelberg | Die Struktur des Aluminiumcarbids Al4C3[END_REF]. Al4C3 crystallizes in the R-3m space group (166). This was confirmed by Cox and Pidgeon [START_REF] Cox | THE X-RAY DIFFRACTION PATTERNS OF ALUMINUM CARBIDE Al 4 C 3 AND ALUMINUM OXYCARBIDE Al 4 O 4 C[END_REF] using powder diffraction. The structure was further refined by Jeffrey et al. [START_REF] Jeffrey | The structure of the aluminum carbonitrides. II[END_REF] and more recently by Gesing and Jeitschko [START_REF] Gesing | The Crystal Structure and Chemical Properties of U2Al3C4 and Structure Refinement of Al4C3[END_REF] using single crystals. The experimental lattice parameters and the atomic positions are summed up in Table 3a and 3b. DFT calculations using PBE were performed by Suetin et al. [START_REF] Suetin | Structural, elastic and electronic properties and formation energies for hexagonal (W0.5Al0.5)C in comparison with binary carbides WC and Al4C3 from first-principles calculations[END_REF], Sun et al. [START_REF] Sun | Prediction on structural, mechanical and thermal properties of Al 4 SiC 4 , Al 4 C 3 and 4H-SiC under high pressure by first-principles calculation[END_REF] and in the frame of the materialsproject.org initiative [START_REF] Jain | Commentary: The Materials Project: A materials genome approach to accelerating materials innovation[END_REF]. We also report results obtained using SCAN. We note that SCAN improves slightly structural values with respect to PBE and the agreement between DFT calculations and experimental values is satisfactory. Again, the improvement mainly concerns the caxis length which is closer to the experimentally observed value when using SCAN.

Heat of formation and heat of capacity of Al4C3

The formation energy of Al4C3 at 0K is obtained by subtracting the weighed sum of the total energy of constituting elements from the total energy of the compound:

∆𝐻 = 𝐸 𝑇 (𝐴𝑙 4 𝐶 3 ) -4/7𝐸 𝑇 (𝐴𝑙) -3/7𝐸 𝑇 (𝐶).
This formation energy is then corrected with respect to the Zero Point Energies (ZPE) estimated from the frequency integration over the vibrational density of states [START_REF] Togo | First principles phonon calculations in materials science[END_REF]:

∆𝐻 𝑐𝑜𝑟𝑟 (𝐴𝑙 4 𝐶 3 ) = ∆𝐻(𝐴𝑙 4 𝐶 3 ) + ∆𝑍𝑃𝐸(𝐴𝑙 4 𝐶 3 )
for which the ZPE correction is calculated as follows:

∆𝑍𝑃𝐸 = 𝑍𝑃𝐸(𝐴𝑙 4 𝐶 3 ) -4/7𝑍𝑃𝐸(𝐴𝑙) -3/7𝑍𝑃𝐸(𝐶)
The heat capacity at constant volume can be determined from the calculated Helmholtz free energy F in the harmonic approximation (HA) using the vibrational density of states as a function of frequency q of the band s:

𝐹(𝑉 0 , 𝑇) = 1 2 ∑ ℏ𝜔(𝑞, 𝑠) + 𝑘 𝐵 𝑇 ∑ 𝑙𝑛 [1 -𝑒𝑥𝑝 ( -ℏ𝜔(𝑞,𝑠) 𝑘 𝐵 𝑇 )] 𝑞,𝑠 𝑞,𝑠
The vibrational entropy and the heat capacity at constant volume are then given by:

𝑆 = -( 𝜕𝐹 𝜕𝑇 ) 𝑉 ; 𝐶 𝑉 = -𝑇 ( 𝜕 2 𝐹 𝜕𝑇 2 )
𝑉 When repeating the HA calculations at several different volumes V to obtain a minimum value of F(V,T), the heat capacity at constant pressure in the quasi-harmonic approximation, QHA, is obtained from:

𝐶 𝑃 = -𝑇 ( 𝜕 2 𝐺(𝑇, 𝑃) 𝜕𝑇 2 ) with G(T,P) = minV[E(V)+F(V;T)+PV].
The standard heat of formation at 298K of Al4C3 is then obtained from 0K calculations by including the ZPE correction and by subtracting the weighed sum of the heat content of constituting elements from the heat content of the compound, namely

∆𝐻 0𝐾 298𝐾 = ∆𝐻 0𝐾 298𝐾 (𝐴𝑙 4 𝐶 3 ) -4/7∆𝐻 0𝐾 298𝐾 (𝐴𝑙) - 3/7∆𝐻 0𝐾 298𝐾 (𝐶).
In Figure 3 we present the calculated phonon Density of States (DOS) of Al4C3 used to determine the temperature dependence of its thermodynamic properties. The derived heat capacity at constant volume (harmonic approximation) and at constant pressure (quasi-harmonic approximation) are plotted in Figure 4 together with the available experimental data (present contribution and data from [START_REF] Furukawa | Heat Capacity and Enthalpy Measurements on Aluminum Carbide (A14C3) From 15 to 1173 OK. Thermodynamic Properties From 0 to 2000[END_REF]). The calculated heat capacity at constant pressure agrees well with the measured data from Furukawa et al. [START_REF] Furukawa | Heat Capacity and Enthalpy Measurements on Aluminum Carbide (A14C3) From 15 to 1173 OK. Thermodynamic Properties From 0 to 2000[END_REF] and from the present experimental determination up to 500K. At higher temperatures, anharmonic contributions become more important and the measured values are higher than the calculated ones.

The calculated ground state energies for Al-fcc, C-graphite and Al4C3 with the SCAN functional are summed up in Table 4. Using the calculated Zero Point Energies (ZPE) and the heat contents from 0K to 298K, the standard heat of formation of Al4C3 is -0.250 eV/atom. Note that the value obtained at T=0K using SCAN, namely -0.236 eV/atom, is much more negative than that obtained with PBE, -0.092 eV/atom. The comparison of our calculated value with the literature information allows now a critical assessment of the available experimental data since the SCAN value is well located in the experimental range. To illustrate this, we report all the experimental data and our theoretical one in Figure 5. We find that the calculated standard heat of formation shows a good agreement with the calorimetric data from Meichsner & Roth [START_REF] Meichsner | Beiträge zur Thermochemie des Aluminiums[END_REF][START_REF] Roth | Die Bildungswärmen von Calciumaluminaten[END_REF] and Rinehart & Behrens [START_REF] Rinehart | Vaporization thermodynamics of aluminum carbide[END_REF] . All other values are considerably more or less exothermic than the calculated value, not to mention again the value using PBE [START_REF] Jain | Commentary: The Materials Project: A materials genome approach to accelerating materials innovation[END_REF] that is off the experimental range. Understanding the origin of the experimental discrepancies is beyond the scope of this work. However, incomplete reaction or metastable reaction products in the combustion or dissolution calorimetry experiments or kinetic blockage of the graphite layer at the sample surface during vaporization may explain the dispersion in experimental values obtained from vapor pressure measurements as indicated by Rinehart & Behrens [START_REF] Rinehart | Vaporization thermodynamics of aluminum carbide[END_REF]. We can also use calculated heat capacity values from the quasi-harmonic approximation to determine the standard entropy. We find a value of 0.924 meV/atom/K (89.12 J/mol/K), which is in excellent agreement with the experimental one [START_REF] Furukawa | Heat Capacity and Enthalpy Measurements on Aluminum Carbide (A14C3) From 15 to 1173 OK. Thermodynamic Properties From 0 to 2000[END_REF][START_REF] Saba | Aluminium carbide. Measurments of the low-temperature heat capacity and correlation of the results with high-temperature enthalpy data[END_REF], i.e. 0.922 meV/atom/K (88.97 J/mol/K).

From the band structure of Al4C3 computed at T=0K, we obtain an indirect gap with a calculated value of 1.342 eV which is equivalent to 1.315 eV using PBE [START_REF] Jain | Commentary: The Materials Project: A materials genome approach to accelerating materials innovation[END_REF]. It is important to mention that the use of SCAN still underestimates band gap energies as this is generally the case for PBE [START_REF] Isaacs | Performance of the strongly constrained and appropriately normed density functional for solid-state materials[END_REF].

Finally, there is a plausible additional correction in relation to the ground state properties of Cgraphite. The heat of reaction for the diamond to graphite transition in carbon was measured by

Rossini & Jessup [START_REF] Rossini | Heat and free energy of formation of carbon dioxide, and of the transition between graphite and diamond[END_REF] using combustion calorimetry and a value of 0.020±0.01eV/atom is reported at 298K. More recently, Kleppa and Hong [START_REF] Kleppa | New applications of high-temperature solution calorimetry 1. Enthalpy of the diamond-to-graphite transformation and enthalpy of formation of Mn5C2 at 1320 K[END_REF] obtain a standard heat of reaction of 0.018±0.010 eV/atom using high temperature solution calorimetry. In addition, Hultgren et al. [START_REF] Hultgren | Selected Values of the Thermodynamic Properties of the Elements[END_REF] published an assessed value of 0.020 eV/atom at 298K while in the compilation of Glushko [START_REF] Glushko | Thermal Constants of Substances[END_REF], a value of 0.019±0.01 eV/atom is given. All these values are in close agreement. At 298K, the standard heat of reaction calculated by SCAN including ZPE corrections and the heat content is equal to 0.070 eV/atom, as shown in Table 5. & Jessup [START_REF] Rossini | Heat and free energy of formation of carbon dioxide, and of the transition between graphite and diamond[END_REF] and the calculated one. Note that the difference with PBE calculations is much more important, namely 0.114 eV/atom.

If we attribute this difference to some limitations of SCAN in the accurate description of the ground state properties of C-graphite as shown in Table 3, the heat of formation of Al4C3 has to be corrected by a value equal to-0.022 eV/atom. The corrected value is then -0.272 eV/atom, which does not modify our conclusions on the assessment of experimental values.

Conclusion

We have shown that the recently proposed SCAN functional including vdW interactions considerably improves the theoretical prediction of the heat of formation of Al4C3 as well as the ground state properties of C-graphite and C-diamond with respect to PBE calculations. We also demonstrated that the calculated heat capacity at constant pressure using the quasi-harmonic approximation is in good agreement with experimental data obtained by DSC up to 500K. Our theoretical approach is also confirmed by the good agreement between the calculated standard entropy and the experimental one. All these results allow us to critically assess the experimental heats of formation of Al4C3 which is characterized by a wide range of available values using different experimental techniques. We thus conclude that SCAN would provide an improved description of the thermodynamic properties of layered materials more generally, at a cost comparable to PBE.

Acknowledgement

Fruitful discussions within the French collaborative network in high temperature thermodynamics GDR CNRS n°3584 (TherMatHT) are acknowledged. Carine Perrin-Pellegrino (IM2NP) is gratefully acknowledged for the XRD analysis of the Al4C3 powder.

Data availability

The raw/processed data required to reproduce these findings were added as supplementary information.

Figure 1 :

 1 Figure 1: Crystal structure of Al4C3

  [20-85°] 2 range, with a step size of 0.013° and a scan step time of 196.7 s. The experimental diffractogram is compared in Figure1to the powder data file 00-035-0799 reference pattern of the Al4C3 aluminum carbide calculated for the copper K( wavelength. An angular correction has been applied to the experimental diffractogram to take into account a slight difference between the z position of the powder diffracting surface and the 2 rotation axis of the goniometer. As already mentioned above, impossibility to level the powder surface during holder filling is the main reason explaining the necessity of this correction.

Figure 2 :

 2 Figure 2: Corrected diffractogram of the as received Al4C3 powder compared to the 00-035-0799

Fig. 3 :

 3 Fig.3: Calculated phonon Density of States (DOS) of Al4C3

Fig. 4 :

 4 Fig.4: Comparison of calculated heat capacity at constant volume (harmonic) and constant pressure

Fig. 5 :

 5 Fig.5: Comparison of heat of formation of Al4C3 data from literature and this work

Table 1 .

 1 Characteristics of the Al4C3 powder according to Alfa Aesar's datasheet, in wt.% .

	Product No	14038
	Purity	99 %+ metals basis
	Lot No	X22B028
	C total	24.6 %
	Al	72.0 %
	Fe	0.1 %
	N	0.8 %
	O	1.4 %

Table 2 :

 2 Lattice parameters and cell volumes for C-graphite and C-diamond.

			a [Å]	c[Å]	V [Å 3 ]	Temperature Reference
	C-graphite	Experimental 2.4612	6.7079	35.189	285K	Nelson & Riley [40]
			2.4579	6.6720		4.2K	Baskin & Meyer [41]
			2.4589	6.7076		297K
			2.464(2)	6.711(4) 35.294	300K	Trucano & Chen [42]
			2.462	6.707	35.207	298K	Zhao & Spain [43]
			2.4617(2) 6.7106(4) 35.219	298K	Howe et al. [44]
		Calculated	2.482	6.60		0K	LCAO, Zunger [48]
			2.468	8.685	45.803	0K	PBE, Jain et al. [23]
			2.450	6.670	34.673	0K	SCAN [This work]

Table 3 :

 3 Cohesive energy (in eV/atom) for C-graphite and C-diamond

		Ecoh (PBE)	Ecoh (SCAN)	Ecoh (Exp.) [50]
	C-diamond	7.737	7.456	7.361
	C-graphite	7.860	7.610	7.374
	From

Table 3

 3

, we can see that PBE cohesive energies are significantly larger than the experimental values in both graphite and diamond while this overestimation is corrected by SCAN.

Table 3a :

 3a experimental and calculated lattice parameters and cell volume for Al4C3 in the hexagonal

	setting (R-3m space group)
		a [Å]	c [Å]	V [Å 3 ]	Temperature Reference
	Experimental 3.325	24.94	238.79	298K	v. Stackelberg &
						Schnorrenberg [51]
		3.30	24.89	234.74	298K	Jeffrey et al. [53]
		3.3355(1)	24.967(3) 240.56	298K	Gesing and Jeitschko [54]
	Calculated	3.352	25.104	244.28	0K	PBE, Suetin et al. [55]
		3.354	25.117	244.69	0K	PBE, Jain et al. [23]
		3.335	24.967	240.50	0K	PBE, Sun et al. [56]
		3.3256	24.844	237.95	0K	SCAN [This work]

Table 3b :

 3b Experimental and calculated atomic positions for Al4C3 in the hexagonal setting

Table 4 :

 4 Calculated ground state properties and derived heat of formation for solid Al4C3.

		ET(0K)	ZPE (0K)	H(298K)-H(0K)	H(298K)
		[eV/atom]	[eV/atom]	[eV/atom]	[eV/atom]
	Al-fcc	-7.699	0.037	0.047
	C-graphite	-10.052	0.175	0.011
	Al4C3	-8.943	0.089	0.024	-0.250

Table 5 :

 5 Calculated ground state properties of C-graphite and C-diamond using the SCAN Then, we get a difference equal to 0.050 eV/atom between the experimental value obtained Rossini

	functional.
		E(0K)	ZPE (0K)	H(298K)-H(0K)	rH(298K)
		[eV/atom]	[eV/atom]	[eV/atom]	[eV/atom]
	C-graphite	-10.052	0.175	0.011
	C-diamond	-9.987	0.186	0.005	0.070


Dr Alexander Pisch
email: alexander.pisch@simap.grenoble-inp.fr

Alain Pasturel

On the heat of formation of ye'elimite Ca4Al6O12.SO4 using density functional theory

Keywords: Modelling, Crystallography/crystal structure, Thermodynamics, Sulfate-based cements

The heat of formation of ye'elimite Ca4Al6O12.SO4 was calculated using density functional theory (DFT). The crystallographic ground state of ye'elimite is the orthorhombic Pcc2 modification in agreement with experimental information. Three different many-body interaction functionals were used for DFT calculations: the generalised gradient approximation (GGA), the Armiento/Mattsson 2005 functional (AM05) and the strongly constrained and appropriately normed function (SCAN). The best results were obtained with the Scan functional and a heat of formation at 298K of -35.5 kJ/mol referred to pure CaO, Al2O3 and CaSO4 was found. All other reported modifications of ye'elimite (cubic, tetragonal, monoclinic) are unstable with endothermic heat of transformations ranging from +18 to +78 kJ/mol.

Introduction

The mineral ye'elimite Ca4Al6O12.SO4 has its name from its natural occurrence in the Hatrurim basin, west of the Dead Sea near Har Ye'elim and Nahal Ye'elim (Israel). The artificial version of this compound is used in calcium-sulfo-aluminate (CSA) cements due to its high reactivity with water. Key characteristics of CSA binders are rapid setting, high early strength and low shrinkage [START_REF] Mehta | Concrete Microstructure, Properties, and Materials[END_REF]. CSA clinkers and cements were developed by Alexander Klein at UC Berkley in the 1960's which explains the second name of ye'elimite -Klein's compound [START_REF] Klein | Calciumaluminosulfate and Expansive Cements Containing Same US Patent[END_REF]. Ye'elimite gained growing interest in recent times due to its use in a novel class of low CO2 binders called BYF (belite-ye'elimite-ferrite) cements. Overall emissions of CO2 are reduced by 20 to 30% in these binders due to a reduced level of limestone in the clinker, a lower burning temperature (typically not exceeding 1573K) and a lower grinding energy [START_REF] Gartner | Industrially interesting approaches to "low-CO2" cements Cem[END_REF].

The crystal structure of ye'elimite belongs to the sodalite family with a generic formula M4(T6O12)X. The sodalite mineral, Na8[Al6Si6O24]Cl2 is the prototype of a series with related structures. The symmetry relations within the sodalite family was reviewed by Fischer and Baur [START_REF] Fischer | Symmetry relationships of Sodalite (SOD)type crystal structures[END_REF]. Ye'elimite exists in two distinct crystal structures with an allotropic phase transition at 737K [START_REF] Depmeier | Aluminate Sodalites -A Family with Strained Structures and Ferroic Phase Transitions[END_REF] or 711K [START_REF] Kurokawa | Phase transformation of Ca4[Al6O12]SO4 and its disordered crystal structure at 1073K[END_REF]. The stable modification at low temperature is orthorhombic of the Pcc2 space group type (Calos et al,1995;[START_REF] Cuesta | Structure, Atomistic Simulations, and Phase Transition of Stoichiometric Yeelimite[END_REF]. Additional crystal structures are reported in the literature. [START_REF] Peixing | The crystal structure of C4A3S Proceedings of 9 th[END_REF] proposed a tetragonal cell of the P-4c2 space group type. [START_REF] Cuesta | Structure, Atomistic Simulations, and Phase Transition of Stoichiometric Yeelimite[END_REF] in their theoretical study using DFT calculations reported that this is not a stable modification of ye'elimite. In fact, it transforms into the stable Pcc2 modification during the relaxation step in the calculation. The high temperature modification is of the I-43m space group type [START_REF] Cuesta | Pseudocubic Crystal Structure and Phase Transition in Doped Ye'elimite Crystal[END_REF][START_REF] Kurokawa | Phase transformation of Ca4[Al6O12]SO4 and its disordered crystal structure at 1073K[END_REF] with some statistically disordered SO4 groups. Both authors claim, that in fact the I-43m modification is only a sub-cell of the I4132 space group of the full system. However, due to the large number of atoms in the potential unit cell, the experimental data do not allow to fit the observed reflections to this space group and to derive corresponding atomic positions. [START_REF] Cuesta | Structure, Atomistic Simulations, and Phase Transition of Stoichiometric Yeelimite[END_REF] in their theoretical study proposed therefore two alternative approximations of the high temperature modification: a tetragonal I-4 cell and a monoclinic Cm cell. Both of them are less stable than the stable orthorhombic Pcc2 cell but present interesting approximations for a deeper understanding of the high temperature disordered cubic modification reported at a later stage.

Thermodynamic data for the heat of formation of ye'elimite are rare and scarce. The knowledge of the Gibbs energy, the enthalpy and the entropy as a function of temperature is important to calculate optimum chemical and mineralogical compositions in industrial clinkers. In addition, high quality thermodynamic data is mandatory to be able to perform high quality heat and mass balance simulations of the clinker production process [START_REF] Meyer | Computation of steady state thermochemistry in rotary kilns: Application to the cement manufacturing process[END_REF] and for modeling the hydration behaviour [START_REF] Martin | Contribution of limestone to the hydration of calcium sulfoaluminate cement Cem[END_REF]. [START_REF] Costa | Heats of formation of 4CaO.3Al2O3.SO3, 4SrO.3Al2O3.SO3 and their mixed crystals[END_REF] determined the heat of reaction for 3 CaO + 3 Al2O3 + CaSO4 = Ca4Al6O12.SO4 using acid solution calorimetry. Pure calcium oxide (CaO), calcium sulfate (CaSO4) and ye'elimite (Ca4Al6O12.SO4) were dissolved in a 4N hydrochloric acid (HCl) aqueous solution at 25°C and the observed heat of dissolution was recorded in an adapted calorimetric set-up. Pure Al2O3 could not be dissolved for kinetic reasons and its heat of dissolution was estimated. The heat of reaction of ye'elimite from calcium carbonate (CaCO3), aluminium oxide (Al2O3) and calcium sulfate dehydrate (CaSO4.2H2O) was measured at 1623K by [START_REF] Ayed | Calcium sulfoaluminate formation Proceedings of the 14 th international conference on cement microscopy[END_REF] using drop calorimetry.

Unfortunately, the authors did not perform a drop of the final product to determine the heat content from room temperature to 1623K. Therefore, the heat of formation at room temperature can only be approximated using tabulated or estimated heat capacity data for all compounds participating in the reaction. [START_REF] Hanein | Thermodynamic data of ye'elemite (C4A3S) for cement clinker equilibrium calculations[END_REF] derived the heat of formation of ye'elimite from sulphur dioxide SO2(g) vapour pressure measurements in the 1290K-1675K temperature range [START_REF] Choi | The sulphur cycle in cement kilns: vapour pressures and solidphase stability of the sulphate phases Cem[END_REF]. The most recent experimental determination of the heat of formation of ye'elimite uses hydrations experiments [START_REF] Skalamprinos | Enthalpy of formation of ye'elimite and ternesite[END_REF]. They derived their value from measurements of the heat of hydration with pure ye'elimite and ye'elimite + anhydrite mixtures. In addition to the experimental values, different estimations exist in the literature. [START_REF] Sharp | Calcium sulfoaluminate cements -low-energy cements, special cements or what?[END_REF] published a value without any further explanation on how it was obtained. [START_REF] Wang | Calculation and verification for the thermodynamic data of 3CaO• 3Al2O3• CaSO4[END_REF] used two different empirical methods based on relative contribution of ions and functional groups to derive the heat and entropy of formation as well as the heat capacity at constant pressure.

A modern way to determine the standard heats of formation of crystalline solids are DFT calculations. Examples for the validity of this approach in the field of oxides are publications from Wattez et al [START_REF] Wattez | First principles calculations of anhydrous cement clinker phases[END_REF] for the heat of formation of C2S, C3S and C4AF using the GGA functional and more recently from Kitchaev et al (Kitchaev et al, 2017) for the enthalpy difference in MnO2 polymorphs using the SCAN functional. The potential of DFT calculations for crystal structure determinations was shown by [START_REF] Cuesta | Structure, Atomistic Simulations, and Phase Transition of Stoichiometric Yeelimite[END_REF] in their theoretical study of the different modifications of ye'elimite using the CASTEP code and the General Gradient Approximation (GGA) functional. They confirmed the ground state orthorhombic structure and reported relative energies for the different modifications. Unfortunately, the authors of the study did not calculate the energies of formation with respect to pure CaO, Al2O3 and CaSO4.

The purpose of this contribution is to determine the heat of formation of ye'elimite with respect to CaO, Al2O3 and CaSO4 using DFT and three different functionals.

Theory

The ground state properties of all reported modifications of ye'elimite were determined using DFT [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF]Kohn, 1964, Kohn and[START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]. The VASP software package [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] in its most recent version (5.4.4) was used for the calculations.

Three different many-body exchange correlation functionals were used for this study: the wellestablished General Gradient Approximation (GGA) [START_REF] Perdew | Accurate and simple analytic representation of the electron-gas correlation energy[END_REF], the hybrid AM05 functional (Armiento andMattsson, 2005, Mattson et al, 2008) and the recently developed Strongly Conditioned and Appropriately Normed (SCAN) semi-local density functional [START_REF] Sun | Strongly constrained and appropriately normed semilocal density functional[END_REF]. For Ca, the 3s, 3p and 4s orbitals, for Al, the 3s and 3p orbitals, for S, the 3p and 3s orbitals and for O the 2s and 2p orbitals are considered as valence states in the calculations.

The cut-off energy for the projector augmented plane-wave bases was set to 800 eV. The gamma centred grid of k-points in the irreducible part of the Brillouin zone was automatically generated following the Monkhorst-Pack scheme [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF].

The lattice parameters as well as the internal atomic coordinates of all modifications of ye'elimite were fully relaxed. The linear tetrahedron method with Blöchl corrections [START_REF] Blöchl | Improved tetrahedron method for Brillouin-zone integrations[END_REF] was used to calculate the electronic density of states (DOS). The relaxations were performed with a convergence criterion of 10 -6 eV/Å for the total energy.

Finite temperature properties such as the Helmholtz free energy F and heat capacities at constant volume / pressure Cv / Cp can be approximated using lattice dynamics theory [START_REF] Ziman | Electron and Phonons Oxford Classic Series[END_REF]. The phonon spectrum of the stable orthorhombic ye'elimite was determined using the frozen phonon (supercell) method. The vibrational modes were calculated using the phonopy code [START_REF] Togo | First principles phonon calculations in materials science Scr[END_REF] coupled to VASP. The convergence criteria for the Hellman-Feynman forces was set to 10 -5 eV/Å to avoid residual strain in the lattice.

The total energy of formation at 0K in kJ/mol is obtained by subtracting the weighed sum of the simple constituents from the calculated ground state energy of the ternary compound:

H(Ca4Al6SO12.SO4) = E(Ca4Al6SO12.SO4) -(3 E(CaO) + 3 E(Al2O3) + E(CaSO4))
The Helmholtz free energy F at constant volume can be calculated in the harmonic approximation (HA) using their phonon density of states as a function of frequency q of the band s:

The vibrational entropy and the heat capacity at constant volume are then given by:

The standard heat of formation at 298K of Ca4Al6SO12.SO4 was calculated using the Cv data of ye'elimite and the tabulated heat content data for CaO, Al2O3 from Janaf [START_REF] Chase | NIST-JANAF Thermochemical Tables[END_REF] and for CaSO4 from [START_REF] Robie | heat capacities and entropies from 8 to 1000K of Langbeinite (K2Mg2(SO4)3), anhydrite (CaSO4) and of Gysum (CaSO4.2H2O) to 325[END_REF].

Results and Discussion

The calculated ground state properties of Ca4Al6SO12.SO4 in the Pcc2 space group are summed up in Table 1 (lattice parameters, cell volumes), Table 2 (atomic positions) and Table 3 (energies of formation). The calculated lattice parameters and cell volume with the GGA functional are close to results reported in the frame of the "materialsproject.org" consortium [START_REF] Jain | The Materials Project: A materials genome approach to accelerating materials innovation[END_REF]. The calculated values from [START_REF] Cuesta | Structure, Atomistic Simulations, and Phase Transition of Stoichiometric Yeelimite[END_REF] differ considerably despite the fact that the authors used also GGA for their calculations. A possible explanation may be that the ground state was not reached in the total energy calculations or that the underlying atomic potentials for Ca, Al, O and S in CASTEP differ from those used in the VASP code. By comparing our calculated data with the experimental values, it is clear that the data calculated with GGA and AM05 are larger than the experimentally reported lattice parameters while the values calculated with SCAN gives more satisfactory results. The main difference to the experimental lattice parameters lies in the fact that the calculated parameters are more orthorhombic-type with an a to b difference of 0.0145 Å instead of 0.0009 nm / 0.00006 nm for the experimental values. As DFT values were calculated at 0K, the difference to the experimental values seems to indicate that the orthorhombic-type character becomes smaller with increasing temperature. This evolution makes sense because the structure tends to higher symmetry with increasing temperature and becomes cubic in the high temperature modification above 737K / 711K. The fully relaxed atomic positions as presented in Table 2 correspond to the calculations with SCAN. The agreement between the experimentally refined positions and the calculated ones is good. This is also the case if one compares mean bond energies between the experimental and calculated crystal structure for which differences lower than 2% are observed. It is important to mention, that such a comparison is biased by the fact, that the experimental data are obtained at room temperature and the calculated one at 0K. The calculated formation energies at 0 K with respect to CaO, Al2O3 and CaSO4 (Table 3) depend also considerably on the used functional. The most negative value was observed with GGA. Our calculated GGA value is in agreement with data reported by [START_REF] Jain | The Materials Project: A materials genome approach to accelerating materials innovation[END_REF] using the same dataset and software for the calculations. We find that AM05 and SCAN yield less negative values. In order to be able to compare the DFT calculated data with the experimental values, one must correct the value from 0K to 298K. To do this, the heat capacity at low temperature is used to calculate the heat content between these two temperatures. Unfortunately, there are no experimental data available in the literature for the heat capacity at constant pressure of orthorhombic ye'elimite at low temperature. However, the heat capacity at constant volume can be derived from phonon calculations in the harmonic approximation. The calculated phonon density of states of orthorhombic Ca4Al6SO12.SO4 (Pcc2) is presented in Figure 1. It shows that there are more vibrational states at lower frequency (acoustical branch) which indicates that the transport properties are preferentially induced by collective displacements with same direction, amplitude and phase. It would be interesting to compare this theoretical results to experimental data using neutron diffraction. The derived Helmholtz free energy, entropy and heat capacity at constant volume in the harmonic approximation are plotted in Figure 2. The heat content between 298K and 0K is 75.7 kJ/mol. This value presents a lower boundary for the heat content due to the fact that it was calculated with heat capacity data at constant volume. The true value at constant pressure may therefore be slightly higher. Combined with the heat contents for CaO, Al2O3 from the JANAF thermodynamic data compilation [START_REF] Chase | NIST-JANAF Thermochemical Tables[END_REF] and CaSO4 from [START_REF] Robie | heat capacities and entropies from 8 to 1000K of Langbeinite (K2Mg2(SO4)3), anhydrite (CaSO4) and of Gysum (CaSO4.2H2O) to 325[END_REF] a correction of 8.1 kJ/mol is applied. The standard heat of formation of stable orthorhombic Ca4Al6SO12.SO4 (Pcc2) at 298K becomes -35.5 kJ/moles. This calculated value is compared to the experimental data in Table 4. Our value agrees with the calorimetric one from [START_REF] Costa | Heats of formation of 4CaO.3Al2O3.SO3, 4SrO.3Al2O3.SO3 and their mixed crystals[END_REF] and the estimations from [START_REF] Wang | Calculation and verification for the thermodynamic data of 3CaO• 3Al2O3• CaSO4[END_REF] and (Sharp et al, 1988). All other experimental values are considerably more negative. However, these values are either measured at high temperature and extrapolated to room temperature with estimated values for the heat capacity of Ca4Al6SO12.SO4 (drop reaction, 2 nd law analysis) or use secondary data for the reaction calculations (hydration). The relative heats of transformation at 0 K of the other structural modifications of Ca4Al6SO12.SO4 as compared to the stable orthorhombic modification are reproduced in Table 5. Cubic (I-43m), Tetragonal (P-4c2, I-4) and monoclinic (Cm) Ca4Al6SO12.SO4 are less stable than orthorhombic (Pcc2) Ca4Al6SO12.SO4. The calculated transformation energies are more endothermic than the observed value of 0.4 kJ/mol for the orthorhombic / cubic transition at 737K [START_REF] Depmeier | Aluminate Sodalites -A Family with Strained Structures and Ferroic Phase Transitions[END_REF].

Conclusions

The structural and thermodynamic properties of various structural modifications of Ca4Al6SO12.SO4 were calculated using DFT. The stable modification at 0K and room temperature is the orthorhombic Pcc2 modification in agreement with experimental information.

Three different functionals were used for the calculations: GGA, AM05 and SCAN. The latter one gives the closest agreement with the experimental data and was used to derive the heat of formation, the entropy and the heat capacity at constant volume. A heat of formation of -35.5

kJ/mol at 298K (with respect to CaO, Al2O3 and CaSO4) and -8401.1 kJ/mol (with respect to the stable elements) was found in agreement the calorimetric value from [START_REF] Costa | Heats of formation of 4CaO.3Al2O3.SO3, 4SrO.3Al2O3.SO3 and their mixed crystals[END_REF]. The reported cubic, tetragonal and monoclinic modifications are all unstable with respect to orthorhombic Ca4Al6SO12.SO4. The calculated transformation energies are more positive than the measured value of 0.4 kJ/mol and therefore none of these structures corresponds to the high temperature modification in agreement with the most recent results from [START_REF] Cuesta | Structure, Atomistic Simulations, and Phase Transition of Stoichiometric Yeelimite[END_REF][START_REF] Kurokawa | Phase transformation of Ca4[Al6O12]SO4 and its disordered crystal structure at 1073K[END_REF].

Figure captions (

 captions Figure captions (images as individual files separate to your MS Word text file).

Figure 1 .Figure 2 .

 12 Figure 1. Calculated phonon DOS for orthorhombic Ca4Al6SO12.SO4 Figure 2. Calculated heat capacity at constant volume and derived standard heat of formation and standard entropy at 298K for orthorhombic Ca4Al6SO12.SO4

Figure 2 .

 2 Figure 2.

Table 1 :

 1 Ground state properties at 0K-lattice parameters and cell volumesof orthorhombic

	Ca4Al6SO12.SO4 (Pcc2). Comparison of experimental and calculated values.
	Experimental	Theoretical
	(Calos et	(Cuesta et	GGA	GGA	GGA	AM05	SCAN
	al , 1995)	al, 2013)	(Cuesta et	(Jain et	(This	(This	(This
			al, 2013)	al, 2013)	work)	work)	work)

Table 2 :

 2 Fully relaxed atomic positions of orthorhombic Ca4Al6SO12.SO4 (Pcc2) as compared to experimental information from[START_REF] Cuesta | Structure, Atomistic Simulations, and Phase Transition of Stoichiometric Yeelimite[END_REF]. Experimental data are printed in italics.

	Atom	Wyckoff position	x	y	z
	Ca1	4e	0.0627(11) 0.2495(14) 0.1554(17)
			0.0586	0.2502	0.1483
	Ca2	4e	0.2549(14) 0.9417(12) 0.2945(16)
			0.2546	0.9354	0.2919
	Ca3	4e	0.2591(16) 0.5301(11) 0.2494(15)
			0.2585	0.5238	0.2399
	Ca4	4e	0.4910(9)	0.2548(13) 0.2196(13)
			0.4926	0.2546	0.2202
	Al1	2d	0.5	0.5	0.2242(15)
			0.5	0.5	0.2287
	Al2	2c	0.5	0	0.2466(12)
			0.5	0	0.2461
	Al3	2a	0	0	0.2292(15)
			0	0	0.221
	Al4	2b	0	0.5	0.2011(11)
			0	0.5	0.1961
	Al5	4e	0.6269(7)	0.1202(7)	0.9963(90)
			0.6277	0.1216	0.9924
	Al6	4e	0.1238(6)	0.6282(7)	0.9519(8)
			0.1197	0.6301	0.9493
	Al7	4e	0.3708(6)	0.6293(7)	0.9760(13)
			0.3703	0.6311	0.9772
	Al8	4e	0.1233(7)	0.8793(7)	0.9787(13)
			0.1215	0.8785	0.9694
	S1	4e	0.2691(7)	0.2693(7)	0.9654(10)
			0.2774	0.266	0.9696
	O1	4e	0.6790(12) 0.2645(16) 0.6111(12)
			0.6664	0.26722	0.6073
	O2	4e	0.8260(8)	0.2042(9)	0.4666(15)

Table 3 :

 3 Calculated energies of formation at 0K referred to CaO, Al2O3 and CaSO4 using GGA,

	AM05 and SCAN functionals.
	Functional	Energy of formation at 0K referred	Reference
		to CaO, Al2O3 and CaSO4
		[kJ/mol]
	GGA	-167.56	(Jain et al, 2013)
		-176.00	(This work)
	AM05	-93.54	(This work)

Table 4 :

 4 Calculated heats of formation at 298K for orthorhombic Ca4Al6O12.SO4 (Pcc2) as compared to the available experimental information. The data are referred to CaO, Al2O3 and

	CaSO4 (column 2) and Ca-fcc, Al-fcc, S-orthorhombic and ½ O2(g) (column 3). All
	supplementary thermodynamic data for the calculation are taken from JANAF (Chase, 1998)
	Technique	Heat of formation	Heat of formation	Reference
		with respect to	with respect to
		CaO, Al2O3 and	Ca, Al, S and ½
		CaSO4 (298K)	O2(g) (298K)
		[kJ/mol]	[kJ/mol]
	Acid solution calorimetry	-49.37	-8415.8	(Costa et al, 1972)
	Drop reaction calorimetry	-81.22	-8447.7	(Ayed et al, 1992 )
	(1623K)
	Heat of hydration	-156.52	-8523.0	(Skalamprinos et al, 2018)
	2 nd law analysis (1290-	-174.52	-8541.0	(Hanein et al, 2015)
	1675K)
	Estimation (Mostafa)	+66.08	-8300.4	(Wang et al, 2011)
	Estimation (Wen)	-26.72	-8393.2	(Wang et al, 2011)
	Estimation (Sharp)	-48.52	-8415.0	(Sharp et al, 1999)
	DFT calculation corrected	-35.48	-8401.1	SCAN (This work)
	with heat content

Table 5 :

 5 Calculated structural modifications of Ca4Al6SO12.SO4 using SCAN compared to the reported experimental data including the heat of transformation at 0K. The heats of transformation are given with respect to the stable Pcc2 modification.

	Structure	Lattice	Cell volume	Heat of	Reference
	type	parameters		transformation at
				0K vs. Pcc2
				[kJ/mol]
	I-43m	a = 0.920 nm	0.77869 nm 3		Experimental (Calos et al, 1995)
		a = 0.9385 nm 0.82661 nm 3		Calculated (Cuesta et al, 2013)
		a = 0.9076 nm 0.74762 nm 3	+75.5	Calculated (This work)
	P-4c2	a = 1.3031 nm	1.55379 nm 3		Experimental (Peixing et al, 1992)
		c = 0.9163 nm
		a = 1.3108 nm	1.58366 nm 3
		c = 0.9217 nm

Acknowledgements

The authors wish to thank the GDR CNRS n°3584 (TherMatHT) for fruitful discussions on this project.


Christian El Hajj

Saïd Moussaoui

Guylaine Collewet

Maja Musse

SPATIALLY REGULARIZED MULTI-EXPONENTIAL TRANSVERSE RELAXATION TIMES ESTIMATION FROM MAGNITUDE MAGNETIC RESONANCE IMAGES UNDER RICIAN NOISE

Keywords: Maximum-Likelihood, Spatial Regularization, Majoration-Minimization, Rician noise, Multi-T2

The extraction of multi-exponential decay parameters from multi-temporal images corrupted with Rician noise and with limited time samples proves to be a challenging problem frequently encountered in clinical and food MRI studies. This work aims at proposing a method for the estimation of multiexponential transverse relaxation times from noisy magnitude MRI images. A spatially regularized Maximum-Likelihood estimator accounting for the Rician distribution of the noise is introduced. To deal with the large-scale optimization problem, a Majoration-Minimization approach coupled with an adapted non-linear least squares algorithm is implemented. The proposed algorithm is numerically fast, stable and leads to accurate results. Its effectiveness is illustrated by an application to a simulated phantom and to magnitude multi spin echo MRI images acquired from a tomato sample.

Introduction

Multiple Spin Echo (Multi-SE) is a technique used in MRI transverse relaxation (T 2 ) relaxometry. Multiple T 2 -weighted images, at a fixed sampling rate, are acquired and an exponential decay curve is fitted to the signal at each voxel, in order to extract the relaxation times T 2 (time constants of the decay curve) and their corresponding amplitude A 0 [START_REF] Adriaensen | Mse-mri sequence optimisation for measurement of bi-and tri-exponential T2 relaxation in a phantom and fruit[END_REF]. The measured decay curve in each voxel can be associated, either to a mono-exponential decay, or, in more general cases, to a multi-exponential decay. In the latter case, extracting all the parameters describing a multi-exponential signal inside each voxel can provide relevant information on the micro-structure or composition of the tissue. In plant tissues, the components of the multi-exponential signal are assigned to principal cell compartments (vacuole, cytoplasm and wall) [START_REF] Van | Mri of plants and foods[END_REF]. T 2 and A 0 maps can therefore be used to access, respectively, information about water status and distribution at the sub-cellular level, representing a smaller scale than the voxel size. This is Thanks to Region Bretagne, École centrale Nantes, Irstea Rennes and GDR ISIS for funding. useful for evaluating cell and tissue structure and functioning, as for example in the study of the mechanical properties of the apple fruit [START_REF] Musse | Monitoring the postharvest ripening of tomato fruit using quantitative MRI and NMR relaxometry[END_REF]. In the brain [START_REF] Laule | Magnetic resonance imaging of myelin[END_REF], three major components can be identified within one voxel: water trapped between the myelin, the intra/extracellular water and the cerebrospinal fluid.

The T 2 values inside each voxel can be modeled either by a discretization of a continuous T 2 distribution [START_REF] Kumar | Bayesian algorithm using spatial priors for multiexponential T2 relaxometry from multiecho spin echo MRI[END_REF], or by imposing a defined number of T 2 components [START_REF] Bouhrara | Incorporation of rician noise in the analysis of biexponential transverse relaxation in cartilage using a multiple gradient echo sequence at 3 and 7 Tesla[END_REF][START_REF] Hwang | Improved myelin water quantification using spatially regularized non-negative least squares algorithm[END_REF]. We will focus our study to the second case, assuming that the number of T 2 decay components is known.

The acquisition of MRI data is carried out in the complex domain, where the real and the imaginary signal components are corrupted by an independent and identically distributed Gaussian noise. Due to various factors, such as magnetic field inhomogeneity, thermal noise, eddy currents , the phase of the complex data presents both time and space dependence. As a consequence, parameters estimation in the complex domain becomes a challenging problem. The common approach used for estimating T 2 and A 0 maps consists of taking the magnitude of the complex data. In the case of single coil acquisitions [START_REF] Henkelman | Measurement of signal intensities in the presence of noise in MR images[END_REF], the noise in the magnitude data is described by a Rician distribution, which, when not properly accounted for, introduces a bias on the estimated parameters [START_REF] Raya | T2 measurement in articular cartilage: impact of the fitting method on accuracy and precision at low SNR[END_REF].

A maximum Likelihood (ML) estimator based on the Rician probability density function can be adopted for the estimation [START_REF] Yokoo | Quantitative R2* MRI of the liver with rician noise models for evaluation of hepatic iron overload: Simulation, phantom, and early clinical experience[END_REF]. Moroever, spatial regularization has been recently proposed to stabilize the solution by taking advantage of low voxel to voxel variability between different image structures [START_REF] Kumar | Bayesian algorithm using spatial priors for multiexponential T2 relaxometry from multiecho spin echo MRI[END_REF][START_REF] Hwang | Improved myelin water quantification using spatially regularized non-negative least squares algorithm[END_REF]. However, direct implementation of spatial regularization has not been possible because of the large amount of memory needed and the considerable computation time. It has been carried out only under Gaussian noise assumption either by dividing the image into overlapping patches [START_REF] Kumar | Bayesian algorithm using spatial priors for multiexponential T2 relaxometry from multiecho spin echo MRI[END_REF] either by carrying the estimation on regions of interest.

In this paper, we propose an efficient algorithm which allows to carry out the estimation on the whole image at once under Rician noise assumption. In order to minimize the effect of the Rician noise we employed a ML estimator and added a penalizing function allowing us to take advantage of spatial regularity and further reduce the effect of the noise.

The first contribution is to carry out the estimation using a Majorization-Minimization (MM) approach coupled with a Levenberg-Marquardt (LM) algorithm adapted to the estimation of all the image parameters simultaneously. The second contribution is to adopt the Armijo linesearch technique in order to ensure the convergence of the algorithm and a bounded maximal step-length to guarantee the parameters positivity. Tests were conducted on simulated phantom and on MRI images of fruit. The performance of the method is illustrated through numerical and visual results.

2 Problem statement MRI data, measured with a Multi-SE sequence, is composed of N τ images, of size [N x × N y ] (N v voxels). The minimum time ∆T E between two successive signal samples τ t and τ t+1 , also called echoes, depends on the image acquisition sequence. The measured signal y j (τ t ) in each voxel j, can be represented by a multi-exponential decay model:

M (τ t , θ j ) = Nc c=1 A 0 (c,j) e - τ t T 2 (c,j) , for t = 1, ..., N τ (1)
where N τ represents the number of echoes and the vector

θ j = [A 0 (1,j) , T 2 (1,j) ... A 0 (Nc ,j) , T 2 (Nc ,j)
] the vector of unknown parameters of length N p = 2N c . N c represents the number of T 2 relaxation times that corresponds to different signal components inside each voxel. By accounting for the Rician probability density function P R , the Maximum Likelihood estimator is obtained by maximizing the following likelihood function:

L (θ) = Nv j=1 Nτ t=1 P R (y j (τ t ) | M (τ t , θ j ) , σ) = Nv j=1 Nτ t=1 y j (τ t ) σ 2 e - y j (τ t ) 2 +M ( τ t ,θ j ) 2 2σ 2 I 0 (∆ jt ) , (2)
where θ = [θ 1 ... θ j ], I 0 is the modified Bessel function of the first kind and order 0 and σ is the standard deviation of the Gaussian noise in the real and imaginary parts.

∆ jt = yj (τt)M (τt,θj ) σ 2
. At high SNR, the Rician probability function approaches the Gaussian distribution. However, low SNR inevitably occurs in relaxometry for high values of τ t which implies the need to take the Rician distribution into account. In the background region no signal is present, the Rician distribution can be approximated by the Rayleigh distribution and σ can thus be estimated using the mean value of the voxels in this region [START_REF] Nowak | Wavelet-based rician noise removal for magnetic resonance imaging[END_REF].

Maximizing the log of the likelihood function expressed by equation ( 2) comes down to minimizing the ML criterion given by:

J M L (y, θ) = Nv j=1 Nτ t=1 M (τ t , θ j ) 2σ 2 2 -log (I 0 (∆ j )) , ( 3
)
where y is the vector containing the measured signals for all the voxels. We can consider that there exists some similarity from one voxel to its neighbor. Spatial regularity can be accounted for by adding a penalizing function to the J M L criterion. For each voxel j we define a neighboring voxels set V j (for example using a window of 3 × 3 voxels) and impose a penalty function ψ on the difference between the parameters:

J reg (θ) = Np p=1 β(p) Nv j=1 k∈Vj ψ (θ j (p) -θ k (p)) . ( 4
)
Here β is a weight vector attributed to the parameters of the penalty function to balance the minimization with the general estimation criterion. By adding the penalizing term to the ML criterion we obtain the following penalized ML criterion:

J P M L (y, θ) = J M L (y, θ) + J reg (θ). (5)

Proposed algorithm

There exists no analytical solution for estimating the parameters using the penalized ML criterion, hence, iterative approach must be adopted. The algorithm that we are proposing in this section allows the use of the Levenberg-Marquardt (LM) algorithm in order to minimize J P M L criterion by adopting a MM approach. The MM algorithm is based on a technique that can be used in order to solve differential optimization problems by using a quadratic surrogate function. Let us consider a convex objective function C(x), at each iteration i of the MM algorithm, there exists a function M aj x|x (i) tangent to C(x) at point x (i) such as C(x) is smaller than M aj x|x (i) for all the values of x [START_REF] Hunter | A tutorial on MM algorithms[END_REF]. M aj x|x (i) is called a tangent majorant of C(x). By using this approach, the minimization of a non quadratic, non convex and multidimensional criterion can be transformed into quadratic criteria. Actually, the J M L criterion is the sum of two terms, a strictly convex quadratic function and a function which is proven to be strictly concave -log (I 0 (∆ j )) [START_REF] Varadarajan | A majorize-minimize framework for Rician and non-central chi MR images[END_REF]. A strictly concave function can be majorized by its tangent. Thus:

-log I 0 ∆ (i) jt - y j (τt) σ 2 R ∆ (i) jt M τt, θ j -M τt, θ (i) j ≥ -log (I 0 (∆ jt )) , (6)
where R(x) = I1(x) I0(x) . The quadratic ML criterion to be minimized at step i of the MM algorithm can be expressed as:

S (i) M L (θ) = 1 2σ 2 Nv j=1 Nτ t=1 yj (τt) R ∆ (i) j -M (τt, θj) 2 (7)
Moreover, by adding the spatial regularization, the voxel wise independence of the solution is lost. We propose to use the MM algorithm in order to obtain a voxel wise separable regularization criterion. Based on Erdogan and Fesslers's work [START_REF] Erdogan | Ordered subsets algorithms for transmission tomography[END_REF], if we consider a convex penalizing function ψ, we can compute a majorant function that reestablishes the variable separability. For each iteration i of the MM algorithm we can establish the following relation [START_REF] Erdogan | Ordered subsets algorithms for transmission tomography[END_REF][START_REF] Pierro | A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography[END_REF][START_REF] Lange | Globally convergent algorithms for maximum a posteriori transmission tomography[END_REF]:

ψ (zj -z k ) = ψ 1 2 2zj -z (i) j + z (i) k + 1 2 -2z k + z (i) j + z (i) k ≤ 1 2 ψ 2zj -z (i) j + z (i) k + ψ 2z k -z (i) j + z (i) k (8)
for z ∈ [θ(1), .., θ(N p )]. The majorant function is thus a sum of the two initial convex functions ψ that are now voxel wise independent at each iteration i of the MM algorithm. We can thus construct a new penalizing criterion by taking advantage of the symmetry of the majorant function:

S (k) reg (θ) = np p=1 β(p) Nv j=1 k∈Vj ψ 2θ j (p) -θ (i) j (p) + θ (i) k (p) (9)
We propose to minimize S (k)

P M L (θ) = S (k) M L (θ) + S (k)
reg (θ) criterion by using the LM algorithm. It is an iterative estimation method that interchanges between the gradient descent method and the Gauss-Newton method. At each iteration l of the LM algorithm the parameters are updated in the direction d LM that is a function of the residual r and the Jacobian matrix J given by :

d (l) LM = J (l) T J (l) + λ (l) -1 J (l) T r (l) (10)
r is expressed such that S (i)

P M L (θ) = 1 2 r 2 2 .
λ is a diagonal matrix that contains the LM parameters for each voxel. It is updated at each iteration, with an initial value that is large enough to take a small step in the steepest descent direction. As the residue r gets smaller, λ is decreased and the LM update direction approaches that of the Gauss-Newton algorithm. Once the direction of descent is obtained, we adopt a backtracking technique based on the Armijo condition in order to compute the step size per voxel j that ensures a sufficient descent and the positivity of the parameters. In the LM algorithm the heaviest computation step is the inversion of the first term in equation [START_REF] Yokoo | Quantitative R2* MRI of the liver with rician noise models for evaluation of hepatic iron overload: Simulation, phantom, and early clinical experience[END_REF].

Actually, by applying the LM algorithm on the whole image at once, we obtain a Jacobian of size N τ × (N p N v ), thus a matrix of size N p N v × N p N v to inverse. However, by using the MM approach, the parameters are independent from voxel to another, thus, this matrix is block diagonal. By taking advantage of this structure, we adopt the inversion of N v blocks of size N p × N p . The proposed Minimization-Majorization for the Penalized Maximum Likelihood (MM-PML) approach is summarized by the following algorithm: Initialize θ (0) ; for i ← 0 to i max or until convergence do Compute θ (i)

j (p) + θ (i) k (p) for all j and k ; Compute M (τ t ) = M (τ t )R D (i) jt ; for l ← 0 to l max do Update λ (l) and Compute d (l) LM ; Compute the maximum step size α M ; set α (0) = α M ; for k ← 0 to k max or until Armijo condition is satisfied do Set α (k+1) ← α (k) * c with 0 < c < 1; Compute the Armijo condition end update θ (l+1) = θ (l) + αd (l) LM ; end Set θ i+1 = θ lmax end Output θ (i+1)
Algorithm 1: The MM algorithm for the penalized approach

Results

The evaluation was divided into two parts, the first one is based on a simulated phantom and the second one consisted in in applying the algorithm on experimental tomato MRI data. In order to simulate a phantom, a circular phantom was constructed, on an image of 128 × 128 voxels, with an outside ring and 9 inner disks with different sizes. A tri-exponential model was generated for the different parts of the phantom with settings that are close to parameters typically found in tomato fruits as seen in figure 5. We then added a Rician noise with a σ = 7 in order to obtain a SNR value close to that found in MRI images.

For the experimental MRI settings, we used a Multi-SE sequence on a 1.5T MRI scanner (Magnetom, Avanto, Siemens, Erlangen, Germany), with inter-echo spacing (∆TE) of 6.5 ms, bandwidth of 260 Hz/pixel, 512 echoes per echo train and a repetition time of 10s. The median planes of fruit (transverse section at middle height of fruit) were imaged with a total of 128 × 128 voxels and a slice thickness of 5 mm, resulting in voxel size of 1.19 × 1.19 × 5 mm 3 . Two scans were made, the first one with 32 accumulations in order to obtain higher SNR and to serve as a reference, and the second one with only 1 accumulation in order to verify the stability of the algorithm at the lowest SNR.

In order to implement the regularization, we had to choose the penalizing function ψ and the penalizing weight vector β. The choice of the penalizing function affects the inter-tissues smoothness, whilst the penalizing weight affects the degree of smoothness on the whole image. In this study, we used an L1-L2 penalizing function with a weight β chosen by trial and error so as to preserve a good separation between edges. As seen in figure 1, the results on the phantom show that the reconstruction was made with satisfying edge preserving and that using the proposed MM-PML method allows to detect the different structures even in low contrast situations. For a sake of comparison the non-regularised version (β=0) is also represented. This illustrates the importance of the regularization without which the maps were very noisy with even no contrast between certain regions. Furthermore, we obtained a normalized root mean square error (NRMSE) of 4.47% per voxel and per parameter compared to a value of 22.85% for the non regularized version. The NRMSE was computed as:

N RM SE = 100 1 NvNp Nv j=1 Np p=1 (θj(p)-θ * j (p)) 2 θ * j (p) 2
, with θ * j (p) the reference parameters at voxel j for parameter p. In figure 2 the scan with 32 acquisitions had an SNR of 687.7 whilst the images acquired with only 1 acquisition had an SNR of 145.74 (here SN R = Nv j=1 yj (τ1) σNv

). We used the parameters estimated from the high SNR image as a reference in order to compute the NRMSE which was equal to 5.56%. Also, visually we were able to validate the robustness of the method by comparing the reconstructed maps from low and high SNR data for the different parameters which were quite similar. The evolution of the J P M L criterion during each step of the MM algorithm is shown in figure, it is clearly shown that by decreasing the surrogate criterion S M L the descent of the majorized criterion is guaranteed which validate the choice of the majorant functions.

Conclusion

In this paper we proposed an efficient method for the reconstruction of T 2 and A 0 maps from noisy MRI magnitude data. We showed that visually and numerically we obtained good results both on a simulated phantom and on real scans of a tomato. We were able, for the first time as far as our knowledge, to reconstruct multi-exponential maps with a voxel level information while incorperating both the Rician noise and the special regularity of the parameters. The proposed MM-PML algorithm showed both stability and accuracy. Furthermore, the estimated parameters can be used as features of a classification algorithm in order to characterize different fruit tissues.

Fig. 4 .

 4 Fig. 4. The evolution of the J P M L criterion (in blue) across the external MM steps (orange points) and the internal LM steps (blue points)

Fig. 5 .

 5 Fig. 5. T 2 and A 0 values used for generating the phantom.

Fig. 1 .

 1 Fig. 1. Reference and estimated A 0 and T 2 maps on the phantom using a non-regularized (second columns) and the proposed regularized algorithm (third columns).

Fig. 2 .

 2 Fig. 2. A 0 and T 2 maps reconstructed from MRI images of a tomato acquired with 32 acquisitions and 1 acquisition.

  A 01 T 21 A 02 T 22 A 03 T 23

	• 95	88 459 356 1015 716
	• 96	76 410 433 1024 870
	• 71	50 339 218 908 627
	• 109 78 375 303 1004 685
	• 108 50 482 202 756 508