# Comprehensive assessment of empirical potentials for molecular dynamics simulations of chromia

Paul Fossati

► To cite this version:

Paul Fossati. Comprehensive assessment of empirical potentials for molecular dynamics simulations of chromia. Solid State Ionics, 2025, 420, pp.116781. 10.1016/j.ssi.2024.116781. cea-04904973

HAL Id: cea-04904973

https://cea.hal.science/cea-04904973v1

Submitted on 21 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Comprehensive assessment of empirical potentials for molecular dynamics simulations of chromia

Paul C.M. Fossati*

Université Paris-Saclay, CEA, Service de Recherche en Corrosion et Comportement des Matériaux, 91190 Gif-sur-Yvette, France

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

Cr₂O₃

Molecular dynamics

Defects

Atomic-scale modeling

Chromia

Diffusion

Empirical potential

## ABSTRACT

The importance of Cr₂O₃ (chromia) lies in its ability to form protective layers on chromium-rich metallic alloys, which is utilised in the industry for constructing corrosion-resistant austenitic steels and nickel-based alloys. A better understanding of large defects in Cr₂O₃ is critical because these defects play a crucial role in the growth kinetics of the protective chromia scale, influencing self-diffusion mechanisms, dominant defect types, and diffusion behaviour, all of which can influence the performance and durability of chromium-based alloys. This study presents a comprehensive evaluation of various empirical potentials for simulating the properties of Cr₂O₃ in order to determine the best model to use to simulate extended defects. The assessment is focused on structural, thermodynamic, elastic, point defect, and grain boundary characteristics. An extensive literature review was conducted to compile a dataset for validating the available empirical potentials for Cr₂O₃ from the literature. The evaluation of these empirical potentials provides valuable insights into their strengths and limitations, enabling researchers to make informed decisions when selecting appropriate potentials for simulating various properties of Cr₂O₃. This study's findings contribute to the ongoing efforts to improve the accuracy and reliability of computational materials science methods for predicting the behaviour of complex oxides.

## 1. Introduction

Chromia (Cr₂O₃) is an industrially and technologically significant oxide since it forms protective layers on chromium-rich metallic alloys. In the nuclear context, this phenomenon is utilised for formulating austenitic steels resistant to corrosion, containing between 17 and 20 mol% chromium, and used as structural materials [1] or in nickel-based chromia-forming alloys such as Inconels used in steam generators [2]. In nature, Cr₂O₃ occurs as the rare mineral eskolaite [3]. Different chromium oxides can form on metallic chromium due to the various degrees of oxidation possible for Cr cations, ranging from 2+ in CrO to 6+ in CrO₃. The thermodynamically stable oxides at ambient pressure are Cr₂O₃ up to 2300 °C, and Cr₃O₄ between 1650 °C and 1705 °C [4]. The other chromium oxides, CrO, CrO₂, and CrO₃, are metastable under ambient conditions. CrO decomposes into Cr₂O₃ at 967 °C [5], CrO₂ decomposes above 400 °C when pure, and at lower temperatures when in contact with Cr₂O₃ [6], and CrO₃ decomposes into various oxides from 250 °C, ultimately forming Cr₂O₃ around 500 °C [5,7].

Chromia is a stoichiometric compound up to its melting point, with minimal deviation. For instance, observations indicate the presence of

10⁻⁴ chromium vacancies per mole of Cr₂O₃ under an oxygen partial pressure of 1 atm and at 1373 K [8]. This manifests as a near-linear appearance in the Cr—O phase diagram in Fig. 1, which was generated using the TAF-ID thermodynamic database [9] incorporating the adjusted Cr—O binary system by Kjellqvist et al. [10].

Non-congruent melting is observed between 2540 K and 2675 K [11–13]. This variation can be attributed to reduction of Cr³⁺ to Cr²⁺ depending on the atmosphere [5]. This makes the melting point dependent on the oxygen partial pressure Pₐ₂ during experiments, which is not controlled consistently across different studies. A melting point of 2603 K [5,14] is used in recent thermodynamic models [10,15,16].

The protective function of the chromia scale on metallic alloys relies on the growth kinetics of this layer [17]. A crucial aspect is self-diffusion in Cr₂O₃ [18], the mechanism of which remains only partially understood [19,20]. While It is generally accepted that growth involves opposing diffusion of chromium from the metal substrate to the atmosphere and oxygen from the atmosphere to the metal [21], the type of defects involved and their diffusion pathways remain uncertain. Specifically, grain boundaries, dislocations, and local stresses are expected to have significant effects, which are not well quantified [22]. Another

<footnote>* Corresponding author.
E-mail address: paul.fossati@cea.fr.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Cr—O phase diagram showing the small composition range for Cr₂O₃.</figcaption></figure>

factor is linked to the gradient of chemical potential and oxygen partial pressure within the oxide layer [21]. It can influence the dominant defect types [23] and diffusion behaviour, thereby affecting the growth kinetics of the chromium oxide layer [24].

Questions about intrinsic defect properties and self-diffusion prompted numerous atomic-scale simulation studies. Calculations were performed to estimate diffusion based on possible paths of point defects [23,25–28]. Electronic structure calculations allow for the study of various charge states possible for point defects and their impact on diffusion paths [23,28]. However, these techniques have limitations, as they are unable to reproduce complex collective movements or anharmonic effects that may be significant at high temperatures. The characteristic dimensions of density functional theory (DFT) simulations also do not enable the study of complex nanostructures such as grain boundaries or dislocations. These challenges can be overcome by the use of empirical potentials, which enable molecular dynamics simulations at larger scales of time and space for a more reasonable computational cost. However, such simulations depend on the availability of high-quality potentials. Several potentials have been published based on pair interactions, sometimes combined with polarisable ion models or many-body terms such as the embedded atom model (EAM). Some studies compare certain potentials [29,30], but these are limited by the properties considered (diffusion over a restricted temperature range and elastic properties) and only consider a select number of potentials available in the literature. The work presented in this paper consists of an in-depth study of the characteristics of Cr₂O₃ predicted by existing empirical potentials to determine their quality and their suitability to study self-diffusion at high temperatures in complex structures including grain boundaries. The characteristics considered for this assessment are crystal structure, thermophysical parameters (enthalpy, heat capacity, thermal expansion) and elastic properties at room temperature, as well as point defect properties (formation energy and diffusion pathways) and grain boundary energies.

This article is organised as follows. Section 2 details the simulation techniques and the empirical potentials used in this work. Then, Section 3 is a literature review that includes a description of which properties chosen to validate the potentials. The complete tables with the data compiled from the literature can be found in the appendix. The results and a discussion of the quality of the different potentials follows in Section 4.

## 2. Methods

### 2.1. Simulations

Atomic-scale calculations in the literature are based on either DFT or empirical potentials. DFT provides a better description of the material’s chemistry and defects at low temperatures, allowing for simulations of electronic defects and magnetic structure that are inaccessible to empirical potentials. However, its high computational cost makes it unsuitable for simulating defects on an extended scale, such as grain boundaries or dislocations, over long simulation times needed to observe diffusion events. For these reasons, empirical potentials continue to be of interest in the study of Cr₂O₃.

The properties of perfect chromia were calculated using molecular dynamics simulations with the potentials described later in Section 2.2. The boxes were orthorhombic supercells, containing 138,240 atoms and measuring approximately 120 × 100 × 110 Å. To cover a temperature range from 300 K to 3000 K, simulations were conducted at successive temperatures separated by 65 K increments. A perfect crystal was relaxed at each temperature and using all potentials. At the end of this relaxation, basic properties (lattice parameters and enthalpy) were measured. Other properties, such as the isobaric heat capacity Cₚ and thermal expansion coefficients α, defined by

$$C_p = \frac{dH}{dT}$$

and

$$\alpha_a = \frac{1}{a} \frac{da}{dT} = \frac{d\ln(a)}{dT}, \quad \alpha_c = \frac{d\ln(c)}{dT}, \quad \alpha_v = \frac{d\ln(V)}{dT}$$

were obtained from the enthalpy and lattice parameters, respectively, using centred finite difference schemes. The LAMMPS code [31] was used for dynamic simulations at different temperatures. The steps to obtain the equilibrium properties of Cr₂O₃ included a fixed-temperature relaxation (NVT) for 10 ps to equilibrate ion velocities. This is followed by a temperature- and pressure-controlled relaxation (NVT) for 250 ps, during which the mean dimensions of the simulation box were calculated. The size of the simulation boxes are adjusted to these averages before a final controlled-temperature relaxation. Temperature and stress were controlled using Nosé-Hoover thermostats [32,33] and a Parinello-Rahman enthalpy-based barostat [34], implemented in LAMMPS as an integrator described by Tuckerman et al. [35]. The time constants are respectively 0.01 ps and 1.0 ps and the integrator time step is 1 fs.

Certain properties, such as the elastic constants and bond lengths, were calculated from optimizations of the geometry of the conventional unit cell. The static simulations were performed using GULP [36], and the energy minimization used the BFGS algorithm. To represent the effect of temperature, the structures were expanded so that their lattice parameters matched the values calculated during dynamic simulations at 300 K and kept fixed during the geometry optimization. Internal parameters, related to the fractional coordinates of chromium and oxygen ions, were also obtained through these static calculations. This method facilitates the measurement of bond lengths compared to using radial distribution functions directly obtained from dynamic simulations. Indeed, as shown in Fig. 2, the peaks corresponding to different bonds overlap, even at 300 K. This is particularly true for O—O bonds, whose first three peaks are difficult to distinguish. However, this figure shows that the bond lengths obtained after energy minimization correspond well to the maxima of radial distribution functions, indicating that the structure of these static calculations is representative of the average structure in dynamic simulations at 300 K.

The same simulations, based on energy minimisation with lattice parameters corresponding to the structure at 300 K to represent the effect of temperature, were used to calculate elastic constants. This setup is designed to produce elastic constants that are directly comparable to
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Pair distribution functions in Cr₂O₃ from MD simulations at 300 K using the Chartier potential; vertical lines indicate the neighbour distances from geometry optimisation of unit cells.</figcaption></figure>

experimental measurements performed at 300 K. Test simulations were done using the Chartier potential to validate this method. Elastic constants calculated using the Born tensor and stress fluctuations over the course of a 0.5 ns NVE simulation typically differed from those obtained from static calculations by less than 2 GPa. This difference is small compared to the typical accuracy of the elastic constants predicted using empirical potentials. Using static calculations therefore causes in this case an acceptable decrease of accuracy in exchange of smaller computational costs and simulation times.

Defects were introduced in orthorhombic supercells of 8640 atoms, measuring approximately 40 × 50 × 40 Å. The geometry of these supercells was optimised at a pressure of 0 GPa in order to reach a stable structure for each potential. Then, the defects were introduced by removing or adding atoms in the centre of the supercells. The structures were relaxed by minimising their energies whilst keeping the shape of the supercells fixed. Frenkel pairs were created with all the possible structures in which the vacancy is in a 5th neighbour position of the interstitial or closer. The Frenkel pair with the lowest energies were considered stable and referred to as the “close” structures, in which both defects interact strongly. Other simulation boxes were set up in which the two defects were created as far from each other as they could be. This is referred to as the “far” structure, in which the separation between the interstitial and the vacancy is large enough for their interaction to be negligible. A similar approach was used for Schottky defects. Because of the fixed-charge nature of the empirical potentials, simulations involving isolated charged defects had a small deviation from electrostatic neutrality. The validity of these calculations was checked by comparing their results with those from neutral supercells containing “far” clusters. Another consequence of the fixed charges is that the standard reference states of metallic chromium and molecular oxygen cannot be used when calculating defects energies. The reference energies chosen here were those of oxygen and chromium in perfect Cr₂O₃. This means that the formation energies of charged defects cannot be compared with DFT values from the literature based on standard reference states. The energies of neutral defects can still be directly compared because they only involve the Cr₂O₃ reference state.

Diffusion pathways were investigated using the climbing image nudged elastic band (CI-NEB) method [37] implemented in LAMMPS, using the defect boxes described previously. In these calculations, series of 27 successive images represented the reaction pathway, with inter-image forces parallel and perpendicular to 30 eV/Å².

The quality of the potentials is evaluated quantitatively for the type of data known with a high accuracy, such as structural parameters. A relative error is calculated for each value obtained by empirical potentials, using the values mentioned in Section 0 as references. Then, the potentials are evaluated separately for each category of data (lattice parameters, elastic constants, etc.) To do this, the mean absolute values of relative errors are calculated. The best potentials for each category are those with the smallest average errors. The potentials are assessed qualitatively for data or properties that are too uncertain or scarce in the literature.

## 2.2. Potentials

Several empirical potentials are available in the literature. However, comparative studies to determine the best potential for studying diffusion in Cr₂O₃ remain rare. To address this issue, I considered most of the published potentials I could find, which are described in this section. A number of these potentials, particularly the older ones, use a polarizable ion model (shell model). This type of model poses technical problems during molecular dynamics simulations, particularly at high temperatures or in the presence of structural defects. They also have an additional cost in terms of computation time compared to rigid ions. For this reason, and because dynamic simulations in temperature are performed without shells, polarization was ignored for all relevant potentials. Removing the shells may significantly change the properties of the simulated material. I checked that the structural parameters and elastic constants were consistent between the same potential with and without shells. However, there were significant discrepancies in the phonon densities of states, and I did not compare temperature-dependent thermodynamic properties. The conclusions reached in this article do not concern the potentials with shells included. In all simulations,
electrostatic interactions were calculated using the Wolf summation method with a damping coefficient of 0.3 Å⁻¹ [38].

### 2.2.1. Atkinson

It is a potential initially fitted to simulate doping in corundum-structured materials Al₂O₃, Cr₂O₃ and Fe₂O₃ [39]. It is based on an older model for the simulation of spinels ZnCr₂O₄ [40] with oxygen description from an earlier model [41]. It is a polarisable ion model, with formal charges and pair interactions represented by a Buckingham function

$$E_{a\beta} = \frac{q_a q_{\beta}}{4\pi\epsilon_0 r_{a\beta}} + A_{AB} \exp\left(-\frac{r_{a\beta}}{\rho_{AB}}\right) - \frac{C_{AB}}{r_{a\beta}^6}$$

where α and β are two distinct atoms of element A and B respectively. The parameters of this expression are given in Table 3, and those of the polarizable ion model in Table 4. This potential describes Cr₂O₃ with lattice parameters of 5.038 Å and 13.328 Å for a and c respectively. The formation energies of Frenkel pairs and Schottky defects calculated using this potential [39] are given in Table 7.

### 2.2.2. Barrera

This potential was designed to simulate the behaviour of Cr₂O₃ at high temperatures [42] and is based on a rigid ion model with pair interactions based on a Pauling-type model of the form

$$E_{a\beta} = \frac{q_a q_{\beta}}{4\pi\epsilon_0 r_{a\beta}} + \frac{e^2}{n(\sigma_A + \sigma_B)} \left(\frac{\sigma_A + \sigma_B}{r_{a\beta}}\right)^n$$

The parameters of this model are given in Table 1. The partial charges were chosen based on Hartree-Fock calculations [43]. This potential does not appear to have been used in other studies.

### 2.2.3. Chartier

This potential was fitted to simulate spinels (Ni,Fe)Cr₂O₄ [44] and is derived from the Morooka potential, which was initially designed for simulating MgO-Al₂O₃ mixtures [45]. It retains the basic features of this potential: rigid ions, partial charges, and pair interactions. Non-Coulombic interactions are represented by Buckingham and Morse potentials, resulting in a general analytical form

$$E_{a\beta} = \frac{q_a q_{\beta}}{4\pi\epsilon_0 r_{a\beta}} + A_{AB} \exp\left(-\frac{r_{a\beta}}{\rho_{AB}}\right) - \frac{C_{AB}}{r_{a\beta}^6} + D_{AB} \left\{ e^{-2\rho_{AB}(r_{a\beta}-r_0)} - 2e^{-\rho_{AB}(r_{a\beta}-r_0)} \right\}$$

The values of the parameters are in Table 3. The charges for chromium and oxygen ions are respectively 1.8 e and -1.2 e. Although this potential contains a description of Cr³⁺ and O²⁻ species and their interactions, its performance and the realism of the properties it predicts for Cr₂O₃ were never verified. In particular, it was not considered in the available potentials comparisons that consider diffusion behaviour [29,30].

### 2.2.4. Grimes

