A New Description of Pure C in Developing the Third Generation of Calphad Databases

In connection to developing the third generation of Calphad databases a new thermodynamic description is presented for unary carbon. Models used in this work have more physical basis and are valid down to 0 K. The anisotropy in graphite, caused by weak Van der Waals inter-plane forces makes it impossible to fit the heat capacity data by a single Einstein temperature for modelling the harmonic vibration of the atoms. By using multiple Einstein temperatures this problem is solved and a good agreement with the experimental data at both low and other temperatures is achieved. Diamond is modeled using multiple Einstein temperatures due to its specific heat capacities at very low temperatures too, and the two-state model is used for modelling the liquid phase.


Introduction
Modelling thermodynamic properties of materials is an important basis to be able to predict phase transformations that occur during materials production, and such simulations can save time and costs. The Calphad method is a very powerful approach for describing thermodynamic properties by fitting model parameters to experimental and calcuated data. Such thermodynamic descriptions are combined to build thermodynamic databases for multicomponent alloys. In the early years of developing databases, phase transformations at high temperatures were the main focus. Thus, the descriptions, also those of the unaries developed previously and currently used by the Calphad community, [1] are valid only down to room temperature (298 K). Necessity to model phase transformations at low temperatures, e.g. martensitic or bainite transformations, has attracted attention to define new models and develop new databases which have more physical meaning and are valid down to 0 K.
To develop new unary databases, it was suggested at the 1995 Ringberg unary workshop [2] that different contributions to the heat capacity, i.e. electronic, magnetic, vibrational, should be taken into account. One approach based on this suggestion, [2] that is currently possible to implement in the Calphad software Thermo-Calc is the one suggested by Chen and Sundman in 2001 [3] for modelling thermodynamic properties of pure Fe. Although this model is simple, it has so far proven to describe Ni, Cr, [4] Mn, [5,6] Co [7] and Pb. [8] In the present article unary carbon, C, is reassessed.
Carbon, is the most important element of the organic structures and has played a vital role in creating life and its evolution on the planet earth. This element is also an essential part of all Fe-based alloys, which build the foundation of our civilization. For these reasons carbon has been a very attractive subject to study both experimentally and theoretically. In studying phase transformations, having a deep knowledge of its properties and possibility to accurately model these properties will improve the productivity and development of new materials enormously. Several experimental investigations (reviewed in detail in the next section) have shown that carbon has two stable solid phases; i.e. graphite and diamond. Diamond is not a stable phase at ambient pressure and temperatures, and even though graphite is the stable polymorph it is rarely seen in alloys due to kinetic difficulties to precipitate. Instead graphite forms metastable carbides in combination with metallic elements, e.g. cementite in combination with Fe. In graphite, there are different types of inter-atomic bonds between carbon atoms in different directions; covalent bonds in the xy-plane and weak Van der Waals forces in the z-direction. This anisotropy in the atomic interplanes, gives specific properties to graphite which will be discussed in section 2.
Regarding modelling and assessment of carbon, Berman [9] assessed the thermodynamic properties of graphite and diamond using the experimental data available until then. Gustafson [10] assessed thermodynamic properties and the phase diagram of pure carbon using the Calphad method. Gustafson [10] modeled all allotropes, solid, liquid and gaseous species and phase diagram of carbon down to room temperature.
Naraghi et al. [11] suggested a new description for pure C down to 0 K in order to assess the binary system Fe-C, using the approach suggested by Vřeštál et al. [12] Vřeštál et al.'s approach [12] provides a simple and easy amendment to the current SGTE description for the temperature interval 0-298 K. The parameters thus obtained are valid only for the low temperature region and do not have any physical significance.
To keep consistency with other elements assessed in developing the third generation of Calphad databases using the same equations as that of Chen and Sundman [3] for not only the temperature range from 0 to 298 K but also beyond, it was necessary to perform a new optimization for the allotropes of pure C using the same model. However, specific atomic bonds in graphite made it impossible to fit experimental data for the heat capacity at low temperatures by using a single Einstein temperature, which was actually the reason that Naraghi et al. [11] did not use this model. Different approaches were attempted in the present work to solve this problem, among which using multiple Einstein temperatures worked out the best, and will be presented here. The details will be discussed in section 3, but first, a detailed literature review of studies on carbon is presented in the next section.
It should be mentioned that to simplify this first step of modelling neither the gaseous species nor the volumepressure phase diagram are considered in the present work.