This potential was designed to study the solubility of Zn in spinels ZnCr₂O₄ [40] and later for CeO₂ doping with certain transition metals [46]. It is very similar to the Atkinson potential, for which it served as a basis [39]. It shares the same analytical form (polarizable ions, formal charges, and Buckingham interactions). The oxygen-oxygen interaction was taken from Gale et al. [47], who adjusted it to simulate Al₂O₃. The parameters of this potential are given in Table 3 and Table 4. This

<table><caption>**Table 1**
: Parameters of the Barrera potential.</caption><thead><tr><th>Element</th><th>σ</th><th>q</th></tr></thead><tbody><tr><td>O</td><td>0.755</td><td>-1.55</td></tr><tr><td>Cr</td><td>1.24</td><td>2.325</td></tr></tbody></table>

potential was used to study diffusion in Cr₂O₃ [25,29]. The Grimes and Atkinson potentials are very similar, which is why the Grimes potential will not be used in the rest of this article.

### 2.2.5. Lawrence

This is a potential derived from the Lewis potential [48]. It is a polarisable ions model with Buckingham interactions. It retains the chromium-oxygen interactions of the Lewis potential [49], but introduces changes intended to improve defect properties. It was used to simulate free surfaces [48] and grain boundaries [50]. The parameters of the Buckingham terms and polarizable ions are given in Table 3 and Table 4 respectively. There is another potential published by Lawrence et al. in 1987 [51]. It is closer to the Lewis potential, the only changes being in the core-shell interactions parameters. This version was not considered in this work because of this similarity.

### 2.2.6. Lewis

This potential allows the simulation of more than twenty different oxides [49]. It is a model of polarizable ions with formal charges and Buckingham interactions. The O—O interactions are taken from an older potential by Catlow [52]. The parameters of the Buckingham terms and polarizable ions are given in Table 3 and Table 4 respectively.

### 2.2.7. Owen

This potential was designed to study chromium diffusion in crystalline and amorphous UO₂ [53]. It was adjusted to reproduce the lattice parameters of Cr₂O₃ and, to a lesser extent, its elastic constants. It is an extension of the CRG potential for actinide oxides [54], which is based on a rigid ions model with partial charges. It relies on pair interactions with Morse and Buckingham terms, whose analytical form is identical to that of the Chartier and Morooka potentials. The parameters for these functions are in Table 3. Their values are similar to those of the Chartier potential, although the O—O interactions are constrained to remain compatible with the CRG potential. In addition to these pair interactions, the Owen potential also includes an EAM-type contribution. The energy due to these interactions is written as

$$E_{a\beta}^{\text{EAM}} = -G_A \sqrt{\sum_{\beta} \sigma_B(r_{a\beta})}$$

where the function σB corresponds to the electronic density induced by the presence of the atom β with the element B. It is written

$$\sigma_B(r_{a\beta}) = \left(\frac{n_B}{r_{a\beta}^6}\right) \times \frac{1}{2} (1 + \text{erf}(20(r_{a\beta} - 1.5)))$$

The parameters of these two functions are given in Table 2. Given that this potential is very recent, it has not been used in any other studies available in the literature.

### 2.2.8. Pedone

This potential was designed to simulate a wide variety of crystalline and amorphous structures, including silicates [55]. It was adjusted to reproduce the equilibrium structure and elastic constants at room temperature of Cr₂O₃. It is a model of rigid ions with partial charges. The pair interactions are of the Morse type with additional short-range repulsive terms in r⁻¹²

$$E_{a\beta} = \frac{q_a q_{\beta}}{4\pi\epsilon_0 r_{a\beta}} + D_{AB} \left\{ e^{-2\rho_{AB}(r_{a\beta}-r_0)} - 2e^{-\rho_{AB}(r_{a\beta}-r_0)} \right\} + \frac{C_{AB}}{r_{a\beta}^{12}}$$

<table><caption>**Table 2**
: Ionic charges and parameters for the EAM contribution in the Owen potential.</caption><thead><tr><th>Ion</th><th>q (e)</th><th>G (eV·√Å)</th><th>n (Å<sup>5</sup>)</th></tr></thead><tbody><tr><td>O</td><td>-1.1104</td><td>0.690</td><td>106.85</td></tr><tr><td>Cr</td><td>1.6656</td><td>1.517</td><td>0</td></tr></tbody></table>
This repulsion allows for the simulation of phenomena involving rapid ion movements such as radiation damage, or high-temperature liquids required for glass preparation. The charges of chromium and oxygen ions are respectively +1.8 e and -1.2 e. The parameters of pair interactions are in Table 3. It does not appear to have been used in any studies on Cr₂O₃.

### 2.2.9. Sun

This potential was specifically designed to study structural defects in Cr₂O₃ and Al₂O₃ [56]. It is derived from the Matsui potential which was developed for the simulation of CaO—MgO—Al₂O₃—SiO₂ compounds [57]. It is a model of rigid ions with formal charges and Buckingham interactions, the parameters of which are in Table 3. This potential describes Cr₂O₃ with lattice parameters of 4.9872 Å and 13.5582 Å for a and c respectively. The formation energies of Frenkel pairs and Schottky defects using this potential [56] are given in Table 7. It was used in other studies on diffusion in Cr₂O₃ [29,30].

## 3. Reference data

### 3.1. Crystal structure

The crystal structure of Cr₂O₃ is of the corundum R̄3c type. This is a trigonal structure (prototype α-Al₂O₃, space group 167) that can be decomposed into two sub-lattices. The first, formed by the oxygen ions that occupy the Wyckoff 18e sites, has a slightly distorted hexagonal close-packed structure [58]. The second contains the chromium atoms that occupy the 12c sites, which correspond to 2/3 of the octahedral interstitial sites in the oxygen sub-lattice. The other octahedral sites (6b) are unoccupied in perfect Cr₂O₃, and are the preferred sites for chromium interstitials [39].

This structure is commonly described by a conventional hexagonal cell (green on Fig. 3) which contains 6 chemical formulas, or by a primitive rhombohedral cell (blue on the same figure) which contains 2. The chromium ions have octahedral coordination environments. The CrO₆ octahedra are distorted and contain three shorter Cr—O bonds of 1.96 Å and three longer bonds of 2.02 Å [59]. The Wyckoff positions corresponding to the conventional cell are given in Table 5.

Numerous values for the different structural parameters are available in the literature, resulting from either experiments or simulations. Some of these values are compiled in Table A.1 in the appendix. Experimental lattice parameters are consistent across studies, and the values of 4.96 Å and 13.59 Å can be retained for a and c, respectively. They are associated

<table><caption>**Table 4**

: Parameters for the polarisable ions in the empirical potentials considered in this study.</caption><thead><tr><th>Potential</th><th>Ion</th><th>Y (e)</th><th>K (eV/Å<sup>2</sup>)</th></tr></thead><tbody><tr><td>Atkinson [39]</td><td>O</td><td>-2.04</td><td>6.3</td></tr><tr><td>Grimes [40]</td><td>O</td><td>-2.04</td><td>6.3</td></tr><tr><td rowspan="2">Lawrence 1987 [51]</td><td>O</td><td>-2.18</td><td>27.29</td></tr><tr><td>Cr</td><td>0.97</td><td>100.0</td></tr><tr><td rowspan="2">Lawrence [48]</td><td>O</td><td>-2.21</td><td>16.00</td></tr><tr><td>Cr</td><td>0.97</td><td>999.9</td></tr><tr><td rowspan="2">Lewis [49]</td><td>O</td><td>-2.21</td><td>27.29</td></tr><tr><td>Cr</td><td>0.97</td><td>67.00</td></tr></tbody></table>

<table><caption>**Table 5**

: Crystallographic positions of atoms in a conventional unit cell of Cr₂O₃.</caption><thead><tr><th rowspan="2">Site</th><th colspan="3">Position</th></tr><tr><th>x</th><th>y</th><th>z</th></tr></thead><tbody><tr><td>O</td><td>18e</td><td>0</td><td>x(O)</td><td>3/4</td></tr><tr><td>Cr</td><td>12c</td><td>0</td><td>0</td><td>z(Cr)</td></tr><tr><td>Int.</td><td>6b</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

with standard deviations of 0.0119 Å (0.24 %) for a and 0.0261 Å (0.19 %) for c. These values will be used as references to determine the validity of empirical potentials. The magnetic transition has a minor effect on lattice parameters, with the difference between antiferromagnetic (AFM) and paramagnetic (PM) states being of the order of 0.005 Å (0.002, 0.003, and 0.013 are reported) for a and for c (0.0008, 0.005, and 0.02 are reported [60–62]). In all cases, this difference of approximately 0.2 % is negligible compared to the error introduced by empirical potentials, and thus the effect of magnetic structure on lattice parameters can be disregarded. The internal x and z parameters, which determine the positions of oxygen and chromium sites respectively, are also known to a good accuracy and presented in Table A.1.

Interatomic distances in Cr₂O₃, namely Cr—Cr, Cr—O, and O—O, are also well known, having been measured experimentally many times. The dispersion of the experimental values is between 0.21 % for the longer Cr—Cr(2) and 1 % for the shorter Cr—Cr(1). Some of the studies did not directly report bond lengths. In these cases, the lengths were calculated from the lattice parameters and atomic internal coordinates. The available literature values are listed in the appendix in Table A.2. Reference bond lengths obtained from the representative lattice

<caption>**Table 3**</caption>

: Parameters for the Buckingham and Morse functions of the potentials considered in this study. *The Pedone potential has interaction terms in r⁻¹² instead of the standard r⁻⁶ of Buckingham potentials.

<table><thead><tr><th>Potential</th><th>Pair</th><th>A (eV)</th><th>ρ (Å)</th><th>C (eV/Å<sup>6</sup>)*</th><th>D (eV)</th><th>β (Å<sup>-1</sup>)</th><th>r<sub>0</sub> (Å)</th></tr></thead><tbody><tr><td rowspan="3">Atkinson [39]</td><td>O—O</td><td>9547.96</td><td>0.21916</td><td>32.0</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cr—O</td><td>1313.18</td><td>0.3165</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>O—O</td><td>560.93434</td><td>0.36</td><td>4.2</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">Chartier [44]</td><td>Cr—O</td><td>698.77575</td><td>0.298944</td><td>—</td><td>1.3202651</td><td>2.0</td><td>1.8</td></tr><tr><td>Cr—Cr</td><td>1288.6597</td><td>0.118180</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">Grimes [40]</td><td>O—O</td><td>9547.96</td><td>0.2192</td><td>32.0</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cr—O</td><td>1204.18</td><td>0.3165</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">Lawrence 1987 [51]</td><td>O—O</td><td>22,764.3</td><td>0.149</td><td>27.88</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cr—O</td><td>1734.1</td><td>0.3010</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">Lawrence [48]</td><td>O—O</td><td>400.4</td><td>0.3994</td><td>27.88</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cr—O</td><td>1734.1</td><td>0.3010</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">Lewis [49]</td><td>O—O</td><td>22,764.3</td><td>0.149</td><td>112.2</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cr—O</td><td>1734.1</td><td>0.3010</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">Owen [53]</td><td>O—O</td><td>830.28</td><td>0.3529</td><td>3.884</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cr—O</td><td>1156.49</td><td>0.2795</td><td>—</td><td>1.348</td><td>2.245</td><td>1.799</td></tr><tr><td rowspan="2">Pedone [55]</td><td>O—O</td><td>2695.52</td><td>0.1179</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cr—O</td><td>—</td><td>—</td><td>22.0 eV.Å<sup>12</sup></td><td>0.042395</td><td>1.379316</td><td>3.618701</td></tr><tr><td rowspan="2">Sun [56]</td><td>O—O</td><td>91.55955</td><td>0.547148</td><td>1.0 eV.Å<sup>12</sup></td><td>0.399561</td><td>1.785079</td><td>2.340810</td></tr><tr><td>Cr—O</td><td>2362.905</td><td>0.270628</td><td>32.32</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>
<table><caption>**Table 6**
: Energy difference between the polymorphs of Cr₂O₃ and the stable structure R̄3c, in eV/atom; the C2/c structure is unstable.</caption><thead><tr><th rowspan="2">Source</th><th>C2/c</th><th>C2/m</th><th>Fdd2</th><th>Ia3̄</th><th>P3̄m1</th><th>Pbcn</th><th>Pccn</th><th>Pna21</th><th>Pnma</th><th>Pnma</th><th>R3̄</th></tr><tr><th>V2O3</th><th>Sm2O3</th><th>Au2O3</th><th>Bixbyite</th><th></th><th>Rh2O3-II</th><th>Sb2O3</th><th></th><th>Perovskite</th><th>Sb2S3</th><th></th></tr></thead><tbody><tr><td>[77]</td><td>0</td><td>0.109</td><td>0.158</td><td>0.034</td><td>0.179</td><td>0.032</td><td>0.363</td><td>0.116</td><td>0.157</td><td>0.237</td><td>-</td></tr><tr><td>[85]</td><td>-</td><td>0.099</td><td>0.156</td><td>0.029</td><td>0.246</td><td>-</td><td>0.466</td><td>-</td><td>0.175</td><td>-</td><td>-</td></tr><tr><td>[66]</td><td>-</td><td>-</td><td>-</td><td>0.047</td><td>-</td><td>0.038</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.06</td></tr></tbody></table>

<table><caption>**Table 7**
: Formation energies in eV of neutral defects in Cr₂O₃ according to the literature.</caption><thead><tr><th rowspan="2">Source</th><th colspan="2">Cr Frenkel</th><th rowspan="2">O Frenkel</th><th rowspan="2">Schottky</th><th rowspan="2">Notes</th></tr><tr><th>[0001]</th><th>Basal</th></tr></thead><tbody><tr><td>[23,28]</td><td>2.6</td><td>5.18</td><td>-</td><td>3.97</td><td>DFT (GGA + U)</td></tr><tr><td>[90]</td><td>4.56</td><td>-</td><td>-</td><td>9.24</td><td>DFT (HSE)</td></tr><tr><td>[90]</td><td>6.71</td><td>-</td><td>11.92</td><td>18.12</td><td>DFT (HSE), unbound</td></tr><tr><td>[27]</td><td>5.07</td><td>2.36</td><td>-</td><td>2.5</td><td>DFT (GGA + U)</td></tr><tr><td>[88]</td><td>2.33</td><td>-</td><td>6.75</td><td>10.72</td><td>DFT</td></tr><tr><td>[91]</td><td>12.72</td><td>-</td><td>12.38</td><td>21.1</td><td>Empirical pot.</td></tr><tr><td>[51]</td><td>13.28</td><td>-</td><td>5.84</td><td>20.75</td><td>Empirical pot.</td></tr><tr><td>[39]</td><td>15.6</td><td>-</td><td>10.68</td><td>27.95</td><td>Empirical pot.</td></tr><tr><td>[56]</td><td>8.46</td><td>-</td><td>12.06</td><td>11.25</td><td>Empirical pot.</td></tr></tbody></table>

<table><caption>**Table 8**
: Energy barriers in eV for the different diffusion mechanisms for chromium vacancies V''Cr in Cr₂O₃ from the literature; the lengths correspond to the distance between the initial and final sites in the perfect structure according to the Chartier potential.</caption><thead><tr><th rowspan="2"></th><th>Pathway 1</th><th>Pathway 2</th><th>Pathway 3</th><th>Pathway 4</th><th rowspan="2">Notes</th></tr><tr><th>(2.7 Å)</th><th>(2.9 Å)</th><th>(3.15 Å)</th><th>(3.59 Å)</th></tr></thead><tbody><tr><td>[91]</td><td>3.52</td><td>2.21</td><td>2.11</td><td>1.75</td><td>Empirical pot.</td></tr><tr><td>[93]</td><td>2.02</td><td>-</td><td>-</td><td>-</td><td>GGA</td></tr><tr><td>[27]</td><td>3.21</td><td>2.57</td><td>-</td><td>2.73</td><td>GGA + U</td></tr><tr><td>[26]</td><td>2.09</td><td>1.91</td><td>2.08</td><td>2.14</td><td>GGA</td></tr><tr><td>[23]</td><td>2.56</td><td>2.02</td><td>2.18</td><td>2.59</td><td>GGA + U</td></tr><tr><td>[25]</td><td>2.27</td><td>2.01</td><td>-</td><td>2.38</td><td>GGA + U</td></tr><tr><td></td><td>2.67</td><td>0.971</td><td>8.19</td><td>1.23</td><td>Empirical pot.</td></tr></tbody></table>