Literature Review
As mentioned earlier, carbon, as one the most important elements of organic compounds, earth's crust and steel alloys, has attracted a large attention to investigate its properties. A detailed literature review for these studies, i.e., experimental measurements and theoretical calculations, is presented in this section.

Heat Capacity and Heat Content
The heat capacity of different allotropes of carbon has been measured by different techniques. Worthing [13] calculated this quantity for diamond and graphite by measuring the electric properties in the temperature range of 200-1200 K and DeSorbo [14] measured it for graphite using calorimetry method in the temperature range 13-300 K. Victor [15] measured the heat capacity of diamond from room temperature to 1100 K using the drop calorimetry method and from that, derived the heat capacity and Debye temperature of this phase.
McDonald [16] measured heat content of the extruded graphite using the drop calorimetry method in the temperature range of 341-1723 K. Ishihara [17] extended his work to the temperature region of 1200-2600 K using an adiabatic calorimetry method. Buyco and Touloukian [18] did a detailed review of the heat capacity theory, different experimental methods for measuring it and experimental results from literature. Kaufman and Nesor [19] optimized the Ni-C system, through which they calculated the enthalpy of formation of different end-members including C, e.g. Ni:C-fcc, bcc and hcp. Cezairliyan and Righini [20] and Matsumoto et al. [21] measured the heat capacity and electrical resistance of graphite samples by the use of heating the sample by passing an electric current through it, while Cezairliyan and Miiller [22] measured the same properties using a high speed pyrometer at 1500-3000 K. Chase et al. [23] reviewed and suggested selected values for thermodynamic properties of all allotropes of carbon up to 6000 K.
Rossini and Jessup [24] evaluated the heat of transformation of graphite to diamond equal to 2417 J/mol and an analytic expression for the heat capacity of these allotropes.

Phase Diagram and Triple Point
The triple point and temperature-pressure phase diagram of carbon has been a matter of huge challenges during the decades. This phase diagram was presented by Rossini and Jessup [24] and Ribaub [25] for the first time. It should be mentioned that although the P À T phase diagram of carbon is not assessed in this work, reviewing literature for this is essential in order to extract the melting point of carbon for modelling liquid. A summary of experimental data for the melting point of graphite is presented in Table 1. In the cases that measurement method is not presented (-sign), the investigation refers to analytical assessment or data analysis.
In addition to experimental data, Ferazz and March [34] investigated the phase diagram of carbon based on phenomenological models and Baker et al., [35] Kelley [36] and Wagman et al. [37] discussed it from a calculation point of view. Bundy et al. [38] reviewed the phase diagram data for pure carbon, especially in the low pressure region in which some phase transitions were reported. Gokcen et al. [39] investigated the phase diagram of carbon and found the triple point of carbon at 120 ± 10 atm and 4130 ± 30 K, and Haaland [40] found it to occur at 107 ± 2 atm. Canon [41] and Young [42] investigated the phase diagram of carbon and they confirmed the results from Leidar et al. [30] Leidar et al. [30] did a detailed review on experimental data for high-temperature properties of carbon and drove mathematical expressions for heat capacity based on the most consistent data.
It can be concluded that there is a large scatter in the experimental data for the thermodynamic properties of liquid C, mainly due to the fact that graphite does not melt at ambient pressures. Nevertheless, as was discussed, most of the studies agree on a melting point in the range of 4000-4500 K (extrapolated from high pressure measurements). However, in the previous assessment of carbon by Gustafson [10] the melting point was selected higher, i.e. 4765 K, and the enthalpy of melting was given a low weight in the optimization, mostly because it was believed that the main part is consumed to superheat the liquid.
The scatter in the melting point of carbon made it very difficult to select a suitable value for this quantity. There is no strong evidence for not believing the extrapolated experimental values, unless one decides to sacrifice the accuracy of the unary assessment for the sake of modelling the higher order systems, as it seems to be the main concern in Ref 10. In the present work the initial values in assessent for melting point and entropy of melting for graphite was selected equal to 4130 K and 24.6 J/mol K respectively, according to Gokcen et al. [39] and the enthalpy of melting was adopted from Leidar et al., [30] equal to 125.520 kJ/mol. Due to the effect of superheating liquid, a large uncertainity equal to 21 kJ/mol was considered for the enthalpy of fusion as recommended by Bundy, [28] and this value was not used as an input in the optimization directly, but just to compare with. The heat capacity data for the liquid were taken from Chase et al. [23]

Phonon Dispersion
Drickamer and Lynch [43] measured the phonon dispersion of diamond by applying an elastic neutron scattering technique. Schwoerer-Böhning et al. [44] and Yanagisawa et al. [45] measured this quantity by using a high resolution inelastic x-ray scattering instrument.
Dolling and Cowley [46] used a shell model for diamond's phonon dispersion and Drickamer and Lynch [43] used the Debye model to calculate it. Pavone et al. [47] calculated the phonon frequencies of diamond using the density functional theory (DFT) and the thermal expansion of this phase by applying the quasi-harmonic approach. Their results show a good agreement with the experimental results.
Nicklow et al. [48] measured the frequency of specific normal modes of graphite's lattice, using the coherent inelastic-neutron-scattering technique. Wilkes et al. [49] and Oshima et al. [50] measured the phonon dispersion of graphite at 77 and 300 K by using Reflection Electron Energy Loss Spectroscopy (REELS).
Krumhansl and Brooks [51] modeled the harmonic lattice vibration of graphite using two different Debye temperatures for different temperature ranges. They could fit the experimental heat capacity of graphite in a better way by this approach, which was also tried in the present work. Al-Jishi and Dresselhaus [52] used a Born-von Karman model for calculating the density of states and their results show a good agreement with experimental data. In more investigations [53][54][55][56][57][58][59] the phonon dispersion of graphite was calculated using DFT-based methods.
Mounet and Marzari [60] and Tohei et al. [61] calculated thermodynamic properties of graphite and diamond up to 3000 K using a quasi-harmonic approach. Recently, it has been a major concern of the Calphad community to use the DFT or experimental phonon frequencies in thermodynamic modelling, since it is believed that the best results can be achieved if one starts the modelling from the atomistic scale. [62] However, as mentioned in Ref 62, DFT data are discrete points at specific volume-temperature space and cannot be used directly in the Calphad software; if they are going to be used in such programs, a parametrization is needed to derive other thermodynamic properties. The computational limitations also make it impossible to use the experimental phonon frequencies as input in the thermodynamic modelling. However, if thermodynamic properties of the materials at finite temperature, e.g. enthalpy, heat capacity, etc., are available from DFT techniques, they can be used as input, in the absence of experimental data or to compare with.
In the present work data from Ref 60, 61 could be used for comparison since they have provided heat capacity data at finite temperatures. The agreement between their results and experimental data is very good. However, only data from Ref 60 are used for comparison in the present work (section 4) at low temperatures, to check whether the DFT results can predict a better agreement with the experimental data than the Einstein model. Tohei et al. [61] have not provided their results in the form of tables, but only figures, and the small magnification at low temperatures makes it almost impossible to extract their data.
It should be mentioned that although some ab-initio calculation results, for example work by Grochala, [63] have shown a very large energy for diamond at 0 K, some other works, like Ref 64-66 predict an energy difference of 400-2000 J/mol between diamond and graphite. This value is in a good agreement with SGTE, [1] i.e. 1895 J/mol at 298 K, and consistent with the experimental data for enthalpy and heat capacity. Thus, this value was used in the optimization.