<table><caption>**Table 9**
: Energy barriers in eV for the different diffusion mechanisms for oxygen vacancies V''O in Cr₂O₃ from the literature; the lengths correspond to the distance between the initial and final sites in the perfect structure according to the Chartier potential.</caption><thead><tr><th rowspan="2"></th><th>Pathway 1</th><th>Pathway 2</th><th>Pathway 3</th><th>Pathway 4</th><th rowspan="2">Notes</th></tr><tr><th>(2.59 Å)</th><th>(2.70 Å)</th><th>(2.82 Å)</th><th>(2.99 Å)</th></tr></thead><tbody><tr><td>[91]</td><td>1.12</td><td>2.49</td><td>2.54</td><td>-</td><td>Empirical pot.</td></tr><tr><td>[93]</td><td>1.15</td><td>-</td><td>-</td><td>-</td><td>GGA</td></tr><tr><td>[27]</td><td>2.21</td><td>2.35</td><td>3.22</td><td>3.65</td><td>GGA + U</td></tr><tr><td>[26]</td><td>1.32</td><td>1.69</td><td>2.36</td><td>2.20</td><td>GGA</td></tr><tr><td>[23]</td><td>1.18</td><td>1.52</td><td>2.62</td><td>2.98</td><td>GGA + U</td></tr><tr><td>[25]</td><td>1.22</td><td>1.46</td><td>2.68</td><td>2.52</td><td>GGA + U</td></tr><tr><td></td><td>0.425</td><td>1.38</td><td>1.63</td><td>2.38</td><td>Empirical pot.</td></tr></tbody></table>

parameters and atomic positions are:

- 2.64 Å and 2.89 Å for Cr—Cr bonds, respectively, for pairs aligned along the [0001] axis and within the basal plane [63,64],

- 1.97 Å and 2.01 Å for Cr—O bonds,

- 2.63 Å and 2.74 Å for O—O bonds.

Most DFT calculations in the literature use the DFT + U formalism to correct the behaviour of chromium's 3d electrons. The U repulsion parameter is typically adjusted to reproduce a combination of crucial values, such as lattice parameters, the magnetic moment of chromium,

<table><caption>**Table 10**
: Surface energies of symmetrical tilt grain boundaries in Cr₂O₃ according to the literature.</caption><thead><tr><th>Surface</th><th>Energy (J/m²)</th><th>Source</th></tr></thead><tbody><tr><td>(1010) prismatic</td><td>0.59</td><td>[50]</td></tr><tr><td>(0001) basal</td><td>1.24</td><td>[50]</td></tr><tr><td>(1011)</td><td>1.56</td><td>[50]</td></tr><tr><td>(1012)</td><td>2.72</td><td>[50]</td></tr><tr><td>(1120)</td><td>3.88</td><td>[50]</td></tr><tr><td>(1010) prS</td><td>0.23</td><td>[96]</td></tr><tr><td>(1012) rS(V)</td><td>0.33</td><td>[96]</td></tr><tr><td>(1010) prG</td><td>0.40</td><td>[96]</td></tr><tr><td>(0001) bR</td><td>0.61</td><td>[96]</td></tr><tr><td>(1012) rG(O)</td><td>1.08</td><td>[96]</td></tr><tr><td>(1014) pyG(Cr)</td><td>1.14</td><td>[96]</td></tr><tr><td>(1014) pyG(O)</td><td>1.14</td><td>[96]</td></tr><tr><td>(0001) bM</td><td>1.35</td><td>[96]</td></tr><tr><td>(0001) bG</td><td>1.59</td><td>[96]</td></tr><tr><td>(1012) rG(V)</td><td>1.93</td><td>[96]</td></tr></tbody></table>