Calphad Modelling
The model suggested by the Ringberg meeting [2] and applied by Chen and Sundman [3] for modelling thermodynamic properties of pure Fe was used in this work. The Gibbs energy is modeled by Eq 1 and 2, for the solid phases, below and above the melting point, respectively: In Eq 1, h E is the Einstein temperature of the solid phase (to model the harmonic vibration of atoms) and is optimized together with the a and b parameters (representing the electronic and anharmonic contributions), E 0 is the total energy at 0 K. Equation 2 is used to properly describe the heat capacity of solid phases above the melting point in the work by Chen and Sundman. [3] In this equation, h E has the same meaning as in Eq 1, while other parameters are calculated in a way that no kink or discontinuity occurs for different thermodynamic properties at the melting point. This treatment also assures that solid phases do not re-stabilize at temperatures above their melting points. It should be mentioned that due to the extremely high melting temperature of graphite, this treatment was not necessary and one single expression, Eq 1, can be used to accurately reproduce the heat capacity of the solid phases.
As was described in the introduction, it is not possible to fit experimental data for the heat capacity of graphite by a single Einstein temperature, due to the different vibration frequencies of atoms in different directions, which give a specific curvature to this property of graphite.
It was already mentioned in the literature review section that Krumhansl and Brooks [51] used two different Debye temperatures for modelling harmonic vibration of atoms in different directions. Their suggested Debye temperatures for graphite, 950 K for displacements normal to the layer planes (h z ) and 2500 K for atom displacements in the x and y planes (h x and h y ), were used as initial values for the Calphad optimization in the present work (in the form of Einstein temperature which is % 0:74h D ). Since the vibration of atoms in the directions x and y are equal, these two can be modelled by the same h E , e.g. the same contribution to the Gibbs energy. However, the harmonic contribution in the z direction is different, as presented in Eq 3-7: By using two Einstein temperatures a good agreement with experimental data can be obtained for the heat capacity, down to 100 K. However, due to the unique density of state (DOS) of this allotrope, more Einstein temperatures are needed to model the harmonic vibration of the atoms over all temperature ranges. It was concluded that by using 5 Einstein temperatures, one can achieve an accurate description for this phase. For that purpose, different Einstein temperatures were weighted differently, for which the weight coefficient was optimized in the fitting against experimental data. These results will be shown and discussed in the next chapter. The same method (multiple Einstein temperatures) was also applied for modelling diamond. Although no such anistropy as for graphite exist in diamond its heat capacity at low temperatures shows a behavior far from the Debye model. Thus, single Einstein or Debye temperature cannot describe properties of this phase. This is discussed in more details in the next chapter too.
For modelling the liquid phase, the two-state model [67,68] was used, shown in Eq 6.
where Gð h am E T Þ describes the harmonic vibration of an ideal amorphous phase. A and B are model parameters related to the zero point energy and other heat capacity contributions of the ideal amorphous phase. DG d models the Gibbs energy difference between liquid-like and solid-like atoms, as presented in Eq 7, in which C, D and E are parameters that can be fitted together with A and B to the experimental data on the enthalpy and entropy of melting as well as the heat capacity. Other terms can be added to this equation, if needed. When the heat capacity data is available for a small temperature range or has large uncertainties, E can be omitted and D constrained to R, the communal entropy.

Results and Discussion
The Calphad description for this optimization is presented in Table 2.
The heat capacity of graphite is shown in Fig. 1(a), and (b). The low temperature region is shown at higher magnification in Fig. 1(b), using 5 (black solid line) and 2 (purple solid line) Einstein temperatures obtained in the present work, compared to the experimental data and the 2D and 3D Debye model, suggested by Krumhansl and Brooks [51] (dotted lines). It should be mentioned that using multiple Einstein temperatures to fit the heat capacity data was shown by Jacobs et al. [69] for different elements and compounds. Arbitrarily weighted multiple Einstein functions have also been used as a mathematical mean by Voronin and Kutsenok [70] to fit the heat capacity data at both low and high temperatures, and in their approach, the ''Einstein'' frequencies are fitting mathematic quantities that have no phyiscal meanings.
The need of at least 5 Einstein frequencies to fit satisfactorily the heat capacity curve over the whole temperature range for graphite suggests that its vibrational density of states is very far from a Debye spectrum. This can of course be confirmed directly by looking at available phonon dispersion relations of graphite or indirectly by examining the dependence of h D ðnÞ on n or h D ðTÞ on T for graphite, where h D ðnÞ is the Debye temperature corresponding to the Debye cut-off frequency that describes the average of the nth moment of the vibrational frequencies over a real phonon spectrum; and h D ðTÞ is the Debye temperature that gives the measured heat capacity at a particular temperature T by applying the Debye model. For an exact Debye solid, all h D ðnÞ are equal to a single Debye temperature, and h D ðTÞ is also a constant. For a real elemental solid, h D ðnÞ and h D ðTÞ vary with n and T, respectively, to a certain extent but normally not so much. [71] Graphite was found to be an extreme exception. First, its h D ðnÞ has no minimum over n, but only increases dramatically from about 490 K to about 2250 K when n goes from À 3 to large values. [72] Second, its h D ðTÞ shows a similar behavior and increases from about 400 K at 0 K to about 1950 K at high temperatures. [61] The huge variations of h D ðnÞ and h D ðTÞ indicate that the actual phonon density of states of graphite must be so complex that the Debye model with one single cut-off frequency is not applicable.
As a matter of fact, using even two Debye temperatures with consideration of the anisotropy of the lattice of graphite [73] has been found impossible to account for the heat capacity data below 100 K by Morgan. [74] The case is almost the same with our trial using two Einstein temperatures.
It might be worth mentioning that Krumhansl and Brooks [51] considered also the elastic anisotropy of graphite. However, instead of using the Debye model for a 3-dimensional lattice, they applied the Debye model for a 2-dimensional lattice. Their approach of a 2D Debye model plus two Debye cut-off frequencies could give a good fit to the heat capacity data below 100 K and above 800 K, but unfortunately not in-between. Figure 2(a) shows the heat capacity for diamond, and Fig. 2(b) shows this quantity compared to the experimental and DFT data. The anisotropy observed in graphite does not exist in diamond. However, although its density of state is not as complex as graphite, the temperature dependence of its Debye temperature below 300 K is far from a Debye model. [61] Thus, the heat capacity cannot be modelled accurately by using either one single Einstein or Debye temperature in this temperature range, while using three Einstein temperatures, a perfect fit can be obtained below 300 K (Fig. 2b). The description of diamond based on using two Einstein temperatures is also shown in this figure, which shows a rather poor agreement with experimental data in this temperature region.
Finally, the heat capacity of liquid, assessed by the twostate model is shown in Fig. 3, compared to the data compiled by JANAF, [23] which are in principle extrapolated from C p data for graphite. Regarding the selected value foe the Einstein temperature of the liquid-amorphous (1400 K), from all experiments available, it is suggested that the Debye temperature of the amorphous state is smaller than that of the crystalline state. The highest Einstein temperature we got for graphite is 1953 K. So, we selected a value smaller than that, i.e. 1953 Â 0:7 ' 1400 K, since the ratio is about 0.7 for known elements.  The liquid description from SGTE [1] is also shown in this figure in red dashed line to compare with the present work. Although both descriptions have the same values above melting point, it can be seen how using two-state model results in a smooth variation with temperature down to room temperature and eventually 0 K in the present work. It is worth mentioning that the liquid two-state model is currently studied carefully to overcome some limitations for modelling liquid phase of different elements. Thus, the current description for C might be revised accordingly in the future, if the model is changed.
The entropy of diamond and graphite from present work at 298.15 K is equal to 2.3 and 5:7 J mol À1 K À1 respectively, which shows a good agreement with the SGTE's descriptions. The enthalpy and entropy of melting from the current description are equal to 101:405:3 kJ mol À1 and 24:1 J mol À1 K À1 as suggested by experimental data in Ref 30 and 39, receptively.

Conclusion
A new Calphad description is suggested for carbon using models that have more physical meaning and are valid down to 0 K. The challenge in modelling graphite is its anisotropy due to the weak Van der Waals inter-plane forces and strong covalent bounds in the basal plane. This phenomena, makes it impossible to fit the heat capacity data by a single Einstein frequency and it is suggested to use multiple (five) Einstein temperatures which is shown to give a reasonable agreement with the experimental data at all temperatures.  3 Heat capacity of liquid, from the present work based on the two-state model compared to the data compiled by JANAF [23] and SGTE [1] (red dashed line) (Color figure online) Diamond, on the other hand, does not have such anisotropy. However, since its Debye temperature has an unusually strong temperature dependence below 300 K, suggesting a severe departure from an ideal Debye behavior, the heat capacity data for this phase cannot be properly modeled by one single Einstein or Debye temperature either. Thus, by using three different Einstein temperatures a good agreement has been achieved at both low and high temperatures.
The two-state model is used for modelling the liquid phase, resulting in a better agreement with experimental data, compared to previous assessments. For the sake of simplicity, the gaseous species of pure carbon and the pressure-temperature phase diagram of carbon are not considered in this work, but can be the subject of future work in developing the third generation of Calphad databases.