<table><caption>**Table 11**
: Lattice parameters and linear thermal expansion coefficients of Cr₂O₃ at room temperature from MD simulations; the values in brackets are the relative error compared to the reference values in %.</caption><thead><tr><th></th><th>a (Å)</th><th>c (Å)</th><th>α<sub>a</sub> (10<sup>-6</sup> K<sup>-1</sup>)</th><th>α<sub>c</sub> (10<sup>-6</sup> K<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Reference</td><td>4.96</td><td>13.59</td><td>7.1</td><td>5.1</td></tr><tr><td>Atkinson</td><td>5.05</td><td>13.35</td><td>5.7</td><td>5.9</td></tr><tr><td></td><td>(1.81)</td><td>(-1.76)</td><td>(-19.71)</td><td>(15.69)</td></tr><tr><td>Chartier</td><td>4.96</td><td>13.42</td><td>7.5</td><td>6.8</td></tr><tr><td></td><td>(0.00)</td><td>(-1.25)</td><td>(5.63)</td><td>(33.33)</td></tr><tr><td>Lawrence</td><td>5.22</td><td>14.13</td><td>6.74</td><td>7.48</td></tr><tr><td></td><td>(5.24)</td><td>(3.97)</td><td>(-5.07)</td><td>(46.67)</td></tr><tr><td>Lewis</td><td>4.91</td><td>12.49</td><td>4.02</td><td>1.10</td></tr><tr><td></td><td>(-1.00)</td><td>(-8.09)</td><td>(-43.38)</td><td>(-78.43)</td></tr><tr><td>Owen</td><td>4.96</td><td>13.59</td><td>8.2</td><td>7.8</td></tr><tr><td></td><td>(0.00)</td><td>(0.00)</td><td>(15.49)</td><td>(52.94)</td></tr><tr><td>Pedone</td><td>5.01</td><td>13.62</td><td>7.6</td><td>6.5</td></tr><tr><td></td><td>(1.00)</td><td>(0.22)</td><td>(7.04)</td><td>(27.45)</td></tr><tr><td>Sun</td><td>5.00</td><td>13.59</td><td>8.33</td><td>9.69</td></tr><tr><td></td><td>(0.80)</td><td>(0.00)</td><td>(17.32)</td><td>(90)</td></tr></tbody></table>

or the bandgap width of the material's electronic density of states. The lattice parameters used for this adjustment are often experimental values measured at 300 K, even in the case of DFT calculations at 0 K. Consequently, simulated lattice parameters are generally less reliable than experimental values and will not be used for empirical potential validation.

## 3.2. Elastic properties

Numerous measurements of the bulk modulus are available, with some gathered in Table A.2. However, individual elastic constants C<sub>ij</sub> are less common. There is only one complete experimental determination of the elastic tensor obtained through acoustic measurements on a Cr₂O₃ monocrystal at 20 °C [60]. These elastic constants are listed in the appendix in Table A.3. This dataset is limited, with only one experimental reference, Alberts & Boeyen [60], and one DFT study using the quasi-harmonic approximation (QHA), Wang et al. [65]. This techniques
<caption>**Table 12**</caption>

Elastic constants of Cr₂O₃ at room temperature calculated using structure parameters from MD simulations; the values in brackets are the relative difference compared to the reference values in %.

<table><caption>**Table 13**</caption><thead><tr><th></th><th><i>B</i><br>(GPa)</th><th><i>C</i><sub>11</sub><br>(GPa)</th><th><i>C</i><sub>12</sub><br>(GPa)</th><th><i>C</i><sub>13</sub><br>(GPa)</th><th><i>C</i><sub>14</sub><br>(GPa)</th><th><i>C</i><sub>33</sub><br>(GPa)</th><th><i>C</i><sub>44</sub><br>(GPa)</th><th>Average difference (%)</th></tr></thead><tbody><tr><td>Reference</td><td>222</td><td>374</td><td>148</td><td>175</td><td>-19</td><td>362</td><td>159</td><td>-</td></tr><tr><td>Atkinson</td><td>295<br>(31)</td><td>564<br>(51)</td><td>222<br>(50)</td><td>165<br>(-5)</td><td>47<br>(346)</td><td>452<br>(25)</td><td>137<br>(-14)</td><td>92</td></tr><tr><td>Chartier</td><td>235<br>(6)</td><td>391<br>(5)</td><td>174<br>(18)</td><td>165<br>(-5)</td><td>2<br>(111)</td><td>333<br>(-8)</td><td>132<br>(-17)</td><td>30</td></tr><tr><td>Lawrence</td><td>275<br>(24)</td><td>513<br>(37)</td><td>217<br>(47)</td><td>138<br>(-21)</td><td>52<br>(372)</td><td>474<br>(31)</td><td>90<br>(-43)</td><td>99</td></tr><tr><td>Lewis</td><td>310<br>(40)</td><td>589<br>(57)</td><td>281<br>(90)</td><td>281<br>(61)</td><td>24<br>(229)</td><td>361<br>(0)</td><td>147<br>(-8)</td><td>68</td></tr><tr><td>Owen</td><td>262<br>(18)</td><td>418<br>(12)</td><td>195<br>(31)</td><td>188<br>(7)</td><td>3<br>(115)</td><td>385<br>(6)</td><td>129<br>(-19)</td><td>32</td></tr><tr><td>Pedone</td><td>233<br>(5)</td><td>387<br>(4)</td><td>164<br>(11)</td><td>161<br>(-8)</td><td>6<br>(131)</td><td>348<br>(-4)</td><td>126<br>(-21)</td><td>34</td></tr><tr><td>Sun</td><td>375<br>(69)</td><td>683<br>(83)</td><td>294<br>(99)</td><td>204<br>(17)</td><td>62<br>(428)</td><td>620<br>(71)</td><td>127<br>(-20)</td><td>118</td></tr></tbody></table>

Bond lengths in Cr₂O₃ at room temperature obtained after energy minimisations at the lattice parameters from MD simulations; the values in brackets are the relative difference compared to the reference values in %.

<table><caption>**Table 14**</caption><thead><tr><th></th><th>Cr-Cr (1) (Å)</th><th>Cr-Cr (2) (Å)</th><th>Cr-O (1) (Å)</th><th>Cr-O (2) (Å)</th><th>O-O (1) (Å)</th><th>O-O (2) (Å)</th><th>Average difference (%)</th></tr></thead><tbody><tr><td>Reference</td><td>2.65</td><td>2.88</td><td>1.97</td><td>2.02</td><td>2.63</td><td>2.74</td><td>-</td></tr><tr><td>Atkinson</td><td>2.90<br>(9.43)</td><td>3.00<br>(4.16)</td><td>1.94<br>(-1.52)</td><td>2.10<br>(3.96)</td><td>2.63<br>(0.00)</td><td>2.70<br>(-1.46)</td><td>2.93</td></tr><tr><td>Chartier</td><td>2.71<br>(2.26)</td><td>2.94<br>(2.08)</td><td>1.94<br>(-1.52)</td><td>2.03<br>(0.49)</td><td>2.62<br>(-0.38)</td><td>2.68<br>(-2.19)</td><td>1.27</td></tr><tr><td>Lawrence</td><td>3.05<br>(15.09)</td><td>3.09<br>(7.29)</td><td>2.01<br>(2.03)</td><td>2.18<br>(7.92)</td><td>2.70<br>(2.66)</td><td>2.84<br>(3.65)</td><td>5.52</td></tr><tr><td>Lewis</td><td>2.77<br>(4.53)</td><td>2.91<br>(1.04)</td><td>1.89<br>(-4.06)</td><td>1.97<br>(-2.48)</td><td>2.44<br>(-7.22)</td><td>2.54<br>(-7.30)</td><td>3.80</td></tr><tr><td>Owen</td><td>2.68<br>(1.13)</td><td>2.95<br>(2.43)</td><td>1.95<br>(-1.02)</td><td>2.038<br>(0.89)</td><td>2.66<br>(1.14)</td><td>2.69<br>(-1.82)</td><td>1.20</td></tr><tr><td>Pedone</td><td>2.76<br>(4.15)</td><td>2.98<br>(3.47)</td><td>1.95<br>(-1.02)</td><td>2.08<br>(2.97)</td><td>2.69<br>(2.28)</td><td>2.71<br>(-1.09)</td><td>2.14</td></tr><tr><td>Sun</td><td>2.86<br>(7.92)</td><td>2.99<br>(3.82)</td><td>1.94<br>(-1.52)</td><td>2.094<br>(3.66)</td><td>2.65<br>(0.76)</td><td>2.70<br>(-1.46)</td><td>2.74</td></tr></tbody></table>

Thermodynamic properties of Cr₂O₃ at room temperature calculated using the empirical potentials. The values for Cₚ³₀ were obtained from MD simulations and those for CₚV³₀ from QHA calculations. The values in brackets are the relative difference compared to the reference values in %.

<table><thead><tr><th>Potential</th><th>C<sub>p</sub><sup>300</sup> (10<sup>-4</sup> eV/atom<br>K)</th><th>C<sub>v</sub><sup>300</sup> (10<sup>-4</sup> eV/atom<br>K)</th><th>S<sup>300</sup> (10<sup>-4</sup> eV/atom<br>K)</th></tr></thead><tbody><tr><td>Reference</td><td>2.172 [73]</td><td>1.7961 [65]</td><td>1.665 [74]</td></tr><tr><td>Atkinson</td><td>2.5318<br>(16.57)</td><td>1.698<br>(-5.46)</td><td>1.269<br>(-23.78)</td></tr><tr><td>Chartier</td><td>2.6336<br>(21.25)</td><td>1.779<br>(-0.95)</td><td>1.410<br>(-15.31)</td></tr><tr><td>Lawrence</td><td>2.5941<br>(19.43)</td><td>2.071<br>(15.30)</td><td>1.781<br>(6.97)</td></tr><tr><td>Lewis</td><td>2.6154<br>(20.41)</td><td>1.420<br>(-20.94)</td><td>1.186<br>(-28.77)</td></tr><tr><td>Owen</td><td>2.6180<br>(20.53)</td><td>1.670<br>(-7.02)</td><td>1.288<br>(-22.64)</td></tr><tr><td>Pedone</td><td>2.6392<br>(21.51)</td><td>1.774<br>(-1.23)</td><td>1.429<br>(-14.17)</td></tr><tr><td>Sun</td><td>2.7465<br>(26.45)</td><td>1.654<br>(-7.91)</td><td>1.221<br>(-26.66)</td></tr></tbody></table>

approximates temperature effects, therefore its results are directly comparable to experimental measurements at ambient temperature. The difference between DFT and experimental values can be quantified by computing the mean of the absolute value of relative differences for each elastic constant

<caption>**Table 15**</caption>

Formation energy in eV of chromium defects in Cr₂O₃ using empirical potentials; the Frenkel energies with a star correspond to the structure with both defects in the same basal plane.

<table><caption>**Table 16**</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Cr<sup>***</sup> Octa.</th><th rowspan="2">V<sup>**</sup><sub>Cr</sub></th><th rowspan="2">Triple defect</th><th colspan="2">Frenkel</th></tr><tr><th>Far</th><th>Close</th></tr></thead><tbody><tr><td>Reference</td><td>-</td><td>-</td><td>-</td><td>9.63-10.05</td><td>2.33-2.6</td></tr><tr><td>Atkinson</td><td>4.95</td><td>7.89</td><td>7.81</td><td>12.81</td><td>10.49*</td></tr><tr><td>Chartier</td><td>-1.21</td><td>7.54</td><td>7.54</td><td>6.32</td><td>2.08</td></tr><tr><td>Lawrence</td><td>1.70</td><td>8.05</td><td>7.78</td><td>9.74</td><td>8.08*</td></tr><tr><td>Owen</td><td>-1.37</td><td>5.26</td><td>5.18</td><td>3.74</td><td>1.14</td></tr><tr><td>Pedone</td><td>-0.51</td><td>8.03</td><td>7.76</td><td>7.53</td><td>5.65*</td></tr><tr><td>Sun</td><td>-0.12</td><td>7.42</td><td>6.75</td><td>7.33</td><td>1.53</td></tr></tbody></table>

Formation energy in eV of oxygen defects in Cr₂O₃ using the empirical potentials.

<table><thead><tr><th rowspan="2"></th><th rowspan="2">O<sub>i</sub><sup>o</sup> octa.</th><th rowspan="2">V<sup>**</sup><sub>O</sub></th><th rowspan="2">Doublet</th><th colspan="2">Frenkel</th></tr><tr><th>Far</th><th>Close</th></tr></thead><tbody><tr><td>Reference</td><td>-</td><td>-</td><td>-</td><td>7.71-11.85</td><td>6.75</td></tr><tr><td>Atkinson</td><td>8.03</td><td>2.49</td><td>7.09</td><td>9.59</td><td>-</td></tr><tr><td>Chartier</td><td>7.21</td><td>1.64</td><td>6.73</td><td>8.37</td><td>-</td></tr><tr><td>Lawrence</td><td>11.42</td><td>0.48</td><td>10.47</td><td>9.74</td><td>-</td></tr><tr><td>Owen</td><td>6.24</td><td>1.38</td><td>6.12</td><td>7.50</td><td>-</td></tr><tr><td>Pedone</td><td>7.64</td><td>2.59</td><td>7.22</td><td>9.81</td><td>-</td></tr><tr><td>Sun</td><td>11.63</td><td>-1.71</td><td>11.69</td><td>9.68</td><td>9.28</td></tr></tbody></table>
<table><caption>**Table 17**
: Formation energy in eV of Schottky defects in Cr₂O₃ using empirical potentials.</caption><thead><tr><th></th><th>Far</th><th>Vertical</th><th>Most stable</th></tr></thead><tbody><tr><td>Reference</td><td>19.92</td><td>-</td><td>3.97-10.72</td></tr><tr><td>Atkinson</td><td>23.23</td><td>13.06</td><td>10.61</td></tr><tr><td>Barrera</td><td>19.25</td><td>-</td><td>7.95</td></tr><tr><td>Chartier</td><td>20.20</td><td>10.69</td><td>10.04</td></tr><tr><td>Lawrence</td><td>17.03</td><td>10.11</td><td>7.34</td></tr><tr><td>Owen</td><td>14.71</td><td>6.66</td><td>6.66</td></tr><tr><td>Pedone</td><td>24.07</td><td>13.43</td><td>12.64</td></tr><tr><td>Sun</td><td>17.10</td><td>10.57</td><td>6.48</td></tr></tbody></table>

<table><caption>**Table 18**
: Energy barriers in eV of the different diffusion pathways for V<sub>Cr</sub><sup>cr</sup> vacancies in Cr₂O₃.</caption><thead><tr><th>Potential</th><th>Pathway 1</th><th>Pathway 2</th><th>Pathway 3</th><th>Pathway 4</th></tr></thead><tbody><tr><td>Reference</td><td>2.56</td><td>2.02</td><td>2.18</td><td>2.59</td></tr><tr><td>Chartier</td><td>2.45</td><td>2.23</td><td>2.02</td><td>1.49</td></tr><tr><td>Owen</td><td>1.23</td><td>1.20</td><td>1.32</td><td>1.09</td></tr><tr><td>Pedone</td><td>2.14</td><td>1.86</td><td>1.77</td><td>1.17</td></tr></tbody></table>

<table><caption>**Table 19**
: Energy barriers in eV of the different diffusion pathways for V<sub>O</sub><sup>cr</sup> vacancies in Cr₂O₃.</caption><thead><tr><th>Potential</th><th>Pathway 1</th><th>Pathway 2</th><th>Pathway 3</th><th>Pathway 4</th></tr></thead><tbody><tr><td>Reference</td><td>1.18</td><td>1.52</td><td>2.68</td><td>2.52</td></tr><tr><td>Chartier</td><td>0.85</td><td>1.26</td><td>2.26</td><td>2.83</td></tr><tr><td>Owen</td><td>0.49</td><td>1.13</td><td>1.67</td><td>2.51</td></tr><tr><td>Pedone</td><td>1.10</td><td>2.36</td><td>1.29</td><td>2.67</td></tr></tbody></table>

<table><caption>**Table 20**
: Surface energy in J/m² of some symmetrical tilt grain boundaries in Cr₂O₃ predicted by empirical potentials. The values in brackets are relative differences in % with reference values [96].</caption><thead><tr><th>Grain boundary</th><th>Reference</th><th>Chartier</th><th>Owen</th><th>Pedone</th></tr></thead><tbody><tr><td>(1010) prS</td><td>0.23</td><td>0.32 (34)</td><td>0.28 (18)</td><td>0.30 (25)</td></tr><tr><td>(1012) rS(V)</td><td>0.33</td><td>0.21 (-37)</td><td>0.58 (75)</td><td>0.411 (25)</td></tr><tr><td>(1010̄) prG</td><td>0.40</td><td>0.53 (33)</td><td>0.45 (13)</td><td>0.73 (82)</td></tr><tr><td>(0001) bR</td><td>0.61</td><td>0.70 (15)</td><td>1.57 (157)</td><td>1.85 (203)</td></tr><tr><td>(1012) rG(O)</td><td>1.08</td><td>0.21 (-81)</td><td>0.58 (46)</td><td>0.41 (62)</td></tr><tr><td>(101̄4) pyG(Cr)</td><td>1.14</td><td>2.43 (113)</td><td>2.80 (145)</td><td>2.80 (146)</td></tr><tr><td>(1014) pyG(O)</td><td>1.14</td><td>1.29 (13)</td><td>1.54 (35)</td><td>1.79 (57)</td></tr><tr><td>(0001) bM</td><td>1.35</td><td>1.46 (8)</td><td>1.00 (26)</td><td>1.98 (47)</td></tr><tr><td>(0001) bG</td><td>1.59</td><td>0.91 (-43)</td><td>0.59 (63)</td><td>0.63 (60)</td></tr><tr><td>(1012) rG(V)</td><td>1.93</td><td>2.28 (18)</td><td>2.65 (37)</td><td>2.93 (52)</td></tr><tr><td>Average difference (%)</td><td>-</td><td>40</td><td>62</td><td>76</td></tr></tbody></table>

$$ \Delta = \sum_{ij} \left| \frac{C_{ij}^{\text{DFT}} - C_{ij}^{\text{exp}}}{C_{ij}^{\text{exp}}} \right| $$

This average difference between the QHA calculations [65] and the experimental reference [60] is 4 %. Lastly, elastic constants obtained from automated DFT calculations are available in the Materials Project database [66]. However, these values are derived at 0 K without correction to reproduce finite temperature. They are less consistent with experimental data than Wang et al.'s calculations, with an average absolute difference of 20 %. Experimental values at ambient temperature [60] serve as references for evaluating the quality of empirical potentials.

### 3.3. Thermophysical properties

Several studies provide measurements of thermodynamic properties of Cr₂O₃, such as heat capacity, shown in Fig. 4, enthalpy increment, or

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Crystal structure of Cr₂O₃, showing the hexagonal conventional unit cell in green and primitive rhombohedral cell in blue. Chromium ions appear as grey balls and oxygen ions as smaller red balls.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Constant-pressure heat capacity of Cr₂O₃ as a function of temperature according to the literature.</figcaption></figure>

entropy. An important feature of the heat capacity is the sharp peak associated with the AFM → PM [67–69]. This peak is clearly defined and allows for precise determination of the Néel temperature T_N at 305.5 K [70]. However, because the empirical potentials considered here cannot reproduce magnetic effects, this peak must be ignored to validate the potentials. Some data is also available for the PM state beyond T_N [70–74], which will be the bases for the potentials assessment.

The thermodynamic functions from the different authors generally agree well within the temperature range of interest (from 300 K to 1600 K), except those of Ziemniak et al. and Wang et al. [65,74] which underestimate C<sub>P</sub>. In the case of quasi-harmonic calculations [65],
systematic underestimation can be attributed to anharmonic effects such as phonon interactions or defects present in physical samples that cannot be represented in the model. The expression by Klemme et al. reaches a maximum around 1300 K, which appears unrealistic. However, it is in very good agreement with the expressions of Gurevich et al. [72], Moore and Kelley [73], and the JANAF tables [71] below 1200 K.

Selecting reference values for heat capacity among the literature values presented in Table A.4 is challenging due to the impact of magnetic transition on experimental results at 300 K, which do not correspond to the PM state as seen in molecular dynamics simulations. The experimental value $C_p^{300}$ is thus overestimated. On the other hand, results from quasi-harmonic DFT calculation are likely underestimated. Therefore, the selected reference value is from the Moore and Kelley polynomial expression with adjusted coefficients fitted for the PM state behaviour at high temperature [73], extrapolated to 300 K.

Several values for the thermal expansion coefficients are available in the literature, summarised in Table A.5. They are typically based on measurements of length changes between two temperatures, from which an approximate coefficient of thermal expansion can be calculated using the relation $\Delta L/L = \alpha_L \times \Delta T$. However, this relation tends to overestimate $\alpha_L$ when $\Delta L/L$ is non-linear. There are two studies demonstrating the temperature dependence of thermal expansion coefficients: one experimental by Huntz et al. [75] and one based on quasi-harmonic DFT calculations by Wang et al. [65]. In all cases, there are no high-temperature data available (beyond 1200 K for Huntz and 1500 K for Wang). Both studies provide the same thermal expansion coefficient at 500 K, but with significantly different temperature trends. The former exhibits a linear trend with a small slope, while the latter varies more rapidly.

Thermal expansion in $Cr_2O_3$ is anisotropic, with $\alpha_a > \alpha_b$ for all sources that calculate these two coefficients separately. A quantitative discussion of the anisotropic character of thermal expansion is challenging due to the scatter in experimental results. However, a good potential is expected to reproduce this behaviour qualitatively.

### 3.4. Polymorphs

A phase transition to a *Pbcn* structure (prototype $Rh_2O_3$-II) is predicted under pressure by DFT calculations [76–78]. The calculated pressure for this transition is 30 GPa [78] for the AFM state and around 14 GPa and 16 GPa for the PM state [76]. An experimental observation of this transition was reported at a similar pressure of 13 GPa [79], marked by changes in the Raman spectrum of the material. However, high-pressure X-ray diffraction analysis did not provide evidence of the structural change, which was attributed to the similarity between the diffraction patterns of the *R*3c and *Pbcn* structures. Another study [80] suggested that the structure appearing at 15 GPa is the monoclinic *I*2/*a* structure (space group *I*5, equivalent to *C*2/*c* ) and corresponding to the *V*2*O*3 prototype [81]. The *R*3c → *I*2/*a* phase transition is described as displacive and consisting of a shift in chromium ions from their off-centred position in *R*3c to the centre of their octahedra in *I*2/*a*. This change is visible in Raman spectra, with the appearance of two bands not present in *R*3c, but not in X-ray diffraction due to the similarity of the two structures. The same effect explains why no DFT calculation was able to simulate the *I*2/*a* structure, the crystal reverting to its stable state during geometry optimisation [66,77]. This is consistent with other experiments where no structural change was observed [82–84], but which showed an anomaly at 20 GPa resulting from magnetic structure changes accompanied by a change in the trend of the *c/a* ratio [84]. Furthermore, even if this transition was not described in $Cr_2O_3$, it is documented in the case of $V_2O_3$, for which the *I*2/*a* structure is stable below 155 K but *R*3c is stable at room temperature [81].

Another transition to an orthorhombic structure was observed experimentally at higher pressures, around 30 GPa [80], consistent with initial DFT estimations [78]. However, the exact structure of this

orthorhombic phase remains unclear, as it may be *Pbcn* or a *Pnma* perovskite structure (prototype $GdFeO_3$). Simulations of the *Pnma* structure in DFT showed it to be mechanically unstable, with imaginary vibrational modes after geometry optimisation [77]. These DFT calculations suggest a spontaneous transformation of *Pnma* into *Pbcn*, which is only slightly energetically unfavourable compared to *R*3c. This supports the *R*3c→*Pbcn* transition over the *Pnma*→*Pnma* transition documented in the literature. The high-pressure phase transition is reconstructive and requires some energy, for example provided by laser heating [80]. The *I*2/*a* structure remains metastable up to at least 60 GPa [80].

Other DFT calculations also indicate that an *Ia*3 structure (bixbyte prototype) is nearly energetically degenerate with *Pbcn* and more favourable than *Pnma* [77]. The energy differences between these polymorphs and the ground state *R*3c structure are presented in Table 6.

Numerous structures for $Cr_2O_3$ are available in databases such as the Materials Project [66] and OQMD [85], with DFT calculations allowing for an assessment of the predictive capability of empirical potentials. However, these calculations were performed automatically and do not offer the same level of reliability as found in the literature. For instance, according to the Materials Project, the ground state of $Cr_2O_3$ is ferromagnetic (FM), which contradicts the entire body of literature that describes it as antiferromagnetic (AFM). Similarly, OQMD identifies the *C*2/*c* structure (equivalent to *I*2/*a*) as the ground state, with a 0.004 eV/atom difference compared to *R*3c. Therefore, these database calculations serve more as qualitative indicators rather than quantitative criteria for estimating potential quality and should be used with caution. The reference energies employed for potential validation are those from Wessel et al., except for the *R*3 structure, where the energy chosen is from the Materials Project database.

A phonon spectrum for the *R*3c structure was reported in the literature, derived from DFT calculations [65], which provides a qualitative point of comparison. This spectrum features notable characteristics such as a peak at 7 THz and a band gap between 14.2 THz and 16.2 THz. However, this spectrum alone is insufficient for rigorous validation of spectra obtained from empirical potentials.

### 3.5. Point defects

The dominant defects at high temperatures (beyond 1100 °C) are intrinsic defects [86]. Chromia exhibits *n*-type semiconductor behaviour at low $p_{O_2}$, which favours chromium interstitials or oxygen vacancies, and *p*-type behaviour when $p_{O_2}$ is near ambient conditions, promoting chromium vacancies and oxygen interstitials [86]. Experimental evidence suggests that the dominant defect is the chromium interstitial $Cr_i^{***}$ for $p_{O_2} < 10^{-6}$ atm in the *n*-type regime [87], and $V_{Cr}^{'''}$ for $p_{O_2}$ around 1 atm in the *p*-type regime [8].

The formation energies of charged defects depend on the chemical potential of electrons and are relative to reference states (generally metallic chromium and molecular oxygen) that cannot be simulated using empirical potentials. Therefore, DFT results on charged defects cannot be used for validating empirical potentials. However, simulations of isolated point defects can be used to verify the most favourable structure for each charged defect (e.g., tetrahedral or octahedral sites for interstitials). The formation energies of electrically neutral defect clusters, namely Frenkel pairs and Schottky defects, can also be directly compared between simulations based on the DFT and empirical potentials.

#### 3.5.1. Chromium defects

Chromium interstitials occupy 6b octahedral sites of the *R*3c structure preferentially to tetrahedral sites shown in Fig. 5. The coordination environment of these interstitials in Fig. 6a is octahedral and similar to that of chromium ions in the perfect crystal. All cation sites in the perfect structure are equivalent by symmetry, so there is only one possibility for
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Successive O-Cr-O basal planes showing an octahedral (grey cube) and a tetrahedral (red cube) interstitial sites. Chromium ions appear as grey balls and oxygen ions as smaller red balls.</figcaption></figure>

chromium vacancies. However, depending on the charge state of the defect, a $V_{Cr}^{a}$ vacancy may be less favourable than a pair of vacancies surrounding an interstitial aligned along the c-axis, in the configuration shown in Fig. 6b [23]. This structure is known as a triple defect or split vacancy and represents an intermediate state involved in one of the chromium vacancy diffusion mechanisms [23,27].

The interstitial in a chromium Frenkel pair is always located on an octahedral site. The vacancy can be situated on any of the neighbouring sites, numbered 1 to 5 in Fig. 6c. Calculations in the literature consider only configurations 1 and 2 [27]. These cases correspond to instances where: 1) the vacancy and interstitial are aligned along the [0001] direction with the interstitial as a nearest neighbour of the vacancy; 2) the vacancy and interstitial are in the same basal plane (0001), with the interstitial as a second-nearest neighbour.

Results from DFT calculations in the literature are compiled in Table 7. The formation energies for the most stable Cr Frenkel pair are consistent, between 2.3 eV and 2.6 eV. However, there is some disagreement regarding the structure, one source [27] finding that the most favourable is in the basal plane (configuration 2 in Fig. 6c), and others [23,88] showing that the formation energy is lower if the defects

are aligned along the [0001] axis (configuration 1 in Fig. 6c). Defect simulations based on empirical potentials and the Mott-Littleton method [89] produced higher formation energies, between 8 eV and 16 eV depending on the potential used. The specific configuration considered was not specified in these studies.

### 3.5.2. Oxygen defects

The most favourable site for oxygen interstitials is the octahedral site (Fig. 5), with the interstitial displaced from the centre of the octahedron along the [0001] axis [88,92]. A peroxide structure was also observed, which is a pair of oxygen ions occupying a single oxygen site. In this configuration, the involved oxygen nuclei relax towards adjacent octahedral interstitial sites, but remain close with a separation distance of 1.48 Å [88]. This structure is more favourable than the interstitial on the octahedral site [28,88,92], with a difference in energy of approximately 1 eV [88]. The preferred Frenkel pair for oxygen is the structure with the interstitial at an octahedral site adjacent to the vacancy [88]. The available values from the literature, including those obtained using empirical potentials, are in Table 7.

### 3.5.3. Schottky defects

Schottky defects in $Cr_2O_3$ consist of two chromium vacancies and three oxygen vacancies. One study [23] find the first structure to be more favourable (3.98 eV versus 4.77 eV). The other finds that the most favourable Schottky defect has the two $V_{Cr}$ as second neighbours, in the same (0001) plane, and the oxygen vacancies distributed between the O layers immediately above and below the chromium vacancies [88]. The other studies do not discuss the structure of Schottky clusters. The formation energies of Schottky defects available in the literature are presented in Table 7, which includes results from DFT simulations and empirical potential calculations.

## 3.6. Diffusion pathways

Several studies on the diffusion pathways of point defects are available in the literature. These are atomistic simulations based on the DFT [23,26,27,93] or empirical potentials [25,91] and use the nudged

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Chromium defects in $Cr_2O_3$: a) chromium interstitial on a octahedral site; b) chromium triple defect; c) possible structures for chromium Frenkel pairs. Chromium Frenkel pairs can be: 1) aligned along [0001]; 2) with the vacancy as a second neighbour of the interstitial in a basal plane (0001); 3) with the vacancy as a third neighbour; 4) with the vacancy as a fourth neighbour along [0001]; 5) with the vacancy as a fifth neighbour. Chromium ions appear as grey balls, chromium interstitials as blue balls, and oxygen ions as smaller red balls. Chromium vacancies appear as grey cubes.</figcaption></figure>
elastic band (NEB) method [37].

Four pathways are generally considered for the diffusion of chromium vacancies, corresponding to the four nearest neighbour spheres around a chromium ion (Fig. 7a):

1. Between two nearest neighbour sites aligned along the [0001] axis, where the coordination polyhedra share a face;

2. Between two nearest neighbour sites in the basal plane (0001);

3. Between two third-nearest neighbour sites in two adjacent basal planes;

4. Between two fourth-nearest neighbour sites along the [0001] axis, through an interstitial site.

Fig. 8 shows examples of NEB calculation results, taken from a study by Medasani et al. [23]. According to the literature, the most favourable pathway is typically within the basal plane [23,26] or along the [0001] axis [27]. A study on the effect of chromium vacancy charge states on diffusion pathways has shown that charged vacancies $V''_{Cr}$ have lower energy barriers than other charge states [23]. The vacancy pair in Fig. 6b is an intermediate state in the diffusion mechanism through pathway 4. According to calculations, this state is either metastable [27] or slightly more stable than an isolated vacancy [23,26].

All DFT calculations in Table 8 indicate that pathway 2 is the most favourable, with energy barriers between 1.91 eV and 2.57 eV. However, the migration energies for all mechanisms are relatively similar and below 3.21 eV. The variations between authors are greater than the differences between paths for a single author. This dispersion is also of the same order as the formation energies of defects in DFT from the literature [94]. This makes it impossible to use this data to validate the potentials quantitatively. Works in the literature using empirical potentials [25,91] provide reasonable migration energies, but find that the most favourable mechanism is pathway 4.

Oxygen vacancies also have four primary diffusion pathways shown in Fig. 7b. These pathways are easier to visualise when considering CrO₆ coordination octahedra as in Fig. 7c. They are, from the shortest to the longest:

1. Along an edge shared between two CrO₆ octahedra in a basal plane (basal short, 2.7 Å);

2. Along an edge shared between two CrO₆ octahedra in a prismatic plane (prismatic short, 2.9 Å);

3. Along an edge in a prismatic plane that is not shared with another octahedron (prismatic long, 3.15 Å);

4. Along an edge in a basal plane that is not shared with another octahedron (basal long, 3.59 Å).

All references in Table 9 indicate that the first pathway, in the basal plane, is the most favourable, with an energy barrier of the same order of magnitude as those for chromium vacancy migration. The second most favourable pathway is in the prismatic plane. Together, these two mechanisms enable three-dimensional diffusion of oxygen vacancies. The dispersion of the results about diffusion pathways and energy barriers between authors does not allow quantitative validation of the potentials.

### 3.7. Grain boundaries

The corrosion layers are clearly polycrystalline [95] and grain boundaries are known to affect diffusion and crystal growth [95]. However, while there are numerous observations of grain boundaries, they are rarely characterised. Only a small number of simulations on grain boundaries in Cr₂O₃ have been published [48,50,96]. The literature is much more extensive for Al₂O₃, which has the same structure and for which numerous high-resolution TEM observations are available. Four types of low-index planes are associated with symmetric tilt grain boundaries in Al₂O₃: basal, pyramidal, rhombohedral, and prismatic [97–101]. Of these, rhombohedral [102] and prismatic [98] grain boundaries have been experimentally observed and characterised in Al₂O₃ using high-resolution TEM. Some of these grain boundaries in Cr₂O₃ were simulated in DFT [96] and with the empirical potential of Lawrence [50]. The grain boundaries reproduced with the empirical potentials are those of van der Geest et al., which follow the following nomenclature:

- The first lowercase letter corresponds to the type of plane, b: basal (0001), pr: prismatic (1010), py: pyramidal (1014), r: rhombohedral (1012);

- The second capital letter indicates the symmetry of the grain boundary, M: mirror, G: shear, S: screw, R: rotation;

- For some grain boundaries, the symbol in parentheses indicates the species occupying the grain boundary plane.

The energies of these grain boundaries obtained by DFT are presented in Table 10, along with those obtained using the empirical potential of Lawrence et al. There is no other source for this data, so any quantitative comparison is limited. However, the low-energy structures,

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Diffusion pathways for vacancies in Cr₂O₃: (a) Vcr; (b) Vo; (c) CrO₆ octahedron showing the Vo diffusion pathways Chromium ions appear as grey balls and oxygen ions as smaller red balls, and the diffusion pathways are coloured from the shortest (in yellow) to the longest (in brown).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 8. Migration energies of vacancies in Cr₂O₃ from DFT simulations [23]: a) V Cr; b) V O.</figcaption></figure>

which correspond to highly symmetric grain boundaries, should also be the most favourable with the empirical potentials.

The energies of the grain boundaries are not precisely known and are sourced only from one location. Therefore, they cannot be used for the quantitative validation of empirical potentials. However, several points

can be noted to qualitatively assess the plausibility of the potentials. Specifically, the prS grain boundary is found to be the most stable, followed by rS(V), prG, and bM.

<figure><img src="image_12.png" /><figcaption>Fig. 9. Thermal expansion of Cr₂O₃ from MD simulations using the different potentials: a) and b) lattice parameters; c) linear thermal expansion coefficient.</figcaption></figure>
# 4. Results

## 4.1. Structural and thermophysical properties

The lattice parameters at room temperature, calculated for the different potentials using molecular dynamics simulations, are presented in Table 11. Four potentials show good agreement with reference values, with a relative difference of less than 1.5 % for both *a* and *c*: those of Chartier, Owen, Pedone, and Sun. This difference is close to the dispersion in the experimental results, which is 1.19 % for *a* and 2.61 % for *c*. The Owen potential is particularly close to the experimental values for both *a* and *c*.

In general, the thermal expansion coefficients are markedly less accurate than the lattice parameters, with no relative error lower than 5 %. Experimental values are obtained by differentiating lattice parameters, which are measured more directly, resulting in increased uncertainty and dispersion of results. For empirical potentials, thermal expansion coefficients are rarely taken into account during their design, which is not surprising given the higher relative errors compared to lattice parameters. Therefore, it is more important to verify qualitative properties such as the order of magnitude of linear thermal expansion coefficients and the qualitative behaviour of $\alpha_a > \alpha_c$ rather than quantitative agreements with reference values. This is true for all potentials except those of Atkinson, Lawrence, and Sun. While no potential is very accurate for the thermal expansion parameters, those of Chartier and Pedone are the closest to experimental values.

The evolution of lattice parameters with temperature is shown in Fig. 9a and b, along with reference values for comparison [103–105]. The *a* parameter is well described by the Chartier and Owen potentials, while the *c* parameter is best reproduced by the Owen, Pedone, and Sun potentials. The Owen potential most accurately describes lattice parameters and their temperature dependence. Other potentials tend to behave poorly at high temperatures. For example, it was impossible to simulate the crystal above 1500 K with the Sun potential, as the initial structure melted during relaxation. Similarly, the structure becomes unstable beyond 770 K with the Lewis potential. The Barrera potential did not allow for dynamic simulation even at 300 K and was therefore disregarded for further calculations.

Literature data on thermal expansion often do not consider thermal expansion coefficients $\alpha_a$ and $\alpha_c$ separately, but instead provide a linear expansion coefficient. To obtain comparable values from molecular dynamics calculations, an effective linear thermal expansion coefficient of $\alpha_L = (2\alpha_a + \alpha_c)/3$ was calculated and is shown in Fig. 9c. Most potentials are relatively close to reference the data. The exceptions are the Lewis potential, which underestimates significantly thermal expansion, while the Owen potential tends to overestimate it.

The elastic constants and interatomic bond lengths were obtained from static calculations using the lattice parameters from the MD simulations, as described in Section 2. They are presented in Table 12 and Table 13, respectively. Their values can be compared to experimental measurements at 300 K for elastic constants (Alberts, 1976) and chosen reference values in Section 1 for bond lengths. The relative differences compared to the experimental values are more significant for the elastic constants than for lattice parameters, as expected when using empirical potentials. In particular, no potential predicts a negative value for constant $C_{14}$, which is the case experimentally. Three potentials have average errors around 30 %: those of Chartier, Owen, and Pedone. The Atkinson, Barrera, and Sun potentials predict poor elastic constants, with $C_{11}$ and $C_{14}$ in particular being significantly overestimated. The $C_{14}$ constant tends to be badly reproduced by all potentials and increases the average difference percentage because of its small magnitude. Ignoring it does not change the relative accuracy of the potentials: in this case, the Chartier, Owen and Pedone potentials have respective average errors of 7 %, 11 %, and 7 % respectively, whereas the other potentials all have average relative errors between than 22 % (Atkinson) and 43 % (Sun).

The interatomic distances calculated from static simulations using

lattice parameters fixed by molecular dynamics simulations at 300 K are presented in Table 13. The Owen potential best predicts these distances (with an average error of 1.20 %), followed by those of Chartier (1.27 %). All other potentials have errors greater than 2 %. Based on elastic properties and structural characteristics, these two potentials are significantly better in quality than the others.

The enthalpy increments of all potentials are very similar, particularly below 2000 K (Fig. 10a) and in good agreement with experimental values [71–73]. The potentials are very similar and in good agreement with experimental references for temperatures above approximately 800 K.

However, molecular dynamics simulations are less suitable for the calculation of thermodynamic properties at room temperature. Indeed, Fig. 10b suggests that all potentials significantly overestimate the specific heat below 700 K and remain at a plateau around $2.6 \times 10^{-4}$ eV/atom K at low temperatures. Lattice dynamics calculations in Fig. 10 with the same potentials yield much lower values, closer to literature values obtained by DFT [65], and tend towards 0 near absolute zero. Even when taking as a reference the extrapolation of Moore and Kelley's formula, which ignores the peak of the magnetic transition [73], MD values are overestimated by approximately 20 %. On the other hand, $C_V^{300}$ values obtained from lattice dynamics correspond better to those of Wang et al. obtained in the same manner but using DFT instead of empirical potentials. This makes any validation of potentials based on thermodynamic data difficult. Nevertheless, the thermodynamic functions from MD simulations are in good agreement with the literature in the 700 K-2000 K range. Most potentials are very similar with regard to specific heats and enthalpy increments obtained from MD calculations and this criterion cannot be used to distinguish them. However, the Atkinson potential gives a higher isochoric heat capacity in lattice dynamics simulations than the other potentials or the literature. The thermodynamic properties at 300 K, which can be validated by quantitative comparison with experimental results, are reported in Table 14. As expected, all the potential significantly over-estimate $C_p^{300}$. The Chartier and Pedone potentials predict values for $C_V^{300}$ that are within 2 % of the literature, but all potentials significantly under-estimate the vibrational entropy $S^{300}$.

The predicted phonon spectra used to calculate the specific heat $C_V^{300}$ are shown in Fig. 11. The curves are quite different, but those for the Chartier, Owen, and Pedone potentials have similar peaks: the first around 0.028 eV (or 7 THz) and the second around 0.045 eV. The first is directly visible in the DFT spectrum at 7 THz [65], while the second resembles the edge of the band gap, which appears at 14.2 THz (0.059 eV) in DFT and is shifted to lower energies with empirical potentials. No potential shows a bandgap like DFT calculations [65].

## 4.2. Polymorphs

The errors in predicted relative stabilities using the potentials compared to reference DFT values are in general around 0.2 eV/atom or less. The exception is the Pccn structure, with all potentials significantly underestimating its energy (Fig. 12). Some potentials predict structures more stable than $R\bar{3}c$, such as $Ia\bar{3}$, which is more stable according to the Atkinson, Barrera, Grimes, Lawrence and Sun potentials. The latter two also predict other structures more stable than $R\bar{3}c$: Fdd2, Pccn and Pna2$_1$ for Potential Sun, and $P\bar{3}m$1 for the Barrera potential.

None of the potentials were able to stabilise the C2/c structure, which transformed into $R\bar{3}c$ during energy minimization. This is consistent with DFT calculations [77] and experimental observation that this structure, which appears in compression at 15 GPa, disappears as pressure returns to ambient conditions.

The Chartier, Owen, and Pedone potentials provide relative stabilities in good agreement with DFT values. This is surprising given the variety of local environments and different coordination numbers across various structures. The Atkinson, Barrera, Lawrence, and Sun potentials
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Thermodynamic properties of Cr₂O₃ from MD simulations using the different potentials: a) enthalpy increment; b) isobaric heat capacity; c) isochoric heat capacity.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Vibration spectra of Cr₂O₃ calculated using the different empirical potentials.</figcaption></figure>

all predict polymorphs that are more stable than R$\bar{3}$c. Consequently, these potentials are not recommended for simulations where phenomena could result in a phase transition or local structural change, such as high-temperature simulations, those with numerous defects, or significant disturbances like irradiation. The instability of the R$\bar{3}$c structure is a plausible explanation for the temperature behaviour of the Sun potential noted in Section 4.1.

### 4.3. Point defects

#### 4.3.1. Chromium defects

The energy of Cr interstitials on tetrahedral sites could not be determined for any potential, as the interstitial moved to an adjacent octahedral site during energy minimization. The defect formation energies in Table 15 show that the triple defect (vacancy doublet with an interstitial, as depicted in Fig. 6b) is consistently more favourable than an isolated vacancy, in agreement with DFT calculations in the literature [23]. The Chartier, Owen, and Sun potentials predict that the most stable structure for Cr Frenkel pairs is with both defects aligned along
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Relative stability of the polymorphs of Cr₂O₃ calculated using the different empirical potentials (points) compared to reference values from DFT calculations (bars); the lines are a guide to the eye showing the Chartier and Owen potentials.</figcaption></figure>

the [0001] axis. This is consistent with the majority of studies in the literature. In contrast, the Atkinson, Lawrence, and Pedone potentials predict this structure to be unstable, with both defects recombining. They predict instead that the most stable structure is with the defects in the same basal plane. While this is consistent with one published study [27], the survey of the literature tends to indicate that the right structure is along the [0001] axis. For all potentials, the Frenkel pair formation energy is close to the value corresponding to infinitely separated defects beyond the first neighbour. This suggests that although there is significant attraction between interstitials and vacancies at very short distances, it quickly diminishes as the distance between defects increases.

### 4.3.2. Oxygen defects

Similarly to chromium interstitials, it was impossible to stabilise interstitials at tetrahedral sites, with observed interstitials occupying only octahedral sites. The interstitial dipole configuration is more favourable than the octahedral interstitial for all potentials except Sun. This finding aligns with most DFT calculations [28,88,92]. The formation energies of oxygen defects are in Table 16. The formation energy for Frenkel pairs where the interstitial and the vacancy are infinitely separated (noted “Far” in the table) fall in the range of the values from the literature for all potentials. Configurations where the interstitial is close to the vacancy (in the 5th nearest neighbour positions, noted “Close” in the table) were unstable for all potentials except Sun. In all other cases, the defects recombined during geometry optimisation. The formation energy for a close Frenkel pair using the Sun potential is still higher than the reference value by about 2.5 eV.

### 4.3.3. Schottky defects

The formation energies for Schottky defects are in Table 17. The most stable configuration is one with both chromium vacancies in a basal plane, except for Owen's potential, where the stable configuration has the chromium vacancies aligned along the [0001] axis. Given the discrepancies small number of DFT calculations in the literature and the spread of their results, it is impossible to determine which potential best represents these defects.

## 4.4. Diffusion pathways

Energy landscape comparisons for chromium vacancies show good qualitative agreement (Fig. 13a, c, and e). All three potentials accurately reproduce the intermediate state along pathway 4 and the shape of DFT curves, with one exception: Owen's potential exhibits a shoulder on

pathway 2, which remains unexplained. Migration energies are reported in Table 18. Empirical potentials suggest that the most favourable mechanism for chromium vacancy diffusion is pathway 4, contradicting DFT literature findings pointing to pathway 2. Energy barriers are of the same order of magnitude as references for Chartier and Pedone potentials but are underestimated by approximately 1 eV for Owen's potential. Error quantification is not provided in tables due to limited reference values and their relative disagreement.

The shape of the energy along the oxygen vacancies diffusion pathways are similar for the three potentials (Fig. 13b, d, and f). Contrary to chromium vacancies, oxygen vacancy energy barriers in Table 19 show good agreement in magnitude between potentials and DFT calculations from the literature, except for the Pedone potential, which finds path 2 more favourable than pathway 3. In line with chromium vacancies, the Chartier potential is closest to DFT energies, while Owen's potential tends to underestimate them.

## 4.5. Grain boundaries

Grain boundary formation energies, shown in Table 20, are of the correct order of magnitude, with average differences compared to DFT reference values [96] ranging between 40% (Chartier potential) and 76% (Pedone). The three most favourable structures in DFT also rank highly with empirical potentials. All three potentials seem reasonable given the lack of reference data available for more accurate validation.

# 5. Conclusion

The literature review in Section 1 was compiled into a dataset for empirical potential validation covering structural, thermodynamic, and elastic properties, as well as point defects and grain boundaries. Potentials originated from various literature sources with varying quality levels. Some were quickly ruled out due to poor performance during molecular dynamics simulations, such as those by Barrera, Lewis, and Sun. None of the three allowed for simulation over the entire temperature range of 300 K to the melting point of 2600 K: Barrera's potential could not be used for simulations at 300 K, while those of Lewis and Sun resulted in R3c structure instability beyond 770 K and 1400 K, respectively. Among the remaining potentials, those by Atkinson and Lawrence were discarded due to difficulties reproducing structural parameters or elastic constants. Some of these potentials were designed with polarisable shells, which were ignored in this study for practical reasons. Their performance with the shells enabled would be different than what
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Energy along the diffusion pathways of $V_{Cr}^{3-}$ and $V_{O}^{2+}$ vacancies in Cr$_2$O$_3$. The left column shows chromium vacancies and the right column shows oxygen vacancies. Each row corresponds to a different empirical potential, from top to bottom: Chartier, Owen, and Pedone.</figcaption></figure>
was presented in this article.

The remaining potentials, those of Chartier, Owen and Pedone, have proven satisfactory for most studied properties. Among these, Owen's potential has a slight advantage regarding structural properties but tends to underestimate migration energies of defects and grain boundary formation. The Chartier potential predicted slightly better defects formation and migration energies. The Pedone potential, whilst being close for all properties, did not show any particular strong point. Although it is not better than the two others for chromia, it is much more versatile and can be used to simulate many complex oxide crystals and glasses. Thus, it is advisable to consider the potentials of Chartier and Owen for molecular dynamics simulation of point and extended defect behaviour and self-diffusion.

### CRediT authorship contribution statement

Paul C.M. Fossati: Writing – review & editing, Writing – original

draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

### Declaration of competing interest

Comprehensive assessment of empirical potentials for molecular dynamics simulations of chromia.

The author declare that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

## Appendix A. Data tables

<caption>Table A.1</caption>

Lattice and internal parameters of Cr₂O₃ from the literature. The parameters are those of the 30-atom conventional hexagonal cell (6 formula units); values in brackets were calculated from the lattice parameters in the reference.

<table><thead><tr><th>Source</th><th>a (Å)</th><th>c (Å)</th><th>V (Å<sup>3</sup>)</th><th>x(O)</th><th>z(Cr)</th><th>Notes</th></tr></thead><tbody><tr><td>[106]</td><td>4.933</td><td>13.586</td><td>286.31</td><td>–</td><td>–</td><td>Exp.</td></tr><tr><td>[59]</td><td>4.941</td><td>13.578</td><td>287.08</td><td>–</td><td>–</td><td>Exp.</td></tr><tr><td>[3]</td><td>4.958</td><td>13.60</td><td>289.45</td><td>–</td><td>–</td><td>Powder</td></tr><tr><td>[107]</td><td>4.96</td><td>13.59</td><td>(289.54)</td><td>0.303</td><td>0.345</td><td>Exp.</td></tr><tr><td>[63]</td><td>4.9607</td><td>13.599</td><td>289.82</td><td>0.306</td><td>0.3475</td><td>Exp.</td></tr><tr><td>[58]</td><td>4.954</td><td>13.584</td><td>(288.72)</td><td>0.306</td><td>0.3475</td><td>Exp.</td></tr><tr><td>[108]</td><td>4.952</td><td>13.566</td><td>(288.10)</td><td>0.306</td><td>0.3476</td><td>Single crystal</td></tr><tr><td>[109]</td><td>4.961</td><td>13.600</td><td>(289.87)</td><td>0.306</td><td>0.3475</td><td>Exp.</td></tr><tr><td>[110]</td><td>–</td><td>–</td><td>288.6</td><td>–</td><td>–</td><td>EoS fit, powder</td></tr><tr><td>[64]</td><td>4.9507</td><td>13.5656</td><td>287.98</td><td>0.3051</td><td>0.3477</td><td>Single crystal</td></tr><tr><td>[111]</td><td>4.9570</td><td>13.5923</td><td>(286.10)</td><td>–</td><td>–</td><td>300 K</td></tr><tr><td>[61]</td><td>4.96004</td><td>13.5982</td><td>(289.72)</td><td>0.30590</td><td>0.34756</td><td>298 K</td></tr><tr><td>[61]</td><td>4.96195</td><td>13.5974</td><td>(289.93)</td><td>0.30583</td><td>0.34757</td><td>318 K</td></tr><tr><td>[82]</td><td>4.9566</td><td>13.595</td><td>289</td><td>–</td><td>–</td><td>Exp.</td></tr><tr><td>[112]</td><td>4.939</td><td>13.627</td><td>(287.88)</td><td>–</td><td>–</td><td>Thin film</td></tr><tr><td>[104]</td><td>4.95721</td><td>13.5917</td><td>(289.25)</td><td>0.30569</td><td>0.34734</td><td>300 K</td></tr><tr><td>[113]</td><td>4.93</td><td>13.507</td><td>288.6</td><td>0.3065</td><td>0.3475</td><td>EoS fit, 3.62 GPa</td></tr><tr><td>[114]</td><td>–</td><td>–</td><td>288.7</td><td>–</td><td>–</td><td></td></tr><tr><td>[115]</td><td>4.9526</td><td>13.5923</td><td>288.73</td><td>0.3058</td><td>0.3475</td><td>HP-HT synthesis</td></tr><tr><td>[115]</td><td>4.9398</td><td>13.5366</td><td>286.06</td><td>–</td><td>–</td><td>Powder</td></tr><tr><td>[116]</td><td>4.948</td><td>13.545</td><td>(287.19)</td><td>–</td><td>–</td><td>Thin film</td></tr><tr><td>[117]</td><td>4.919</td><td>13.589</td><td>(284.76)</td><td>–</td><td>–</td><td>Thin film</td></tr><tr><td>[118]</td><td>4.958</td><td>13.590</td><td>(287.31)</td><td>–</td><td>–</td><td>Exp.</td></tr><tr><td>[43]</td><td>5.048</td><td>13.735</td><td>(303.10)</td><td>0.30060</td><td>0.35055</td><td>Hartree-Fock</td></tr><tr><td>[82]</td><td>–</td><td>–</td><td>286.2</td><td>0.552</td><td>0.350</td><td>FPLMTO</td></tr><tr><td>[76]</td><td>4.970</td><td>13.600</td><td>(290.93)</td><td>0.557</td><td>0.347</td><td>LSDA, EoS fit</td></tr><tr><td>[119]</td><td>4.941</td><td>13.829</td><td>(292.38)</td><td>–</td><td>–</td><td>DFT (GGA)</td></tr><tr><td>[119]</td><td>5.073</td><td>13.839</td><td>(308.44)</td><td>–</td><td>–</td><td>DFT (GGA + U)</td></tr><tr><td>[120]</td><td>4.9598</td><td>13.5894</td><td>(289.51)</td><td>–</td><td>–</td><td>DFT (LDA)</td></tr><tr><td>[65]</td><td>–</td><td>–</td><td>289.323</td><td>–</td><td>–</td><td>DFT + U, QHA (300K)</td></tr><tr><td>[66]</td><td>5.00</td><td>13.51</td><td>293.02</td><td>–</td><td>–</td><td>DFT (GGA + U)</td></tr><tr><td>[77]</td><td>5.038</td><td>13.625</td><td>299.70</td><td>0.3014</td><td>0.3503</td><td>DFT (GGA)</td></tr><tr><td>[27]</td><td>5.074</td><td>13.850</td><td>(308.80)</td><td>–</td><td>–</td><td>DFT (GGA)</td></tr><tr><td>[121]</td><td>5.08</td><td>13.99</td><td>(312.66)</td><td>–</td><td>–</td><td>DFT (GGA + U)</td></tr><tr><td>[23]</td><td>–</td><td>13.680</td><td>–</td><td>–</td><td>–</td><td>DFT (GGA + U)</td></tr><tr><td>[122]</td><td>5.07</td><td>13.85</td><td>(308.31)</td><td>–</td><td>–</td><td>DFT (GGA + U)</td></tr></tbody></table>
<caption>**Table A.2**</caption>

Interatomic distances and bulk modulus of Cr₂O₃ from the literature; values in brackets were calculated from the structural parameters in the reference.

<table><caption>**Table A.3**</caption><thead><tr><th>Source</th><th>Cr-Cr (Å)</th><th>Cr-O (Å)</th><th>O-O (Å)</th><th>B (GPa)</th><th>Notes</th></tr></thead><tbody><tr><td>[59]</td><td>–</td><td>1.96</td><td>2.01</td><td>2.57</td><td>–</td><td>Exp.</td></tr><tr><td>[107]</td><td>(2.58)</td><td>(2.88)</td><td>(1.98)</td><td>(1.99)</td><td>(2.60)</td><td>Exp.</td></tr><tr><td>[63]</td><td>2.70</td><td>2.88</td><td>1.97</td><td>2.02</td><td>2.63</td><td>Exp.</td></tr><tr><td>[108]</td><td>2.66</td><td>(2.89)</td><td>1.97</td><td>2.02</td><td>2.64</td><td>Single crystal</td></tr><tr><td>[58]</td><td>(2.65)</td><td>(2.89)</td><td>1.97</td><td>2.02</td><td>(2.63)</td><td>Exp.</td></tr><tr><td>[109]</td><td>(2.65)</td><td>(2.89)</td><td>1.965</td><td>2.016</td><td>(2.63)</td><td>Exp.</td></tr><tr><td>[123]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>Exp.</td></tr><tr><td>[110]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>Exp.</td></tr><tr><td>[64]</td><td>2.6486</td><td>2.8848</td><td>1.962</td><td>2.009</td><td>2.617</td><td>231</td><td>EoS fit, powder</td></tr><tr><td>[124]</td><td>2.6451</td><td>2.8854</td><td>1.9617</td><td>2.0116</td><td>2.6252</td><td>238</td><td>Single crystal</td></tr><tr><td>[61]</td><td>2.6533</td><td>(2.89)</td><td>(1.96)</td><td>2.0155</td><td>(2.63)</td><td>2.7319</td><td>Based on [58]</td></tr><tr><td>[61]</td><td>2.6534</td><td>(2.89)</td><td>(1.97)</td><td>2.0157</td><td>(2.63)</td><td>(2.74)</td><td>298 K</td></tr><tr><td>[82]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>(2.74)</td><td>318 K</td></tr><tr><td>[104]</td><td>(2.65)</td><td>(2.89)</td><td>(1.97)</td><td>(2.01)</td><td>(2.62)</td><td>240</td><td>Exp.</td></tr><tr><td>[113]</td><td>(2.63)</td><td>(2.87)</td><td>1.957</td><td>2.007</td><td>(2.62)</td><td>(2.74)</td><td>300 K</td></tr><tr><td>[114]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>(2.62)</td><td>(2.72)</td><td>EoS fit, 3.62 GPa</td></tr><tr><td>[115]</td><td>(2.65)</td><td>(2.89)</td><td>(1.96)</td><td>(2.01)</td><td>(2.62)</td><td>(2.74)</td><td>HP-HT synthesis</td></tr><tr><td>[83]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>245</td></tr><tr><td>[43]</td><td>(2.76)</td><td>(2.95)</td><td>1.990</td><td>2.052</td><td>(2.62)</td><td>(2.76)</td><td>262</td><td>Hartree-Fock</td></tr><tr><td>[78]</td><td>2.745</td><td>2.946</td><td>1.980</td><td>2.048</td><td>–</td><td>–</td><td>LDA</td></tr><tr><td>[82]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>FPLMTO</td></tr><tr><td>[76]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>LSDA, EoS fit</td></tr><tr><td>[119]</td><td>2.642</td><td>4.271</td><td>–</td><td>–</td><td>–</td><td>–</td><td>DFT (GGA)</td></tr><tr><td>[119]</td><td>2.723</td><td>4.271</td><td>–</td><td>–</td><td>–</td><td>–</td><td>DFT (GGA + U)</td></tr><tr><td>[65]</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>238</td><td>DFT + U, QHA (300 K)</td></tr><tr><td>[66]</td><td>–</td><td>–</td><td>1.97</td><td>2.02</td><td>–</td><td>–</td><td>203</td><td>DFT (GGA + U)</td></tr><tr><td>[77]</td><td>(2.73)</td><td>(2.95)</td><td>(1.98)</td><td>(2.04)</td><td>2.63</td><td>(2.75)</td><td>DFT (GGA)</td></tr><tr><td>[27]</td><td>2.717</td><td>4.027</td><td>–</td><td>–</td><td>–</td><td>–</td><td>DFT (GGA)</td></tr><tr><td>[121]</td><td>–</td><td>–</td><td>2.02</td><td>2.05</td><td>–</td><td>–</td><td>DFT (GGA + U)</td></tr></tbody></table>

Elastic constants of Cr₂O₃ from the literature.

<table><caption>**Table A.4**</caption><thead><tr><th>Source</th><th>C<sub>11</sub> [GPa]</th><th>C<sub>12</sub> [GPa]</th><th>C<sub>13</sub> [GPa]</th><th>C<sub>14</sub> [GPa]</th><th>C<sub>33</sub> [GPa]</th><th>C<sub>44</sub> [GPa]</th><th>Notes</th></tr></thead><tbody><tr><td>[60]</td><td>374</td><td>148</td><td>175</td><td>-19</td><td>363</td><td>159</td><td>Exp., 20 °C</td></tr><tr><td>[65]</td><td>373</td><td>160</td><td>178</td><td>-21</td><td>349</td><td>160</td><td>DFT QHA, 300 K</td></tr><tr><td>[66]</td><td>348</td><td>142</td><td>137</td><td>-8</td><td>308</td><td>136</td><td>DFT + U, 0 K</td></tr></tbody></table>

Room-temperature thermodynamical properties of Cr₂O₃ from the literature.

<table><thead><tr><th>Source</th><th>C<sub>p</sub><sup>300</sup> (J/mol Cr<sub>2</sub>O<sub>3</sub> · K)</th><th>S<sub>p</sub><sup>300</sup> (J/mol Cr<sub>2</sub>O<sub>3</sub> · K)</th><th>Notes</th></tr></thead><tbody><tr><td>[73]</td><td>104.8</td><td>–</td><td>Calc.</td></tr><tr><td>[5]</td><td>118.74</td><td>–</td><td>Exp.</td></tr><tr><td>[71]</td><td>121.922</td><td>80.65</td><td>Exp.</td></tr><tr><td>[74]</td><td>127.08</td><td>82.3</td><td>Calc.</td></tr><tr><td>[70]</td><td>112.0</td><td>82.8</td><td>Calc.</td></tr><tr><td>[67]</td><td>121.0</td><td>81.1</td><td>Exp.</td></tr><tr><td>[125]</td><td>–</td><td>85.74</td><td>Exp.</td></tr><tr><td>[72]</td><td>121.5</td><td>80.95</td><td>Exp.</td></tr><tr><td>[65]</td><td>86.65</td><td>–</td><td>DFT QHA</td></tr></tbody></table>
<table><caption>**Table A.5**
Room-temperature thermal expansion of Cr₂O₃ from the literature; values in brackets were calculated from other properties in the source.</caption><thead><tr><th rowspan="2">Source</th><th colspan="2">α<sub>L</sub> (10<sup>-6</sup> K<sup>-1</sup>)</th><th rowspan="2">α<sub>V</sub> (10<sup>-5</sup> K<sup>-1</sup>)</th><th rowspan="2">Notes</th></tr><tr><th>a</th><th>c</th></tr></thead><tbody><tr><td>[126]</td><td></td><td>(8.42)</td><td>-</td><td>Calculated from elongation between 300 K and 800 K</td></tr><tr><td>[127]</td><td></td><td>(8.9)</td><td>(6.31)</td><td>Calculated from elongation between 300 K and 1200 K</td></tr><tr><td>[105]</td><td></td><td>(8.39)</td><td>(7.23)</td><td>Calculated from lattice parameters at 298 K and 1743 K</td></tr><tr><td>[5]</td><td></td><td>9.6</td><td>-</td><td></td></tr><tr><td>[60]</td><td>7.1</td><td></td><td>1.93</td><td>Room temperature, PM</td></tr><tr><td>[128]</td><td></td><td>5.7</td><td>-</td><td></td></tr><tr><td>[103]</td><td></td><td>(9.64)</td><td>(9.33)</td><td>Calculated from lattice parameters at 298 K and 842 K</td></tr><tr><td>[129]</td><td></td><td>-</td><td>2.95</td><td></td></tr><tr><td>[75]</td><td></td><td>6.03</td><td>1.59</td><td></td></tr><tr><td>[114]</td><td></td><td>-</td><td>1.86</td><td></td></tr><tr><td>[105]</td><td></td><td>(7.89)</td><td>(6.40)</td><td>Calculated from lattice parameters at 298 K and 1343 K</td></tr><tr><td>[65]</td><td></td><td>4.12</td><td>-</td><td>DFT QHA</td></tr><tr><td>[118]</td><td></td><td>5.32</td><td>-</td><td></td></tr></tbody></table>

References

[1] U. Ehrnstén, Corrosion and stress corrosion cracking of austenitic stainless steels, in: Comprehensive Nuclear Materials, Elsevier, 2020, pp. 118–128, 10/mc4g.

[2] F. Cattant, D. Crusset, D. Féron, Corrosion issues in nuclear industry today, Mater. Today 11 (10) (2008) 32–37, 10/cm9xw3.

[3] O. Kouvo, Y. Vuorelainen, Eskolaite, a new chromium mineral, Am. Mineral. 43 (1958) 1098–1106.

[4] N.Y. Toker, L.S. Darken, A. Muan, Equilibrium phase relations and thermodynamics of the Cr-O system in the temperature range of 1500°C to 1825°C, Metal. Trans. B 22B (1991) 225–232, 10/b9rfzr.

[5] I.S. Kainarskii, E.V. Degtyareva, Chromic oxide as refractory material, Refractories 18 (1–2) (1977) 42–47, 10/c4phf7.

[6] D. Rodbell, R. DeVries, On the decomposition of CrO₂ in air, Mater. Res. Bull. 2 (4) (1967) 491–495, 10/fmz78b.

[7] W.K. Jóžwiak, W. Ignaczak, D. Dominiak, T.P. Maniecki, Thermal stability of bulk and silica supported chromium trioxide, Appl. Catal. Gen. 258 (1) (2004) 33–45, 10/bf3jtp.

[8] C. Greskovich, Deviation from stoichiometry in Cr₂O₃ at high oxygen partial pressures, J. Am. Ceram. Soc. 67 (6) (1984) C-111-C-112, 10/dpb4kbp.

[9] C. Guéneau, et al., TAF-ID: an international thermodynamic database for nuclear fuels applications, Calphad 72 (2021) 102212, 10/gjfq73.

[10] L. Kjellqvist, M. Selleby, B. Sundman, Thermodynamic modelling of the Cr-Fe-Ni-O system, Calphad 32 (3) (2008) 577–592, 10/df4hf3.

[11] E.N. Bunting, Phase equilibria in the system Cr₂O₃-Al₂O₃, Bur. Stand. J. Res. 6 (6) (1931) 947, 10/gs4z7f.

[12] S.J. Burnett, 'Properties of Refractory Materials', United Kingdom Atomic Energy Authority, Harwell, Technical Report AERE-R-4657, 1964.

[13] P.T.B. Schaffer, Dichromium trioxide, in: Handbook of High Temperature Materials 1, Plenum Press, New York, 1964, p. 331.

[14] R.N. McNally, F.I. Peters, P.H. Ribbe, Laboratory furnace for studies in controlled atmospheres; melting points of MgO in a N₂ atmosphere and of Cr₂O₃ in N₂ and in air atmospheres, J. Am. Ceram. Soc. 44 (10) (1961) 491–493, 10/c5zs97.

[15] M. Kowalski, P.J. Spencer, Thermodynamic reevaluation of the Cr-O, Fe-O and Ni-O systems: remodelling of the liquid, BCC and FCC phases, Calphad 19 (3) (1995) 229–243, 10/df7v4b.

[16] J.R. Taylor, A.T. Dinsdale, A thermodynamic assessment of the Ni-O, Cr-O and Cr-Ni-O systems using the ionic liquid and compound energy models, Z. Metallkd. 81 (5) (1990) 354–366, 10/gtd4gd.

[17] P. Schmuki, From bacon to barriers: a review on the passivity of metals and alloys, J. Solid State Electrochem. 6 (3) (2002) 145–164, 10/bfd86w.

[18] A.C.S. Sabioni, B. Lesage, A.M. Huntz, J.C. Pivin, C. Monty, Self-diffusion in Cr₂O₃ I. Chromium diffusion in single crystals, Phil. Mag. A 66 (3) (1992) 333–350, 10/ cjbwcz.

[19] X. Huang, L. Martinelli, S. Bosonnet, P.C.M. Fossati, L. Latu-Romain, Y. Wouters, Chromium depletion in a Ni-30Cr alloy during high-temperature oxidation, High Temp. Corrosion Mater. 100 (5–6) (2023) 745–773, 10/gtb9zv.

[20] E. Schmucker, C. Petitjean, L. Martinelli, P.-J. Panteix, S. Ben Lagha, M. Vilasi, Oxidation of Ni-Cr alloy at intermediate oxygen pressures. I. Diffusion mechanisms through the oxide layer, Corros. Sci. 111 (2016) 474–485, 10/ f83wws.

[21] P. Kofstad, K.P. Lillerud, On high temperature oxidation of chromium: II. Properties of Cr₂O₃ and the oxidation mechanism of chromium, J. Electrochem. Soc. 127 (11) (1980) 2410–2419, 10/cnsbx3.

[22] A. Holt, P. Kofstad, Electrical conductivity and defect structure of Cr₂O₃. I. High temperatures (>~1000°C), Solid State Ion. 69 (1994) 127–136, 10/bd7sm4.

[23] B.K. Medasani, M.L. Sushko, K.M. Rosso, D.K. Schreiber, S.M. Bruemmer, Vacancies and vacancy-mediated self diffusion in Cr₂O₃: a first-principles study, J. Phys. Chem. C 121 (3) (2017) 1817–1831, 10/f9rzdw.

[24] L. Latu-Romain, et al., The role of oxygen partial pressure on the nature of the oxide scale on a NiCr model alloy, Oxid. Metal. 88 (3–4) (2017) 481–493, 10/ gb93b.

[25] P. Cao, D. Wells, M.P. Short, Anisotropic ion diffusion in α-Cr₂O₃: an atomistic simulation study, Phys. Chem. Chem. Phys. 19 (21) (2017) 13658–13663, 10/ gtd4d9.

[26] C. Gray, Y. Lei, G. Wang, Charged vacancy diffusion in chromium oxide crystal: DFT and DFT+U predictions, J. Appl. Phys. 120 (21) (2016) 215101, 10/gs4v8n.

[27] F. Lebreau, M.M. Islam, B. Diawara, P. Marcus, Structural, magnetic, electronic, defect, and diffusion properties of Cr₂O₃: a DFT+U study, J. Phys. Chem. C 118 (31) (2014) 18133–18145, 10/f6ds4f.

[28] B.K. Medasani, M.L. Sushko, K.M. Rosso, D.K. Schreiber, S.M. Bruemmer, First-principles investigation of native interstitial diffusion in Cr₂O₃, J. Phys. Chem. C 122 (24) (2018) 12984–12993, 10/gdtf6f.

[29] J. Vaari, Molecular dynamics simulations of vacancy diffusion in chromium(III) oxide, hematite, magnetite and chromite, Solid State Ion. 270 (2015) 10–17, 10/ ghfxbh.

[30] J. Wang, D. Shin, S. Shin, Comprehensive evaluation and parametric sensitivity of interatomic potential models for diffusion kinetics of Cr₂O₃ in molecular dynamics, AIP Adv. 9 (1) (2019) 015123, 10/gs4v8c.

[31] A.P. Thompson, et al., LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun. 271 (2022) 108171, 10/gmwrb.

[32] W.G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A 31 (3) (1985) 1695–1697, 10/dw65r7.

[33] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81 (1) (1984) 511–519, 10/bcp3gb.

[34] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (12) (1981) 7182–7190, 10/ bjn99x.

[35] M.E. Tuckerman, J. Alejandro, R. López-Rendón, A.L. Jochim, G.J. Martyna, A Liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal-isobaric ensemble, J. Phys. A 39 (19) (2006) 5629–5651, 10/c9fr6d.

[36] J.D. Gale, A.L. Rohl, The general utility lattice program (GULP), Mol. Sim. 29 (5) (2003) 291–341, 10/cv9tnq.

[37] G. Henkelman, B.P. Überuaga, J.ónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113 (22) (2000) 9901–9904, 10/fmm8n.

[38] D. Wolf, P. Keblinski, S.R. Phillipot, J. Eggebrecht, Exact method for the simulation of coulombic systems by spherically truncated, pairwise $r^{-1}$ summation, J. Chem. Phys. 110 (17) (1999) 8254–8282, 10/cfcxdk.

[39] K.J.W. Atkinson, R.W. Grimes, M.R. Levy, Z.L. Coul, T. English, Accommodation of impurities in α-Al₂O₃, α-Cr₂O₃ and α-Fe₂O₃, J. Eur. Ceram. Soc. 23 (16) (2003) 3059–3070, 10/d98qs9.

[40] R.W. Grimes, D.J. Binks, A.B. Lidiard, The extent of zinc oxide solution in zinc chromate spinel, Phil. Mag. A 72 (3) (1995) 651–668, 10/bsgcnx.

[41] R.W. Grimes, Solution of MgO, CaO, and TiO₂ in α-Al₂O₃, J. Am. Ceram. Soc. 77 (2) (1994) 378–384, 10/dqwg3.

[42] M.A. San Miguel Barrera, J.F. Sanz, L.J. Álvarez, J.A. Odriozola, Molecular-dynamics simulations of premelting processes in Cr₂O₃, Phys. Rev. B 58 (10) (1998) 6057–6062, 10/cw43qw.

[43] M. Catti, G. Sandrone, G. Valerio, R. Dovesi, Electronic, magnetic and crystal structure of Cr₂O₃ by theoretical methods, J. Phys. Chem. Solid 57 (11) (1996) 1735–1741, 10/fhppk7.

[44] A. Chartier, B. Golovchuk, S. Gossé, L. Van Brutzel, Disordering and grain boundaries of (Ni,Fe)Cr₂O₄ spinels from atomistic calculations, J. Chem. Phys. 139 (13) (2013) 134702, 10/gs4v8k.

[45] S. Morooka, S. Zhang, T. Nishikawa, H. Awaji, Potential model parameters for molecular dynamics simulation of alumina-magnesia systems, J. Ceram. Soc. Japan 107 (1252) (1999) 1225–1228, 10/fwdqms.
[46] L. Minervini, Defect cluster formation in $M_2O_3$-doped CeO$_2$, Solid State Ion. 116 (3-4) (1999) 339-349, 10/bq593j.

[47] J.D. Gale, C.R.A. Catlow, W.C. Mackrodt, Periodic ab initio determination of interatomic potentials for alumina, Model. Simul. Mater. Sci. Eng. 1 (1) (1992) 73-81, 10/bn429t.

[48] P.J. Lawrence, S.C. Parker, P.W. Tasker, Computer simulation studies of perfect and defective surfaces in $Cr_2O_3$, J. Am. Ceram. Soc. 71 (8) (1988) C-389-C-391, 10/bwqz7p.

[49] G.V. Lewis, C.R.A. Catlow, Potential models for ionic oxides, J. Phys. C Solid State Phys. 18 (6) (1985) 1149-1161, 10/drspt8t.

[50] P.J. Lawrence, S.C. Parker, Computer modelling of oxide surfaces and interfaces, in: C.R.A. Catlow, S.C. Parker, M.P. Allen (Eds.), Computer Modelling of Fluids Polymers and Solids, Springer Netherlands, Dordrecht, 1990, pp. 219-248, 10/ mc4q.

[51] P.J. Lawrence, S.C. Parker, P.W. Tasker, Computer modeling of the defect properties of chromium oxide, $Cr_2O_3$, in Nonstoichiometric compounds, in: C.R. A. Catlow, W.C. Mackrodt (Eds.), Advances in Ceramics 23, The American Ceramics Society, 1987, pp. 247-256.

[52] C.R.A. Catlow, Point defect and electronic properties of uranium dioxide, Proc. R. Soc. Lond. A 353 (1675) (1977) 533-561, 10/dn4rmr.

[53] M.W. Owen, M.W.D. Cooper, M.J.D. Rushton, A. Claisse, W.E. Lee, S. C. Middleburgh, Diffusion in undoped and Cr-doped amorphous $UO_2$, J. Nucl. Mater. 576 (2023) 154270, 10/gtd4df2.

[54] M.W.D. Cooper, M.J.D. Rushton, R.W. Grimes, A many-body potential approach to modelling the thermomechanical properties of actinide oxides, J. Phys. Condens. Matter 26 (10) (2014) 105401, 10/grsbzr.

[55] A. Pedone, G. Malavasi, M.C. Menziani, A.N. Cormack, U. Segre, A new self-consistent empirical interatomic potential model for oxides, silicates, and silica-based glasses, J. Phys. Chem. B 110 (24) (2006) 11780-11795, 10/ddgp2v.

[56] J. Sun, T. Stirner, A. Matthews, Calculation of native defect energies in α-$Al_2O_3$ and α-$Cr_2O_3$ using a modified Matsui potential, Surf. Coat. Tech. 201 (7) (2006) 4201-4204, 10/bx3v46.

[57] M. Matsui, Molecular dynamics study of the structures and bulk moduli of crystals in the system CaO-MgO-$Al_2O_3$-SiO$_2$, Phys. Chem. Minerals 23 (1996) 345-353, 10/d4g94p.

[58] R.W.G. Wyckoff, Compounds $R_2X_3$, in: Crystal Structures, 2nd ed. 1, Interscience, New York, 1964, pp. 1-44.

[59] W.H. Zachariasen, Untersuchung über die Kristallstruktur von Sesquioxyden und Verbindungen $ABO_3$, in: Skrifter utgitt av Det Norske Videnskaps-Akademi i Oslo 4, Jacob Dybwad, Oslo, 1928.

[60] H.L. Alberts, J.C.A. Boeyen, The elastic constants and distance dependence of the magnetic interactions of $Cr_2O_3$, J. Magnetism Magnetic Mater. 2 (1976) 327-333, 10/fg5tw2.

[61] M. Baster, F. Bourée, A. Kowalska, Z. Latacz, The change of crystal and exchange parameters in the vicinity of $T_N$ in $Cr_2O_3$, J. Alloy. Comp. 296 (1-2) (2000) 1-5, 10/b74crc.

[62] S. Greenwald, Changes in lattice constant of $Cr_2O_3$ near the Curie temperature, Nature 177 (4502) (1956) 286-287, 10/cm2j7w.

[63] E.E. Newnham, Y.M.D. Haan, Refinement of the α-$Al_2O_3$, Ti$_2O_3$, V$_2O_3$ and $Cr_2O_3$ structures, Z. Krist. 117 (1-6) (1962) 235-237, 10/gtd4f2.

[64] L.W. Finger, R.M. Hazen, Crystal structure and isothermal compression of $Fe_2O_3$, $Cr_2O_3$, and $V_2O_3$ to 50 kbars, J. Appl. Phys. 51 (10) (1980) 5362, 10/fd5q4r.

[65] Y. Wang, et al., First-principles lattice dynamics, thermodynamics, and elasticity of $Cr_2O_3$, Surf. Sci. 606 (17-18) (2012) 1422-1425, 10/f35q79.

[66] A. Jain, et al., Commentary: the materials project: a materials genome approach to accelerating materials innovation, APL Mater. 1 (1) (2013) 011002, 10/ gdwpk8.

[67] C.T. Anderson, The heat capacities of chromium, chromic oxide, chromous chloride and chromic chloride at low temperatures, J. Am. Chem. Soc. 59 (3) (1937) 488-491, 10/cx5fvc.

[68] R.H. Bruce, D.S. Gannell, Specific heat of $Cr_2O_3$ near the Néel temperature, Phys. Rev. B 15 (9) (1977) 4451-4459, 10/d6n9sj.

[69] J. Volger, Anomalous specific heat of chromium oxide ($Cr_2O_3$) at the antiferromagnetic Curie temperature, Nature 170 (4337) (1952) 1027, 10/ dn726v.

[70] S. Klemme, H. St, C. O'Neill, W. Schnelle, E. Gmelin, The heat capacity of $MgCr_2O_4$, FeCr$_2O_4$, and $Cr_2O_3$ at low temperatures and derived thermodynamic properties, Amer. Mineral. 85 (11-12) (2000) 1686-1693, 10/gtd4fr.

[71] M.W. Chase, NIST-JANAF themofemical tables, fourth edition, J. Phys. Chem. Ref. Data Monogr. 9 (1998).

[72] V.M. Gurevich, O.L. Kuskov, N.N. Smirnov, K.S. Gavrichev, A.V. Markin, Thermodynamic functions of eskolaite $Cr_2O_3(c)$ at 0-1800 K, Geochem. Int. 47 (12) (2009) 1170-1179, 10/dj4pxw.

[73] G.E. Moore, K.K. Kelly, High-Temperature Heat Content of the Chromium Carbides and of Chromic Oxide, US Bureau of Mines, 1944, Technical paper 662.

[74] S.E. Ziemniak, L.M. Anovitz, R.A. Castelli, W.D. Porter, Thermodynamics of $Cr_2O_3$, FeCr$_2O_4$, ZnCr$_2O_4$, and $CoCr_2O_4$, J. Chem. Thermodyn. 39 (11) (2007) 1474-1492, 10/dq3g5c.

[75] A.-M. Huntz, A. Piant, J.L. Lebrun, S. Daghigh, Evidence of stress relaxation in thermally grown oxide layers—experiments and modelling, Mater. Sci. Eng. A 248 (1-2) (1998) 44-55, 10/bh2qc.

[76] A.Yu. Dobin, W. Duan, R.M. Wentzcovitch, Magnet Structural effects and phase transition in $Cr_2O_3$ under pressure, Phys. Rev. B 62 (18) (2000) 11997-12000, 10/fk92vn.

[77] C. Wessel, R. Dronskowski, A first-principles study on chromium sesquioxide, $Cr_2O_3$, J. Solid State Chem. 199 (2013) 149-153, 10/f4q7dw.

[78] W. Duan, G. Paiva, R.M. Wentzcovitch, A. Fazzio, Optical transitions in ruby across the corundum to $Rh_2O_3$ (II) phase transformation, Phys. Rev. Lett. 81 (15) (1998) 3267-3270, 10/cdv6zs.

[79] J. Mougin, T. Le Bihan, G. Lucazeau, High-pressure study of $Cr_2O_3$ obtained by high-temperature oxidation by X-ray diffraction and Raman spectroscopy, J. Phys. Chem. Solid 62 (3) (2001) 553-563, 10/cx82nz.

[80] S.-H. Shim, T.S. Duffy, R. Jeanloz, C.-S. Yoo, V. Iota, Raman spectroscopy and X-ray diffraction of phase transitions in $Cr_2O_3$ to 61 GPa, Phys. Rev. B 69 (14) (2004) 144107, 10/c4cqd.

[81] P.D. Dernier, M. Marezio, Crystal structure of the low-temperature antiferromagnetic phase of $Cr_2O_3$, Phys. Rev. B 2 (9) (1970) 3771-3776, 10/ bp7zzc.

[82] S. Rekhi, L.S. Dubrovinsky, R. Ahuja, S.K. Saxena, B. Johansson, Experimental and theoretical investigations on eskolaite ($Cr_2O_3$) at high pressures, J. Alloy. Comp. 302 (1-2) (2000) 16-20, 10/b3q79p.

[83] A. Kantor, et al., High-pressure structural studies of eskolaite by means of single-crystal X-ray diffraction, Amer. Mineral. 97 (10) (2012) 1764-1770, 10/f4cn89.

[84] N.O. Golosova, et al., Structural and magnetic properties of $Cr_2O_3$ at high pressure, J. Alloy. Comp. 722 (2017) 593-598, 10/gf2bn8.

[85] J.E. Saal, S. Kirklin, M. Aykol, B. Meredig, C.卢overton, Materials design and discovery with high-throughput density functional theory: the open quantum materials database (OQMD), JOM 65 (11) (2013) 1501-1509, 10/f5h4rh.

[86] A. Holt, P. Kofstad, Electrical conductivity and defect structure of $Cr_2O_3$. II. Reduced temperatures (<~1000°C), Solid State Ion. 69 (1994) 137-143, 10/ bd7sm4.

[87] E.W.A. Young, J.H. Gerretsen, J.H.W. De Wit, The oxygen partial pressure dependence of the defect structure of chromium(III) oxide, J. Electrochem. Soc. 134 (9) (1987) 2257-2260, 10/dbm575.

[88] J.J. Carey, M. Nolan, Influence of trivalent doping on point and Frenkel defect formation in bulk chromium (III) oxide, Solid State Ion. 307 (2017) 51-64, 10/ gs4v8j.

[89] N.F. Mott, M.J. Littleton, Conduction in polar crystals. I. Electrolytic conduction in solid salts, Trans. Faraday Soc. 34 (1938) 485-499, 10/b72t3b.

[90] B.K. Medasani, M.L. Sushko, K.M. Rosso, D.K. Schreiber, S.M. Bruemmer, Temperature dependence of self-diffusion in $Cr_2O_3$ from first principles, J. Phys. Chem. C 123 (36) (2019) 22139-22150, 10/gs4v8v.

[91] C.R.A. Catlow, J. Corish, J. Hennessy, W.C. Mackrodt, Atomistic simulation of defect structures and ion transport in α-$Fe_2O_3$ and α-$Cr_2O_3$, J. Am. Ceram. Soc. 71 (1) (1988) 42-49, 10/fqpkw2.

[92] A.B. Kehoe, E. Arca, D.O. Scanlon, I.V. Shvets, G.W. Watson, Assessing the potential of mg-doped $Cr_2O_3$ as a novel *p*-type transparent conducting oxide, J. Phys. Condens. Matter 28 (12) (2016) 125501, 10/gs4v8s.

[93] H. Yu, C. Chen, R. Jiang, P. Qiu, Y. Li, Migration of ion vacancy in hydroxylated oxide film formed on Cr: a density functional theory investigation, J. Phys. Chem. C 116 (48) (2012) 25478-25485, 10/f4g82c.

[94] X. Huang, P.C.M. Fossati, L. Martinelli, S. Bosonnet, L. Latu-Romain, Y. Wouters, A DFT study of defects in paramagnetic $Cr_2O_3$, Phys. Chem. Chem. Phys. 24 (17) (2022) 10488-10498, 10/gs4v8q.

[95] X. Huang, Défauts ponctuels dans $Cr_2O_3$ et oxydation d'un alliage Ni-30Cr à haute température: approches expérimentale et numérique, Université Grenoble Alpes, Grenoble, 2022 [Online]. Available: https://theses.hal.science/tel-03775143/.

[96] A.G. Van Der Geest, M.M. Islam, T. Couvant, B. Diawara, Energy ordering of grain boundaries in $Cr_2O_3$: insights from theory, J. Phys. Condens. Matter 25 (48) (2013) 485005, 10/gs4v8b.

[97] S. Fabris, C. Elsässer, Σ13 (10-14) twin in α-$Al_2O_3$: a model for a general grain boundary, Phys. Rev. B 64 (24) (2001) 245117, 10/ckwkh4.

[98] S. Fabris, S. Nufer, C. Elsässer, T. Gemming, Prismatic Σ3 (10-10) twin boundary in α-$Al_2O_3$ investigated by density functional theory and transmission electron microscopy, Phys. Rev. B 66 (15) (2002) 155415, 10/c4z62d.

[99] T. Geipel, K.P.D. Lagerl, P. Pirouz, A.H. Heuer, A zonal dislocation mechanism for rhomboedratal twinning in sapphire (α-$Al_2O_3$), Acta Metall. Mater. 42 (4) (1994) 1367-1372, 10/bp78z4.

[100] A.G. Marinopoulos, S. Nufer, C. Elsässer, Interfacial structures and energetics of basal twins in α-$Al_2O_3$: first-principles density-functional and empirical calculations, Phys. Rev. B 63 (16) (2001) 165112, 10/fh5sgx.

[101] A.G. Marinopoulos, C. Elsässer, Microscopic structure and bonding at the rhombohedral twin interface in α-$Al_2O_3$, Acta Mater. 48 (18-19) (2000) 4375-4386, 10/bz82tc.

[102] S. Nufer, et al., Quantitative atomic-scale analysis of interface structures: transmission electron microscopy and local density functional theory, Phys. Rev. Lett. 86 (22) (2001) 5066-5069, 10/cng7hx.

[103] A.M. Dymshits, et al., Thermoelastic properties of chromium oxide $Cr_2O_3$ (eskolaite) at high pressures and temperatures, Phys. Chem. Minerals 43 (6) (2016) 447-458, 10/gtd4ff.

[104] A.H. Hill, A. Harrison, C. Dickinson, W. Zhou, W. Kockelmann, Crystallographic and magnetic studies of mesoporous eskolaite, $Cr_2O_3$, Micropor. Mesopor. Mater. 130 (1-3) (2010) 280-286, 10/cfsbc.

[105] H. Kudielka, Die thermische Ausdehnung der isotypen Mischreihen-Endglieder $Cr_2O_3$ und α-$Al_2O_3$, ermittelt mit einer neuen, lichtstarken-Bohlin-Kammer, Monatsh. Chem. 103 (1972) 72-80, 10/ch22ng.

[106] W.P. Davey, Precision measurements of the crystal structures of $Al_2O_3$, $Fe_2O_3$ and $Cr_2O_3$, Phys. Rev. 21 (1923) 716, 10/d67zvd.

[107] J. Graham, Lattice spacings and colour in the system alumina-chromic oxide, J. Phys. Chem. Solid 17 (1-2) (1960) 18-25, 10/d46fhm.

[108] H. Saalfeld, Strukturuntersuchungen im System $Al_2O_3$-$Cr_2O_3$, Z. Krist. 120 (1-6) (1964) 342-348, 10/gtd4f7.
[109] G.K. Lewis, H.G. Drickamer, Effect of high pressure on the lattice parameters of Cr₂O₃ and α-Fe₂O₃, J. Chem. Phys. 45 (1) (1966) 224-226, 10/dknrws.

[110] Y. Sato, S. Akimoto, Hydrostatic compression of four corundum-type compounds: α-Al₂O₃, V₂O₃, Cr₂O₃, and α-Fe₂O₃, J. Appl. Phys. 50 (8) (1979) 5285, 10/dchpz5.

[111] H. Sawada, Residual electron density study of chromium sesquioxide by crystal structure and scattering factor refinement, Mater. Res. Bull. 29 (3) (1994) 239-245, 10/bsh8vx.

[112] T. Ivanova, M. Surtchev, K. Gesheva, Characterization of CVD chromium oxide thin films, Phys. Stat. Sol. (a) 184 (2) (2001) 507-513, 10/cq6nqf.

[113] P. Dera, B. Lavina, Y. Meng, V.B. Prakapenka, Structural and electronic evolution of Cr₂O₃ on compression to 55 GPa, J. Solid State Chem. 184 (11) (2011) 3040-3049, 10/dr9q8k.

[114] T.J.B. Holland, R. Powell, An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids, J. Metam. Geol. 29 (3) (2011) 333-383, 10/cqwm4g.

[115] S.V. Ovsyannikov, L.S. Dubrovinsky, High-pressure high-temperature synthesis of Cr₂O₃ and Ga₂O₃, High Press. Res. 31 (1) (2011) 23-29, 10/fh8qdp.

[116] T. Larbi, B. Ouni, A. Gantassi, K. Doll, M. Amlouk, T. Manoubi, Structural, optical and vibrational properties of Cr₂O₃ with ferromagnetic and antiferromagnetic order: a combined experimental and density functional theory study, J. Magnetism Magnetic Mater. 444 (2017) 16-22, 10/gtfrgr.

[117] J. Singh, V. Verma, R. Kumar, R. Kumar, Structural, optical and electrical characterization of epitaxial Cr₂O₃ thin film deposited by PLD, Mater. Res. Express 6 (10) (2019) 106406, 10/gtd3xz.

[118] P. Zhao, H. Zhao, J. Yu, H. Zhang, H. Gao, Q. Chen, Crystal structure and properties of Al₂O₃-Cr₂O₃ solid solutions with different Cr₂O₃ contents, Ceram. Int. 44 (2) (2018) 1356-1361, 10/gcr3m2.

[119] A. Rohrbach, J. Hafner, G. Kresse, Ab initio study of the (0001) surfaces of hematite and chromia: influence of strong electronic correlations, Phys. Rev. B 70 (12) (2004) 125426, 10/d94hkt.

[120] C. Chen, H. Yu, S. Zheng, First-principles study of hydrogen diffusion mechanism in Cr₂O₃, Sci. China Technol. Sci. 54 (1) (2011) 88-94, 10/bh3nbb.

[121] J.J. Carey, M. Legesse, M. Nolan, Low valence cation doping of bulk Cr₂O₃: charge compensation and oxygen vacancy formation, J. Phys. Chem. C 120 (34) (2016) 19160-19174, 10/f82g4r.

[122] M.M. Islam, T. Couvant, P. Marcus, B. Diawara, Stress concentration in the bulk Cr₂O₃: effects of temperature and point defects, J. Chem. 2017 (2017) 1-8, 10/gs4v8r.

[123] L.R. Rossi, W.G. Lawrence, Elastic properties of oxide solid solutions: the system Al₂O₃-Cr₂O₃, J. Am. Ceram. Soc. 53 (11) (1970) 604-608, 10/bdjht6.

[124] H. Schober, T. May, B. Dorner, D. Strauch, U. Steigenberger, Y. Morrii, Lattice dynamics of Cr₂O₃, Z. Phys. B 98 (1995) 197-205, 10/fvsgzp.

[125] A. Holzheid, H. St, C. O'Neill, The Cr-Cr₂O₃ oxygen buffer and the free energy of formation of Cr₂O₃ from high-temperature electrochemical measurements, Geom. Cosmochim. Acta 59 (3) (1995) 475-479, 10/dn96k6.

[126] O.H. Krikorian, Thermal expansion of high temperature materials, in: U.S. Atomic Energy Commission, Lawrence Radiation Laboratory, UCRL-6132, 1960, https://doi.org/10.2172/4137098.

[127] J.G. Goedjen, J.H. Stout, Q. Guob, D.A. Shores, Evaluation of stresses in Ni-NiO and Cr-Cr₂O₃ during high temperature oxidation by in situ X-ray diffraction, Mater. Sci. Eng. A 177 (1-2) (1994) 115-124, 10/d3rjm5.

[128] J.J. Barnes, J.G. Goedjen, D.A. Shores, A model for stress generation and relief in oxide-metal systems during a temperature change, Oxid. Met. 32 (5-6) (1989) 449-469, 10/fd7hfp.

[129] B.J. Skinner, S.P. Clark Jr., Thermal expansion, in: Handbook of physical constants 97, Geological Society of America Memoirs, 1966, pp. 75-96, 10/mc4r.